fl...

‘ w
. mg , .
g.“ , a ......
gm“

H.

2 a.
.u‘unhrd

 

 

& .1. ..

1L,

 

 

 

 

 

. . 4%.?

I .il! .7.

 

 

20k“

9

 

LlBRAfifiY
Michigan State
University

 

This is to certify that the
dissertation entitled

HIGH RESOLUTION ANALYSIS OF EXTRAGALACTIC
GLOBULAR CLUSTERS

presented by

Christopher 2 Waters

has been accepted towards fulfillment
of the requirements for the

Ph. D degree in Physics and Astronomy

 

 

 

 

éé 4/"

 

Major Professor’s Signature

Date

 

MSU is an afinnative-action, equal-opportunity employer

— -.-.-.-._.-._ _ L

_ -.-.—.—-.-a--.—

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE
MAR (144‘: 2t. :31 3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:/ClRC/DateDue.indd-p.1

 

 

 

HIGH RESOLUTION ANALYSIS OF EXTRAGALACTIC GLOBULAR
CLUSTERS

By

Christopher Z Waters

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Physics and Astronomy

2007

ABSTRACT

HIGH RESOLUTION ANALYSIS OF EXTRAGALACTIC GLOBULAR
CLUSTERS

By

Christopher Z Waters

Globular clusters are massive compact groups of stars, with masses that range beyond
IOGMQ. Because they are so large, they can remain bound together as they orbit
their host galaxies. They are also very luminous, which ensures that they can be
seen at distances far beyond the point where individual stars are no longer visible.
The combination of these two qualities makes them wonderful test particles to explore
how the dynamical interactions of stars in the cluster change the observed parameters.
The evolution of these clusters has not been very well constrained by observations.
They must lose mass as they orbit, but the exact way that this mass loss changes

their observed properties is not well known.

M87 is a massive galaxy located 16 Mpc away from the Milky Way. This makes
it much farther than other galaxies that have clearly resolved globular clusters. How-
ever as M87 is so massive (M ~ 1012MQ), its globular cluster population is much
more numerous than those of other closer galaxies. Only about 150 globular clusters
have been detected in the Milky Way, whereas M87 should have close to 10000 clus-
ters. This large population allows for any observed relations to be less influenced by

statistical uncertainty.

The core of M87 was imaged with the Advanced Camera for Surveys on the Hubble
Space Telescope as part of a 50 orbit program in 2005 and 2006. During each orbit,
multiple exposures were taken in the infrared F814W and red F606W filters, giving
total exposure times of 73,8003 in F814W and 24,5008 in F606W. These very long

exposures provide some of the deepest data ever taken with HST.

As the data used in this project can resolve the faintest clusters, we can use it to
investigate the luminosity and mass functions of the globular clusters in M87. The
final sample contains 2091 clusters, with a luminosity function that matches well with
previously published results. The mass function generated from these clusters shows
the signature of mass loss from two-body relaxation. Much theoretical work has been
done to investigate this evolution, but since there are few galaxies in which large
numbers of clusters can be observed, these theoretical predictions have been difficult
to test in the past. The change in the mass to light ratio between clusters of different
ages and metallicity is an important complication in the shape of the mass function.
However, by correcting for these changes, this sample shows that the different color
groups of the M87 globular clusters indicate different formation epochs.

These data also provide much higher angular resolution than previously available
for populous extragalactic systems. This resolution ensures that the clusters are
broader than just simple point sources, allowing them to be fit with theoretical models
of the cluster structure. Such fits show that the relations between the cluster structure
and luminosity appear to be universal, as those found for M87 match well with the
Milky Way, the only other complete sample that exists. These structure fits also
show that the probability of the formation of low mass X-ray binaries in a cluster is

influenced by the rate of stellar interactions.

COpyright by
CHRISTOPHER Z WATERS
2007

To the squirrels of Michigan State University

ACKNOWLEDGMENTS

I would like to thank many people for helping with the science of this thesis. My
discussions with Arunav Kundu have been essential in making sure that the analysis
of the images and data have been done correctly. It is also likely that the superking
program would not be in a useful state if not for the comments and suggestions of
Mark Peacock. Tod Lauer provided great insight into how to correctly think about
the formation of images, which has ensured that my fitting code works as well as it
does. I’d also like to thank Ted Baltz for obtaining the telescope time necessary for
this project, which has given me such an excellent dataset to work with.

I’d like to thank my advisor, Steve Zepf, for organizing this project, and for having
the patience to let me work on problems at my own pace. This extends to my thesis
committee, Jack Baldwin, Megan Donahue, Wayne Repko, and Bhanu Mahanti, who
have made my defense of this thesis fairly easy.

I would like to thank the other grad students for many years of bowling, movies,
and board games. Nathan DeLee and Brian Marsteller, in particular, have helped
over the years by talking out various issues with me.

Katie Rabidoux kindly provided the photometric data for the stars in M3, used in
constructing the color magnitude diagram in chapter 1. This is secondary, however,
to her help in reminding me that crab cakes are wonderful. Finally, I’d like to thank
Julie Krugler for suggesting many trips to cofi‘eehouses, without which I never would

have finished writing this thesis.

vi

TABLE OF CONTENTS

List of Tables ........................................................ x
List of Figures ....................................................... xi
List of Symbols ...................................................... xiv
1 Introduction ...................................................... 1
1.1 Globular Clusters ......................... 1
1.2 Summary of Thesis Goals ..................... 3
1.3 Previous Studies of Globular Cluster Systems .......... 4
1.4 Outline ............................... 5

2 Dynamics of Globular Clusters ................................... 9
2.1 King Models ............................ 9
2.2 Tidal Radii ............................. 13
2.2.1 Tidal Radii with Circular Orbits ............... 14

2.2.2 Tidal Radii with Elliptical Orbits ............... 16

2.3 Mass Loss .............................. 16
2.3.1 Dynamical Friction ....................... 17

2.3.2 Stellar Evolution ........................ 18

2.3.3 Gravitational Shocks ...................... 19

2.3.4 Two-Body Relaxation ..................... 21

2.4 Observed Properties of Globular Clusters ............ 24

3 Point Spread Functions ........................................... 27
3.1 Pixel Response Functions and the Effective PSF ........ 29
3.2 PSFs for Extended Objects .................... 31
3.3 Evaluating the PSF for HST ................... 31
3.3.1 TinyTim ............................ 32

3.3.2 Empirical PSFs ......................... 32

3.4 PSF Generation .......................... 34

4 Data Reduction .................................................. 37
4.1 Image Combination ........................ 37
4.1.1 Lauer Fourier Method ..................... 38

4.1.2 Drizzle ............................. 39

4.1.3 Determining Shifts ....................... 42

4.2 Data Summary ........................... 43
4.2.1 WFPC / 2 ............................ 43

4.2.2 ACS ............................... 44

4.3 Image Preparation ......................... 45

vii

4.3.1 Unsharp Masking ........................ 46

4.3.2 Isophote Fitting ........................ 47
4.3.3 Final Method .......................... 47

5 Globular Cluster Luminosity Function ............................ 57
5.1 Object Detection .......................... 57
5.1.1 Instrumental Magnitudes ................... 59
5.2 Completeness Correction ..................... 59
5.3 Photometric Calibration ...................... 60
5.3.1 Aperture Correction ...................... 60
5.3.2 Color Correction and Extinction ............... 62
5.4 Contamination ........................... 65
5.4.1 Noise Objects .......................... 65
5.4.2 Background Galaxies ...................... 66
5.5 Cluster Candidates ......................... 67
5.5.1 Photometric Accuracy ..................... 71
5.6 Luminosity Function ........................ 73
5.6.1 Color Dependence ....................... 73
5.6.2 Radial Dependence ....................... 74
5.6.3 Comparison with Published Results ............. 75
5.7 Blue Tilt .............................. 82
6 Globular Cluster Mass Function .................................. 87
6.1 Initial Mass Function of Globular Clusters ............ 88
6.2 Mass Loss Rates .......................... 89
6.2.1 Fall & Zhang .......................... 90
6.2.2 Baumgardt & Makino ..................... 90
6.2.3 Lamers et a1. .......................... 92
6.3 Mass Function ........................... 93
6.3.1 Comparison to Mass Loss Models ............... 95
6.3.2 Mass to Light Ratio ...................... 96
6.3.3 Comparison to Previous Results ................ 99
6.3.4 Color Dependence of the Mass Loss .............. 100

7 Structural Parameter Fitting ..................................... 104
7.1 Error Analysis ........................... 105
7.1.1 Position Errors ......................... 107
7.1.2 Magnitude and Background Errors .............. 107
7.1.3 Structure Errors ........................ 108
7.2 Fitting Results ........................... 109
8 Globular Cluster Structure Results ............................... 112
8.1 Effect of Cluster Luminosity on Structure ............ 113
8.2 Effect of Distance on Structure .................. 116
8.3 Relations Between Core Parameters ............... 116
8.4 The Fundamental Plane of Globular Clusters .......... 117
8.5 LMXB Probability ......................... 118

viii

9 Conclusions ...................................................... 124

A Detailed Source Code ............................................ 127
A.1 libkingmodel ............................ 127
A.1.1 kingmodel.c ........................... 128

A.1.2 kingimage.c ........................... 131

A.1.3 kingconvolve.c ......................... 132

A.1.4 kingutils.c ............................ 133

A2 Superking .............................. 134
A21 Initialization .......................... 135

A22 Fit Evaluation ......................... 136

A23 Fitting Procedure ....................... 137

B Globular Cluster Catalog ......................................... 140
References ....................................................... 143

ix

LIST OF TABLES

5.1 Source Extractor quantities ......................... 58
5.2 Color correction parameters ......................... 63
5.3 Sample size .................................. 71
5.4 Radial dependence of the GCLF ....................... 75
5.5 GCLF parameters .............................. 81
5.6 Bimodality with cluster luminosity ..................... 82
6.1 Best fitting mass loss fits to ACS data. ................... 94
6.2 Best fitting mass loss fits with variable T. ................. 99
A.1 Search ranges for fitting ........................... 139

1.1

1.2

1.3

2.1

2.2

3.1

3.2

3.3

4.1

4.2

4.3

4.4

4.5

4.6

4.7

4.8

5.1

5.2

5.3

5.4

LIST OF FIGURES

Optical image of galactic globular cluster M3 ............... 6
Globular cluster color magnitude diagram ................. 7
Optical image of Virgo Cluster ....................... 8
King model mass density and projected surface density profiles. ..... 12
Comparison of W0 and c ............................ 24
Comparison of HST PSF to an Airy disk ................... 29
Comparison of TinyTim and empirical PSFS. ................ 33
Comparison between calculated PSF and star. ............... 36
Original raw ACS image ............................ 49
Distortion corrected ACS image. ...................... 50
Final F606W reduced image .......................... 51
Final F814W reduced image .......................... 52
Result of unsharp masking method of galaxy subtraction .......... 53
Result of ELLIPSE method for galaxy subtraction. ............ 54
Image histograms of galaxy subtracted images ................ 55
Result of final galaxy subtraction ....................... 56
Completeness surfaces. ............................ 61
Aperture correction surfaces .......................... 63
Color correction curve relating instrumental to final colors ......... 65
Original F606W UDF image. ........................ 69

5.5

5.6

5.7

5.8

5.9

5.10

5.11

5.12

5.13

6.1

6.2

6.3

6.4

6.5

7.1

7.2

7.3

7.4

8.1

8.2

8.3

8.4

8.5

F606W UDF image with added noise to simulate the characteristics of the

M87 data. ................................. 70
Color magnitude diagrams of detected objects ................ 72
Color magnitude diagram of the final sample of clusters. ......... 77
Final GCLF. ................................. 78
Comparison of the red and blue cluster GCLFs. .............. 79
Radial Dependence of the GCLF ...................... 80
Comparison of the V GCLF to previously published results ........ 84
Comparison of the I GCLF to previously published results ........ 85
Blue tilt .................................... 86
Comparison of M or M ‘0'25 mass loss model to N-body simulations. . . 92
M87 GCMF with best fitting mass loss models. .............. 95
Best fitting mass loss model with variable T ................. 98
Comparison of the new GCMF to those previously published. ...... 101
Comparison of the GCMFs for the red and blue samples ......... 103
Expected errors in position. ......................... 108
Expected errors in magnitude and background level ............. 110
Expected errors in the structural parameters. ............... 111
Fitting residuals. ............................... 111
Parameter correlations with absolute magnitude .............. 115
Parameter trends in averaged bins of luminosity .............. 120
Parameter correlations with projected galactocentric distance. ...... 121
Correlations between core structural parameters .............. 122
LMXB Structure ............................... 123

Images in this dissertation are presented in color.

 

xiii

LIST OF SYMBOLS

rc Core Radius .................................
W Reduced Potential ..............................
T Mass to Light Ratio .............................
rt Tidal Radius .................................
t,— Relaxation time ...............................
W0 Central potential ...............................
c Concentration ................................
Thm Half mass radius ...............................
Rh Half light radius ...............................
#0 Central surface brightness ..........................
p0 Central mass density .............................
trc Core relaxation time .............................
trh Half mass relaxation time ..........................
H Point Spread Function of a telescope. ...................
.@ Pixel Response Function ..........................
Hefl: Effective PSF .................................
MC IMF cutoff mass ...............................
'7 Mass loss decay exponent ..........................
7' Total mass loss from evaporation ......................

xiv

1 1
13
13
22
24
24
24
24
25
26
26
26

27

30
88
89

89

 

CHAPTERI:
INTRODUCTION

 

Globular clusters (GCs) are spherical groupings of many stars that are bound together
by their own self gravity. They appear in every galaxy observers have studied, and
seem to be among the oldest structures in the universe. High resolution images of
galactic globular clusters, such as the Sloan Digital Sky Survey (SDSS) image of
the globular cluster M3 presented in figure 1.1, show an obvious overdensity of stars

compared to the background, with a very dense core of stars blending together.

1.1 GLOBULAR CLUSTERS

The first globular cluster was discovered by Johann Abraham Ihle on August 26, 1665
while he was observing Saturn (Schultz, 1866). This object, now known as M22, was
just one of the 28 galactic globular clusters that appear in Charles Messier’s catalog
of “nebulae,” published roughly a hundred years later. The nature of these objects
was first suggested by William Herschel (Herschel, 1789), who was able to resolve the
individual stars in the cluster. This fact made it clear that globular clusters must be
composed of many stars, and as the stars in the center of the cluster are more tightly
packed, they must be held together by the mutual attractions between the component
stars.

Following their discovery, star counts were used to estimate the shape of the mass
profile, and with the introduction of statistical mechanics in the late nineteenth cen-

tury, globular clusters began to be modeled as a gas of stars (Plummer, 1911). By the

middle of the twentieth century, steady state solutions were found that could satisfac-
torily match the observed density profiles (Hénon, 1961; King, 1966). However, these
solutions indicated that the structure and evolution of these objects are intrinsically
linked, and must be treated simultaneously. This realization has led to recent at-
tempts to explain globular clusters using N-body simulations of the component stars,
which can investigate both the structure and evolution (Baumgardt & Makino, 2003).

Globular clusters are believed to form from the collapse of giant molecular clouds.
This collapse must happen relatively quickly, as the ignition of hydrogen burning in
massive stars is likely sufficient to disperse any remaining gas. Because of this, all
the stars in the cluster can be assumed to have formed at the same time, and as they
all form from gas with the same chemical enrichment, they must all share the same
metal content.

These two facts make globular clusters excellent laboratories to examine the lives
of stars, as the only difference between the individual stars is their masses. The stars
in the cluster therefore fall upon a path that traces the evolutionary history of the
stars. Figure 1.2 illustrates the color magnitude diagram of the galactic globular
cluster M3. The main sequence of stars is clearly visible up to the turnoff around
V N 18.5. This point marks where massive stars are starting to evolve away from the
main sequence, and up the red giant branch. The main sequence lifetimes of stars
at the turnoff provide an estimate of the age of the cluster. These ages indicate that
the globular clusters in the Milky Way are very old, with current ages on the order
of the age of the universe.

Globular clusters are generally more metal poor than the Sun, consistent with their
great ages. Like the metal poor stars of the galactic halo, the Milky Way globular
cluster system is spherically distributed around the center of the galaxy. This fact
was instrumental in one of the great discoveries of the structure of the Milky Way.

Shapley (1921) used the observation that the globular clusters visible from Earth

are not isotropic in the sky, and that they seem to be distributed about a point 8kpc

away. This led to the conclusion that the Earth is not located at the center of the

galaxy.

1.2 SUMMARY OF THESIS GOALS

This thesis examines very deep observations of the central regions of the giant elliptical
galaxy M87. This galaxy is a member of the Virgo cluster of galaxies, and is located
at a distance of 16 Mpc (Macri et al., 1999). Figure 1.3 shows a section of the Virgo
cluster as seen by the SDSS (Stoughton et al., 2002) with the bright object in the
lower left corner being M87. The observations for this project were taken using the
Advanced Camera for Surveys (ACS) aboard the Hubble Space Telesc0pe (HST).
This camera has a field of view of 202'.’ x 202'.’, which covers an area roughly the
size of the bright core of M87 visible in the SDSS image. Although this relatively
small image size does not allow the analysis of the full M87 globular cluster system
(believed to number up to 10000 objects), by taking advantage of very long exposure
time images, the sample of the GCS will be effectively volume limited, with even the
faintest objects measured.

Without a cutoff in luminosity, such a sample can be used to construct a complete
luminosity function for the M87 globular cluster system. This luminosity function can
then be used to constrain the evolution of the clusters. The high angular resolution
that HST provides also allows the clusters to be resolved, which allows their structure
to be examined. Therefore, this data set presents an opportunity to put observational
constraints on the structure and evolution of the cluster properties. As the number
of clusters that comprise the sample is very large (N ~ 2000), such constraints have

the power to be more statistically rigorous than can be provided by the Milky Way.

1.3 PREVIOUS STUDIES OF GLOBULAR CLUSTER SYS-

TEMS

For obvious reasons, the Milky Way GC system was the first to be examined in detail.
The Harris (1996) compilation is the standard catalog of Milky Way clusters, drawing
the best determined values from the published literature. The relations between the
various structural parameters were examined by Djorgovski & Meylan (1994), which
is currently the most reliable set of correlations between the parameters, despite the
low number of clusters. I

Although M31 is an obvious next choice for a target, many factors contribute
to hinder the classification of the GC system. As the galaxy is mostly face on, the
clusters are projected against the galaxy light, which makes them difficult to detect. In
addition, the large angular size of the galaxy means that a complete survey requires
many pointings. The most complete current survey is that of Barmby 85 Huchra
(2001), containing 400 clusters.

Elliptical galaxies are generally considered the best locations to search for GCs as
they have smooth light profiles against which clusters Show up easily. In addition,
elliptical galaxies have much larger numbers of globular clusters per unit luminosity.
M87 has been studied before with the WFPC/ 2 camera on HST (Whitmore et al.,
1995; Kundu et al., 1999; Waters et al., 2006). Recently, it was re-examined as part
of the ACS Virgo Cluster Survey (ACSVCS) (C6té et al., 2004) which studied 100
galaxies in the Virgo cluster. The number of clusters examined per galaxy varies
widely, which limits the conclusions that can be drawn from this sample. Another

issue is that the survey did not have the depth required to detect the faintest clusters.

1.4 OUTLINE

This thesis is organized as follows: Chapter 2 discusses the dynamical processes that
shape globular clusters, and presents the background theory that is used in later
chapters. Chapter 3 discusses the point spread function, and how this affects the
image quality. Chapter 4 outlines the data reduction techniques used. Chapters 5
and 6 cover the globular cluster luminosity and mass functions, respectively, and
discuss their applications to the evolution of the globular cluster system. Chapter
7 presents a new method for measuring the structural parameters of the globular
clusters from high quality data, and chapter 8 presents the results of the measured
parameters. Finally, Chapter 9 summarizes the results of this study, and presents
the final conclusions. Appendix A gives the documentation for the code written for
this thesis that simulates and fits globular cluster profiles and images. Appendix B

discusses the final database of cluster parameters.

 

Figure 1.1 Image of the Milky Way globular cluster M3. This clearly shows that the
core of the cluster contains many stars, but that even at radii much larger than this
core, the cluster contains more stars than the background. This image was taken
from the SDSS image server at httpz/,fwwwwikiskycom.

 

 

 

 

Figure 1.2 Color magnitude diagram of the stars in the globular cluster M3. As the
stars are all of roughly the same age, they trace out the path of stellar evolution.
The main sequence clearly is truncated at V ~ 19, indicating that more massive stars
have evolved, populating the red giant branch. Measurements with V > 17 are taken
from SDSS catalogs (Stoughton et al., 2002). The measurement of the brighter stars
shown in blue were supplied by Katie Rabidoux (private communication), based on
observations taken at the MSU 24" telescope.

 

Figure 1.3 Image of the Virgo cluster of galaxies. MS7, the target galaxy for this
project, is located in the bottom left corner. M86 and M84 are the other two large
elliptical galaxies on the right side of the image along with the stream of galaxies
known as Markarian’s Chain. This image was taken from the SDSS image server at
http:,t'y’www.wikisky.com.

 

CHAPTER 2:
DYNAMICS OF GLOBULAR CLUSTERS

 

2.1 KING MODELS

The standard assumption for modeling the dynamics of globular clusters is that they
are comprised of a large number of identical stars of mass m moving in the potential
created by their own self gravity. We can initially assume that the system is collision-
less, and neglect two-body interactions. The distribution function of the stars can
then be defined to be f(a':', 27, t) 2 0.

Jeans theorem states that the steady-state solution to the collisionless Boltzmann
equation is a function only of integrals of motion in the potential (Binney & Tremaine,
1987). If we assume our system is spherically symmetric, then there are four such
integrals: the total energy E and the three components of the angular momentum

vector 1:. This allows us to write our distribution function as

f = NEE)

For our globular cluster system, the potential that binds the stars together is

generated by the stars themselves. Therefore, from Poisson’s equation
V24) = —47er = —47er/fd327 (2.1)

We can make the assumption that the distribution function is isotropic, so the distri-

bution depends only on the particle energy based on the observations that globular

clusters do not have Significant rotation. Therefore,

ch = f(E)

A first guess for the distribution function for the stars in a globular cluster is that
of the isothermal sphere, in which we assume that the energy of all stars is given by

a Maxwellian distribution:

-—¢I>— 1mv2

_ 2 +
fisothermal(E) = fOe E/a = f06 0 (2'2)
This yields an equation for the density
00 —Q—%1 11202 2 _ (p
pisothermal(q)) = 471/0 f05 0 v d” = p08 E2 (2-3)

Poisson’s equation allows us to rewrite this in terms of the radius from the center of

the cluster instead of the local potential

 

1 d 26“)
—47er — T—QE (7' E)
0.2
pisothermalir) = 277GT2 (2'4)

From this potential, we can define the core radius (called the King radius by some

902
Tc— V 47er0 (2.5)

This radius is the point where the projected density equals roughly one half the central

authors) as

 

value. This radius represents a convenient radius scale for the cluster.
It is clear that this is a poor solution for globular clusters, as the total mass

diverges to infinity when integrated over all radii. We can modify the isothermal

10

 

sphere by noting that globular clusters orbit host galaxies, and as such, experience
an external potential <I>G. This assumption leads to a cutoff in the size of the cluster,
as there will be a point where a test particle feels an equal pull to the cluster and the
galaxy. This point is known as the tidal radius of the cluster (see section 2.2).

We can set the distribution for this model to be a modified “lowered isothermal”

model:

f0 (TE/‘72 — e‘ET/a2 E < ET
ms): ( l (2.6)
0 E 2 ET

We can make the simplification that any particle that “just reaches” the tidal radius
will be stripped from the cluster, so we can set ET = <I>(rt). We can then renormalize

the potential such that <I>(rt) = 0, which gives

<I>+1mv2
f0(e_ 02 —1)E<0

fKing(E) = (2-7)
0 E 2 0
This distribution function yields a density of
”max <1>+ 21mv2
PKing = 47rf0 e_ a 11sz (2.8)
0
where '0me = _%1q_> from the virial theorem. Defining a “reduced potential” W =

—%, allows this to be solved as
0'

mam = A (eW arm/w) — (Ea: (1 + §W)) (2.9)

We can again use Poisson’s equation to find a the density at a given radius
1 rzd—W = —-47rG)017'2 eW erf(\/ W) — (LEW 1+ 2W (2.10)
dr d'r 7r 3

11

Unlike the isothermal sphere, this differential equation must be solved numerically
(see Appendix A.1 for details on fast and accurate evaluation of King models). Figure
2.1 shows how the volume and projected surface density change with radius for a set

of central potentials, W0.

 

 

 

 

 

 

 

 

 

 

 

I I I fir r r I T
0 I: -' ~u-I: -_“-_a-.._ ‘_._ _ ‘
._ .‘.\.\ ~fi 0 _ \
I- .. -\ \\ —I
-1 _ 3‘".- ~ ~ \\ .-
'- ' ‘ A
A x ”5‘1 '
e2 . ' r
\ 8
A V >-
u - \
a“ 1:
__ v
33 ~ w -2 —
V
a . 2
I— ho '-
-4 L- .9.
- . “ _3 __
\ ‘ ‘ I
-5 — ‘_ ‘- II ‘I
_ x '._ - .
-6 . I g l t I ."z ‘ ' _4 . l l I I 1 ii i1
-4 -3 g —1 0 -4 — g -1 0
10810 7‘/7"t) 10810 T/Tt)

Figure 2.1 Mass density and projected surface density for single mass King models
values of W0 = 3 (red), 6 (green), 9 (blue), and 12 (magenta). The highest value of
W0 falls off the fastest in this plot.

Although this model has been used successfully to fit globular cluster surface
brightness profiles, deviations for high quality profiles have led to many attempts
at a more realistic model. The most intuitively reasonable is the multi-mass King
model (Da Costa & Freeman, 1976), which allows the stars in the cluster to have a

distribution instead of just a single value. If we assume energy partition between all

12

2_

the stars, such that midi —

mjajz then the distribution function becomes

<I>+1mv2
Zicif0(e— 02 —1) E<0

0 E20

fMM(E) = (2.11)

where the C, are the fractional masses for each class of stars. Extra care must be
taken when using such models to fit clusters, because although the mass is given by
this distribution function, the fact that the mass to light ratio T may not be the
same for each class of stars means that the observed light distribution may not match
the mass distribution.

Another change to the standard King model that has been proposed is to relax the
requirement of velocity isotropy (Gunn & Griffin, 1979). If the stars have different
masses, then some of them will be scattered onto non radial orbits. Such a model

provides a distribution function of the sort

<I>+lmv2
2(6—22—0 _.)

0 E20

f (E, L) = (2.12)

However, both of these changes to the simple Single mass King model require very
precise star counts to observe any deviations from the predicted projected density of
the standard single mass King model projected density. They also require many more
assumptions about the properties of the stars in the cluster, and because of this, can

be adjusted in many ways to ensure a good fit.

2.2 TIDAL RADII

As globular clusters orbit a host galaxy, their sizes must be constrained in some

fashion. At some radius rt, their density profiles must drop to zero, and all stars

13

 

outside this radius must be considered bound instead to the host galaxy. A simple
estimate of this size can be found by calculating the point at which tidal interactions
with the host galaxy will cause stars to leave the cluster.

The gravitational force keeping a star bound to the cluster is obviously just

 

GM M
Fcluster = ——(2'2—: (2°13)
7'0
Along a line connecting the cluster and the host galaxy, the tidal force will be
QATGMgM
Ftidal = 3 * (2-14)
"G

In this case, setting Ar :—: r0 gives the tidal force felt by a star along this line, and

allows us to estimate the tidal radius as

/ M
rt = r0 3 _2ll/ICG (2.15)

2.2.1 TIDAL RADII WITH CIRCULAR ORBITS

This result neglects the fact that the cluster is not stationary in the galactic potential,
but instead orbits the galaxy. If we take the simple assumption of circular orbits,
then the cluster and host galaxy experience a potential that can be considered time

independent in a frame that rotates about the center of mass with angular velocity

(.02 = G(MG;- MC) (2.16)
’"G

 

14

In such a frame, Jacobi’s integral is constant, and is given by

1
EJ = i112 + CD0 + (DC (217)
1 1
= 579 + -2-r2w2 + <I>G+C (2.18)
1
= -2—7'~2 + <I>eff (2.19)

where (Deff is the effective potential felt by a star:

 

= GMG + GMC + 1720(MG + mg)

2.20
TC TC 2 Ta ( )

(peff

As 1‘2 is always a positive quantity, there exists a range of valid potentials to keep E J
constant, which in turn constrains the region of space in which the stars can occupy
and be part of the cluster. This boundary is the Roche surface for the cluster, and
we can estimate the tidal radius by taking the distance from the center of the cluster

to the Lagrange point between the cluster and the galaxy. This yields a value

/ M
m = 7G 3 _31WCG (2.21)

which is close to the estimate from a static case (Binney & ’Ifemaine, 1987).
Unfortunately, this tidal radius is poorly defined. The surface of zero velocity
is not spherical, and as such, the radius changes at different orientations with the
host galaxy. In addition, stars that pass this radius are not lost with 100% efficiency.
Although they are likely to remain unbound if they pass beyond this surface, they will
still travel along with the cluster for some time. Finally, the most important problem
is due to clusters not orbiting on perfectly circular paths, but rather following elliptical

orbits.

15

 

 

2.2.2 TIDAL RADII WITH ELLIPTICAL ORBITS

A cluster on an elliptical orbit experiences a potential that varies with time. As this
is difficult to calculate (as it also depends on the relaxation time of the cluster), the
standard assumption (King, 1962) is that the tidal radius is set at pericenter, and

stays fixed at that size for the remainder of the orbit. This assumption leads to a

/ M
_ 3 C
Tt — Tc; —————-(3 + e)MG (2.22)

where e is the orbital eccentricity of the cluster. Other models (Innanen et al., 1983)

tidal radius with the form

have modified this result by considering the motion of the cluster through a realistic
mass distribution, and calculating the tidal radius at pericenter. This gives just a
slight change, rt = gr, elliptical-

Recently, Brosche et a1. (1999) have suggested that the ratio of the true observed

tidal radius to the theoretical tidal radius can be parameterized by the form

Tt observed c TA b
—— = 10 (—) (2.23)
rt theoretical 7P

where r A and r p are the apo- and peri-center distances. They provide their preferred
values of b = 0.664 and c = —0.405. This can be constrained from observations, as

the theoretical tidal radius is pr0portional to M 1/ 3 (Baumgardt & Makino, 2003).

2.3 MASS LOSS

Globular clusters lose mass continually over their lifetimes. There are four main
methods by which mass leaves a globular cluster: dynamical friction, stellar evolution,
gravitational shocks, and evaporation due to two-body relaxation. As each of these
mechanisms have differing time scales, the relative contribution from any one method

changes over the cluster’s lifetime.

16

2.3.1 DYNAMICAL FRICTION

Dynamical friction is a slowing felt by a mass as it travels through a population of
other masses. In terms of a globular cluster, dynamical friction serves to degrade the
orbit of the cluster as it passes through the stars of the galaxy.

Following Binney & Tremaine (1987), we can write the force felt by a globular

cluster due to the interactions with the galaxy mass density:

—47r1nAG2M2p(r) 2X _ 2
F dynamical friction = C (erf (X) e X ) (2-24)

2 __
”c

\/7—r

where lnA is the Coulomb logarithm, X = 72%; ~ 1, and a?) = G—Aim. This force

will cause the cluster to lose angular momentum at a rate

dL F 7“ dr
as the cluster speed can be assumed to remain the same. We can solve this differential
equation for the time it takes for all of the angular momentum to be dissipated, i.e.,

r —* 0. Doing this, with the M87 mass distribution presented by Vesperini et a1.
(2003) yields a timescale

3.0905x1013yr 120 MC ‘1 r 2
tdynamical friction ~ In A m E k—pc (2-26)

This very long timescale makes it clear that the effects of dynamical friction are

 

expected to be very small over the lifetime of a globular cluster. In fact, only the
most massive clusters are likely to experience any significant mass loss from this

mechanism.

17

 

2.3.2 STELLAR EVOLUTION

The next source of mass loss in globular clusters comes from the evolution of the stars
that make up the cluster. Taking a simple relation for the main sequence lifetime of

bright, massive stars (where M > M9) (Chernoff & Weinberg, 1990)

—3
tMS ~ 6 x 109 (3%) yr (2.27)

We can see that these lifetimes can be far less than the expected lifetime of the cluster,
especially for the most massive stars. Since all stellar remnants have a smaller mass
than the initial mass of the star, having the highest mass stars fully evolve on very
short time scales (tMS ~ 106) will clearly cause a sudden drop in cluster mass after a
similar time period.

The effect of stellar evolution on the cluster mass can be estimated by taking a
distribution of stars (such as a Kroupa (2001) IMF with NI. ~ 1 x 106), and letting it
evolve with the mass of the star immediately changing to the mass of its remnant after
its main sequence lifetime has elapsed. This model is clearly an oversimplification, as
it ignores mass loss due to stellar winds for high mass stars, and neglects the entire
post—main sequence evolution. Regardless, it will provide an estimate of the relative
importance of stellar evolution on the cluster mass loss.

Based on the main sequence lifetimes listed above, it is clear that the majority
of stellar evolution happens within the first 1 Gyr, after which, the rate should fall
quickly. As the fraction of mass in a cluster due to massive stars is assumed to be
constant (as this is only dependent on the cluster IMF), we can write the mass loss
using the form

dM

_ z _ 2.2
dt VSEM ( 8)

where USE is a time dependent function related to the number of stars leaving the

main sequence. It is clear that USE must have a sudden drop at early times as the

18

most massive stars evolve. This drop should then slow to a small value as lower mass
stars become the “most evolved,” as these stars evolve on much longer timescales.

The fact that the majority of the mass loss due to stellar evolution occurs very
early in the cluster’s lifetime is convenient, as it means we can simply ignore all
consideration of this mass loss by simply scaling all cluster initial masses to their
value after stellar evolution has subsided, and then ignore the effect on the subsequent
evolution. This is reasonable as the fact that the mass loss equation above has the
solution

M(t) = Moefci V8E<t’)dt’ (2.29)

Mass loss with this form will only change the normalization of the mass function,
although the shape will remain the same.

Finally, it is worth noting that the sudden drOp in mass at early times is likely
to be sufficient to disrupt clusters that have low binding energy. This suggests that
clusters that form with very low initial concentrations, or alternatively, with a dis-
proportionately large number of massive stars relative to the total, are unlikely to

survive past the first 1 Gyr, and will be absent from later surveys.

2.3.3 GRAVITATIONAL SHOCKS

Gravitational shocks occur as the cluster passes by a large but finite mass distribution,
such as the disk of a spiral galaxy or a giant molecular cloud in a galaxy halo. The
passage of the cluster by these objects creates gravitational tidal forces that on average
transfer energy to the stars in the cluster, making them more likely to escape the
cluster.

Given the tidal acceleration acting on a point a: in the cluster due to an object at

distance R from the cluster,

2$GM
“tidal = __R—3£ (2-30)

19‘

 

we can substitute R2 = b2 + v2t2, where b is the distance of closest approach for the
cluster and the perturbing mass. Using the impulse approximation for this accelera-

tion yields the expected change in velocity for the test star:

 

2GMPSU
A = 2.31
v 522) ( )
We can then find the average change in energy per unit mass
_ 1 2 _ 1 2 20M}: 2
AE — 2A2) —— 2a: ( b2v ) (2.32)

and, by integrating over the entire cluster, determine the total change in energy in

the cluster stars (Binney & Tremaine, 1987)

M
AEtotal = T391772; (

 

(2.33)

2GMp)
b2?)

The timescale for this process to disrupt a cluster is of the order for this change

in energy to exceed the binding energy of the cluster

tSH N V_1 Ebinding

encounters AEtotal (2'34)

For shocks from the galaxy disk (q.v. Spitzer, 1987; Ostriker et al., 1972; Fall & Zhang,

2001)
3 Gil/[01302]?
tSH = —— (2.35)
20 rgggn
where PC is the orbital period of the cluster around the galaxy, ”Z is the cluster
velocity as it passes through the plane of the disk, and gm is the maximum gravita-

tional acceleration experienced by the cluster. For Milky Way clusters, this gives a

timescale on the order of tSH ~ 6 x 109yr. With this dissolution timescale, we can

20

—--—.4~’ ‘

write the mass loss rate as (Fall 81: Zhang, 2001)

7 A
VSH = —iS—— (2.36)
3 tSH

where the factor of gfl accounts for corrections to the simple treatment of the deriva-
tion using the impulse approximation, and It S is the constant relating how much mass

is lost for a given change in cluster energy

M_K§
M "‘ SE
This directly gives the relation
dM
:lt— — z/SHM (2.37)

As before with stellar evolution, this form has an exponential solution, which again
yields evolution that does not change the shape of the mass function, but only shifts

the normalization over time.

2.3.4 TWO-BODY RELAXATION

The final mass loss process and generally the most important is evaporation from
the cluster due to two body relaxation. In the dense stellar environments of globular
clusters, there will be many interactions between the stars, which serves to turn
the velocity distribution into a Maxwellian. However, the high velocity tail of this
distribution will ensure that some stars in the cluster will have velocities that will
move them beyond the tidal radius of the cluster, at which point they will be lost to
the host galaxy.

This mass loss must then have a timescale related to the time needed for the

cluster to achieve a Maxwellian velocity distribution. This is the definition of the

21

 

relaxation time, which we can write as

 

tr = 2 (2.38)

where v2 is the is kinetic energy per mass for a star, and ((Av)l2l) is the average change
of this energy per unit time. These quantities can be found by solving the Fokker-
Planck approximation for the cluster (Spitzer, 1987; Binney & Tremaine, 1987 ), yield-

ing a relaxation time:
.03

t = 0.065—————
" nm2G'2 lnA

(2.39)

Generally for a globular cluster, the important timescale is the relaxation time within
the half mass radius. In this case n = %%€—‘17 = 8%,??? and assuming the velocities

are circular
M1/27‘Z/2

t —_
r or mGl/2lnA

(2.40)

From the virial theorem, we know that the kinetic energy of the stars in the
cluster is equal to half the cluster potential energy. This leads to the conclusion that
the escape velocity is only twice the average velocity of a star in the cluster. If we
determine how many stars have such velocities over the course of a relaxation time,

we can also find the expected mass loss rate:
uev = — (2.41)
where fie is the probability that a star with escape velocity is able to reach the edge of

the cluster before being scattered back to a lower velocity. Following Spitzer (1987)

00 00
5e = @- (v)v2dv=i§/ e_x23:2d:13 (2.42)
”f 2v.n 7rl/ 2.45

which gives an estimate of {6 ~ 7.4 x 10—3, within an order of magnitude of the value

22

of fie = 0.045 given by more detailed evaluation with realistic mass distributions for
the stars (Hénon, 1961). This gives a mass loss rate of (Fall & Zhang, 2001)
7.2556mGI/21nA

Vev = (2.43)
1 2 3/2
M / Th

 

Substituting in the mean cluster density ,5, which we assume to be constant over the

lifetime of the cluster, we find

Vev = kg. (G',5)1/2 mln AM‘1 (2.44)
and a mass loss equation
dM
= kg. (Gp)1/2 mlnA (2.46)

As this mass loss is independent of mass, the solution is a linear decay, where

M(t) = M0 — uevt (2.47)

This is a different form for the mass loss than for stellar evolution and shocks. It
gives an amount of mass lost per unit time that is independent of the cluster mass
and constant in time. Because of this, small clusters will be completely destroyed
due to evaporation long before the heaviest clusters have lost even a fraction of their
mass. Such mass loss has the effect of quickly depleting the mass function at low
mass. After time t, we would expect that all clusters with initial mass M0 < nevt to
have been fully disrupted. This does not mean, however, that no clusters lower than
this will be found. These new low mass clusters will be the remnants of more massive

clusters that have also lost mass through evaporation.

23

 

 

2.4 OBSERVED PROPERTIES OF GLOBULAR CLUSTERS

The structure of globular clusters is defined in terms of the central potential W0 This
parameter is not however, how the structure is generally defined based on observa-
tions. Instead, the King model concentration is used in its place, and is defined based
on the two main lengths for the cluster

c: loglo (£3) (2.48)

C

Figure 2.2 shows that c can be used as a replacement for W0.

 

 

 

 

3 l l I l l
2.5- 4

Cl

.2 2 — -

E

1%

31.5— ~

0

O

2?

'22 1 " ‘
0.5— _
0 l l J l l 1_

0 2 4 v9 8 10 12
0

Figure 2.2 Comparison of the central potential W0 with the observable concentration
c.

The half mass radius Thm is another commonly measured quantity is defined as
the radius that contains half of the total mass of the cluster. As the mass is generally
an inferred parameter based on the surface brightness, the half light radius, Rh is

often used in its place. This radius contains one half of the total flux of the cluster,

24

and can be used as a suitable substitute for the half mass radius. Both of these radii
are larger than the core radii for most clusters.

A variety of photometric parameters can also be defined to quantify the flux
normalization. The most common of these is the central surface brightness, 440 For
most observations, this is difficult to measure accurately, as it requires a fit to the
surface brightness profile, which is then extrapolated to zero radius. Therefore, it is
often replaced with the average surface brightness within the half light radius, which

is more easily measured. The equation for this is simply
(4),, = V + 2.510g10 (27.3%,) (2.49)

The central surface brightness can be used to find the central luminosity and mass
densities. Converting the central surface brightness in magnitudes into a luminosity

surface density (in L—JEEZ) is reasonably straight forward
log10(20) = 0.4 (26.362 — W(0)) (2.50)

which can be used to determine the luminosity density in the center of the cluster:

. 20
.70:—

2.51
Tcp ( )

In this case, p is a function of the King model that determines how much of the surface
luminosity is contributed from the core, and how much arises from the projection.

This can be calculated from the King model by taking

p~/O Z(r)dA//0 p(r)dV (2.52)

Converting this central luminosity density into a mass density is just a simple multi-

25

plication by the mass to light ratio:

p0 = T30 (2-53)

The relaxation time is the main timescale for the evolution of the cluster. At the

core of the cluster, the relaxation time is defined as (equation 2.3.4)

k _ 1/2
_ 7 1 3
tTC — 1.491X1O E—A-<m*> p0 Tc (2.54)
—1
1 2 4
= 2.5013 x 107p0/ r2 (In (%M)) yr (2.55)

This timescale often overestimates the evolution, so the half mass relaxation time is

generally used instead. This time is

1
tr}, = 8.933 x 105m(m)—1M1/2r2/2 (2.56)
—-1
= 2.6799 x106M1/2rZ/2(1n(3£0-M)) yr (2.57)

Finally, the metallicity of globular clusters can be estimated reasonably well from
the photometric color. This relationship is defined based on the known metallicities
and colors from the Milky Way clusters, and is usually written in the form (Kundu
et al., 1999):

[Fe/H] = —5.89 + 4.72 (V — I) (2.58)

As globular clusters have colors around V -— I ~ 1.0, they have metallicities around
[Fe/ H] ~ —1.2, illustrating that globular clusters are more metal poor than the sun,

which is to be expected based on their ages.

26

 

CHAPTER 3:
POINT SPREAD FUNCTIONS

 

At the distance of M87, the projected sizes of globular clusters are similar to the
resolving limit of HST. Because of this, we have to account for the effects of the point
spread function, If. The point spread function (hereafter PSF) defines how the light
of a point source (such as a star at great distance) is spread over the focal plane of
the telescope by diffraction.

In the simplest model, we can consider a telescope as simply a circular aperture
separated from the focal plane by a distance m. If we only allow monochromatic
light of wavelength A0 to pass through the aperture, then we can model the light as it

passes through the aperture as a series of spherical wavefronts. This yields an image

E = // §2ei(kr'wt)dA (3.1)
A

perture 7'

amplitude of

For the circular aperture, we can break the integral into strips of dA = :cdy where
at = 2m. Since the difference in amplitude due to diffraction is related to the
differences in path length between all contributions, we can rewrite the equation of
the spherical wave at a given time to explicitly contain this difference in path length

k1“ — wt = yk sin 0. Performing the integration, and making further substitutions

27

v = if and '7 == kRsinO

 

2 1 _
E = 239R / mm
0 -1
27TE0R2 J1(’7)
7'0 '7

 

(3.2)

where J1(:2:) is the first Bessel function. The final image intensity is the square of

this, so

_ J1(kRsin6’) 2
I _ IO( kRsinB ) (33)

which is the definition of the Airy disk.

We can compare this simple model to the known properties of HST by noting that
the core size of the PSF for visible light is roughly (1’05. The first minimum of the
Airy disk occurs when 2R Sine = 1.22A0. Setting A0 = 555nm and R = 1.2m, we
find that this gives a size 6 ~ 07058, showing that the HST PSF is dominated by the
diffraction effects of the telescope aperture. Figure 3.1 shows a comparison between
highly sampled HST PSFS for the F814W filter, and an Airy disk generated at the
peak wavelength of this filter A = 814nm.

The extra structure that is visible in the real HST PSF shows that real PSFS are
composed of more than a single simple diffraction pattern. It is clear from inspection
of the formula above that the effect of diffraction is a Fourier transform of the aperture.
The optical path of HST contains many more elements than a simple circular opening,
such as support structures, and the “spider” that holds the secondary mirror. The

addition of these objects creates the added complexity that is seen in the final PSF.

28

 

Figure 3.1 Comparison of the Airy disk for a circular aperture comparable in size to
the HST to a highly sampled real HST PSF. This PSF is generated for the F814W
filter with the Airy disk created for the peak wavelength of that filter. The first zero
of the Airy disk clearly correlates to the Size of the central core of the true PSF.

3.1 PIXEL RESPONSE FUNCTIONS AND THE EFFECTIVE
PSF

A further complication in modeling the PSF arises from the fact that the focal plane
of the telescope is not a perfect imaging plane, but rather an array of sensors that
make up the CCD imager. For an ideal CCD, the measured image of a point source

would be equivalent to

NM) = 11093069111072) (3-4)

where III is the Shah function, a regular grid of impulse functions. Due to imperfec-
tions in the production of CCDs, pixels do not fully sample light equally well when
that light is centered differently. If we define the pixel response function g as the

sensitivity function of a given pixel over its surface, then the actual detected image

29

will be

KM) = H($,y)®9?($,y)®111(i,3') (3-5)

The pixel response represents the changes in the detection efficiency of a given pixel
when light is incident on different portions of the pixel. One of the main sources of this
is due to the finite thickness of the detector material, and the fact that that material
does not absorb photons with 100% efficiency. This leads to internal reflections off
the back surface of the detector, which in turn causes the light to be scattered into
neighboring pixels. Assuming this scattering is uniform, a point source centered near
a pixel edge will have fewer photons detected in the incident pixel compared with a
source at the pixel center. The individual pixels in the detector are also not perfectly
electrically isolated from each other. This can allow captured photons (now present
as a charge in the detector) to bleed into neighboring pixels. This adds a further
contribution to %.

Since we can never truly observe the real PSF, but only the convolution of the PSF

and 32’, we can define a new function Hefl‘, the effective PSF as

Heff = H (8) .9? (3.6)

By defining this, we now have a function that works as the standard PSF would in
a continuous focal plane. The effective PSF must be smoother and broader than the
instrumental PSF, as it convolved with .%7, which has a width similar to the size of a
pixel. However, by switching to Hefl‘, we no longer need to worry about integrating

over the surface of a pixel, as that integration is incorporated already.

30

3.2 PSFS FOR EXTENDED OBJECTS

We have defined the PSF so far in terms of the effect of diffraction on a point object.
For extended objects, the relations are largely the same. For a point source, the pixel
phase (where the center of the PSF falls in the pixel) changes the flux significantly.
This is not necessarily the case for extended objects, as the fact that the illumination
covers the entire pixel reduces the effect of the pixel phase. If we assume that the
pixel is illuminated by roughly the same flux across its entire surface, then the PSF
for an extended object can be evaluated by simply evaluating the Heff at the center
of the pixel. As the light distribution becomes less flat across the pixel, the PSF shifts
to the location of the peak of the flux. However, this will smoothly shift to the point

like case if the distribution becomes significantly peaked.

3.3 EVALUATING THE PSF FOR HST

Since we need the PSFS for HST to accurately model the detected GCs, we can take
advantage of the fact that being in Space makes the PSF generally stable with time.
For ground based observatories, motions of the atmosphere can significantly change
the width of the PSF. These motions create a PSF that is much larger than the
aperture diffraction pattern. For ground based observations, the PSF can then be

assumed to be of the form
_ 2 2
IIground = Heff ‘3’ e (a: +3; ”208mm (3-7)

where Oatm is of the order of 1’.’, almost certainly larger than Heg. By being in space,
HST does not have this added convolution, and so depends only on the telescope
itself, which leads to consistent modeling of the PSF. Creating such a model for the
PSF is essential for HST, as many pointings do not have sufficient stars to create one

directly from the observation.

31

3.3.1 TINYTIM

TinyTim (Krist, 1993) is a program that is designed to model the PSF for HST. It
takes a theoretical model of the telescope aperture and obstructions, and uses Fourier
Transforms to directly estimate the Shape of the PSF. It incorporates the expected
object spectrum and the filter response to determine the relative contributions for all
wavelengths. The pixel response of the detector is added as a convolution at the end
of the evaluation.

This program works well for WFPC / 2, and is the standard method for generating
PSFS for this instrument. However, the ACS detector is far off the central axis of
the telescope, and as such, has serious geometric distortion. The standard correction
for this distortion is to use Drizzle (see section 4.1.2 for details) to correct the ACS
frames. As the parameters that govern this procedure can change, the default Tiny-
Tim distortion correction does not in general provide a PSF that accurately represents

what is actually observed by the detector.

3.3.2 EMPIRICAL PSFS

To resolve this problem of poorly modeled PSFS, Anderson & King (2006) took the
tactic of empirically measuring the PSFS. To do this, they used 126 orbits of observa-
tions with HST to image the same field in the globular cluster NGC 6397. This field
contained roughly 4000 stars, which were imaged at different rotations and shifts, to
ensure that the stars do not fall on identical locations on every image. A model of the
PSF was then created, and fit to the stars, yielding a position and flux for each star
on each frame. These positions are then used to construct a model of the distortion
for the detector. Finally, the true pixel values were used to update the model of the
PSF. This process was iterated until a final solution was found. This iterative process
is important, as asymmetries in the PSF can alter the measured centroid for the star,

which will in turn yield a worse distortion model (Anderson & King, 2003).

32

The final result of this iterative process was an accurate distortion model, and
a highly sampled Hefl'. As there are 4000 stars on 126 images, each of which is
roughly 10 pixels in radius, it is clear that this data provides millions of samples of
the effective PSF. To eliminate changes in the PSF across the surface of the detector,
the models were calibrated over different areas so that any changes in the PSF beyond
the standard distortion corrections would be accounted for in the model.

This iterative process was repeated for a variety of filters, and then tabulated into
reference images that contain highly sampled PSFS at a variety of positions across the
image. By interpolating between the different PSFS, we can generate a high resolution
PSF for any position on the detector, and then by interpolating that highly sampled
PSF, we can make one for a given pixel phase. These PSFS are designed for the raw

flat fielded and distorted “FLT” frames. This choice was made as this image is the

most photometrically accurate, as it has had the least processing or resampling.

 

Figure 3.2 Comparison of TinyTim PSF to the Anderson 85 King (2006) empirical
ACS PSF after applying the distortion correction. Although the general shapes are
similar, the cores clearly differ, with the TinyTim PSF being more oblong than the
empirical PSF.

33

3.4 PSF GENERATION

To generate PSFS for use in modeling the globular clusters, the individual data FLT
frames were zeroed out, and the empirical PSFS were placed at the positions of globular
clusters in the original data. The shifts between frames were included in this, such
that the PSFS on each frame have slightly different pixel phases. This was done with
the goal that once the frames were distortion corrected and aligned, all PSFS for a
given position will also align. These PSF frames are then combined in the same way
as the data, ensuring that the final PSFS are the most accurate representation of a
point source in the data.

Quantifying the error in the final PSFS is very difficult, as there are few bright
stars in the field. As M87 lies out of the plane of the Milky Way, this lack of stars is
not surprising. Complicating the effect is the fact that one of the two obvious stars
is saturated at the core, making it effectively useless in characterizing the quality of
the PSFS we generate. The one remaining star is located at the edge of the bottom
chip. This location is likely to have the largest distortion effects, and so is likely to
be one of the most difficult locations to generate an accurate PSF.

To classify the error, a PSF was generated for the location of the star ((a:,y) =
(3995, 340)), and the star itself was extracted from the background subtracted double
resolution image (see chapter 4.2.2). The centroids were found for both the star and
the PSF, and radial profiles were generated. AS we are only mostly certain that this
star is not saturated, the choice was made to flux calibrate the PSF by forcing the
profile at 5 pixels to match. This choice excludes the central peak, and calibrates
based on the light in the second maximum. Once this was done, the star and flux
calibrated PSF were interpolated to a common grid, and the percent error between

the two was calculated as
_ Star — PSF

E Star

(3.8)

34

Only the central 20 x 20 pixels were checked, as beyond this point, the star image
has dropped to less than 0.1% of the peak, and noise from the galaxy background
begins to be the main source of light. Inside this box, the error is well described by
a Gaussian with (E) ~ 2% and 0E ~ 6%. The error in the central pixel is 5% which
suggests that this star is not in fact saturated.

This result is reasonable, as it is consistent with the nominal uncertainty in the
empirical PSF. Anderson 85 King (2006) provide an error estimate of 5% for their
empirical PSFS, in the case of no additional corrections beyond the tabulated grid.
As their method of correction is based on doing a Similar error analysis to this, but
for multiple stars across the image, we are unable to apply their method to our data.
These errors are used to construct a perturbation to the standard PSFS, which reduces
the error by about a factor of two.

In addition to the simple errors, we can check that the PSF shape is correct.
Although the star and PSF have nearly identical “bumps and wiggles,” the PSF appears
to have a slightly broader core than the star. Approximating the core of both with
a Gaussian Shows that the PSF is indeed ~ 0.15 pixel broader in both directions
(of;tar = 1.64, 0;” = 1.80, a?” = 1.71, 0531? = 1.85). This error likely arises
from some error in the placement of the PSFS onto the individual frames. In that
coordinate system, this error transforms to an interframe scatter of about 0.08 pixel,
a value which would not be surprising. One possible source of this scatter is from
errors in the geometric distortion, which seems to account for the majority of this

scatter (0distortion ~ 0.05 pixels, Meurer (2002)).

35

 

 

Figure 3.3 Comparison of final drizzled ACS PSF generated from the Anderson &
King (2006) empirical PSFS to the star on the final drizzled data image. The PSF is
made for the position of the star to ensure the distortion calculations match and that
the different response across the detector are removed.

36

 

CHAPTER 4:
DATA REDUCTION

 

The data used for this project comes entirely from the Hubble Space Telescope. AS
discussed in Chapter 3, the PSF for ground based telescopes is much too wide for
extragalactic globular clusters to be imaged accurately. Since we wish to measure
the structure of the clusters, this added broadening will wash out the cluster light,
preventing any structure from being visible. In addition, the broadening makes iden-
tification of faint objects difficult, which would limit the depth to which we can probe.
Very high signal to noise data can be created with HST, which ensures that these

faint objects will be well detected.

4.1 IMAGE COMBINATION

As HST is limited in the length of any single image by its orbit time, multiple short
exposures must be taken to get the long total exposure times needed to achieve the
required high Signal to noise. Using multiple exposures of the same image also allows
for problems on the image to be removed. Since the HST detectors are in space,
they are much more susceptible to interference from cosmic rays, which can interact
with the detector and show up as erroneous bright objects. However, as these cosmic
rays are very unlikely to occur on the same pixels in all images, they can be easily
removed by comparing the multiple frames. In addition, the effects of cosmic rays
build up over the exposure, so a single long exposure will have more contamination

from cosmic rays than a short exposure.

37

The HST detectors also have a number of bad pixels and bad columns, which
do not correctly measure the light incident on them. The main cameras are also
comprised of multiple detectors, which have gaps between them. By dithering the
telescope pointing, the location of the scene on the detectors also changes, and both
of these gaps in the data can be repaired using pixels from other images. The final
benefit of multiple dithered images is that the data can be combined to yield higher
resolution combined images than any of the original images. The creation of this
higher resolution image requires that the dither include non-integer shifts. As we
assume that the pixels of the detector sample the light distribution in some regular
fashion, these fractional shifts sample the light in the intervening Space. The simplest
case is one in which four images are combined with relative shifts (0,0), (0,0.5),
(0.5,0), (0.5,0.5). Since these images sample the light regularly on a grid with twice
the original image resolution, a new image can be created by interlacing the pixels
together.

It is important to note that the construction of new higher resolution images has
little effect on the angular size of the PSF. This new image in no way deconvolves
the scene from the PSF, it merely samples the data better. This means that objects
smaller than the PSF will remain unresolved in the new data. In addition, the very
act of co-adding data tends to introduce further blurring, depending on the technique

and sampling pattern used.

4.1.1 LAUER FOURIER METHOD

One method of recovering resolution from multiple frames is presented by Lauer
(1999), using image combination in the Fourier domain to create a Nyquist sampled
image. Briefly, this method expands each data image by interleaving blank pixels to
pad the data to the final resolution. The Fourier transform of these images is taken,

and multiplied by a phase factor to incorporate the interframe shifts. These transform

38

images each contain some information about the image at frequencies higher than the
original sampling frequency. By summing these transform images the contributions
from each dithered image can be incorporated to provide the best estimate for this
high frequency information. The final combined image is then created by taking the
inverse Fourier transform.

In the simple interlacing case as considered above, this method yields identical
results. However, it can also yield reasonable estimates even when the sampling isn’t
ideal. Another benefit is that overconstrained sets of images can be combined in a
least squares method to provide an image that deals with noise on the individual
frames.

Although this method can accurately recreate the underlying scene, it has fairly
strict data requirements. First, the input frames must be cleaned of all cosmic rays
and errors, and must have any geometric distortion due to the instrument optics
corrected. These requirements make it immediately difficult for dealing with ACS
images, which have large distortion effects. Secondly, the input frames should fully
sample the pixel phase space, otherwise, they will create biased images, that depend
on the individual frames unevenly. These images may also contain aliases from Fourier
domain “satellites,” which can arise if the final Fourier transform image does not
properly taper to zero at high frequency. Such satellites will create a blurring of the
high frequency information, limiting the final resolution that can be created. Finally,
there is no packaged form for this method, which prevents widespread adoption,

despite its mathematical elegance.

4.1.2 DRIZZLE

Another method commonly used to combine data is the Multidrizzle package in Pyraf
(Koekemoer et al., 2002). Although less mathematically rigorous, it has much looser

data requirements and is mostly automated. This process not only combines the

39

images, but also cleans image defects and corrects for the geometric distortion.

The Drizzle algorithm operates on the input images by first shrinking the pixels
by a scale factor called the “pixfrac.” This step holds the centers of the pixels fixed,
but decreases the area, leaving gaps between the pixels. These shrunken pixels are
then transformed via a geometric transformation onto the final image grid. The small
pixels are then “drizzled” onto the final image by allowing the pixel values to be added
to the output pixels in proportion to the area covered.

This addition to the output grid can be altered by changing the kernel used for
the addition, with a “square” kernel simply adding the fraction based on the area of
overlap, and a “Gaussian” kernel weighting the contribution from the center of the
input pixel more than the edges. An important way to estimate the image is to use
the “turbo” drizzle, which contributes all of the input pixel’s light only on a point in
the center.

After all the individual frames have been drizzled to the output image, a rescaling
is done by normalizing against the weight in each pixel. Since all pixels are not
guaranteed to have the same input area contributed to them, this rescaling is essential
to ensure that the counts measured on all pixels have the same basis. If few images
are used, this rescaling can add noise to the final image due to non-uniform weights.

Multidrizzle first reads the flags for each pixel in the raw flat fielded “FLT” im-
ages, and creates a set of masks for known bad pixels. The FLT images are then
projected onto a common coordinate system via a “turbo” drizzle. The WCS from
each image, plus an external correction supplied from a shiftfile, is used to create the
transformation from the “FLT” to the “single_sci.” This step corrects the frames for
geometric distortion, and removes the background ofl'sets between images, ensuring
the various frames have the same geometric and photometric calibration.

At this stage, the many “single_sci” images represent multiple realizations of the

true scene. However, cosmic rays and unknown image defects create deviations from

40

this true scene. We can use the fact that we have many images to remove these
defects. A median image is created from the “single_sci” images, with discrepant
pixels clipped out. This median image is “blotted” back to the FLT frame for each
image, to create “BLT” images, which are an expectation of what each raw FLT frame
should look like.

Cosmic ray rejection is accomplished by comparing the FLT to the BLT images.
To avoid improperly masking the bright peaks of real objects (which are most likely
to have the largest scatter between individual images), the masking is weighted by
a gradient image, in which each gradient pixel is set equal to the largest deviation
between the neighbors in the FLT image. A pixel is then flagged as a cosmic ray if

the difl'erence between the FLT and the BLT exceeds:
S
|FLT — BLTI > scale - VFLT + N - ”noise (4.1)

Errors in the image alignment and sky level can easily lead to improper masking
of objects in the CR phase. If the Shifts are incorrect, then the objects will be clipped
from the median image, and hence will likely be flagged as cosmic rays. Errors in
the sky level between frames can do the same, by limiting the number of images that
truly contribute to the median image. Such an error can allow cosmic rays to slip
through by being improperly excluded from the median image.

Once the cosmic rays are found, their locations are added to the static mask
created earlier. This creates a final mask used for the final drizzle. This drizzle uses
a more accurate “square” kernel with the pixfrac set to 0.7 to minimize the size of the
final PSF. Each FLT is then drizzled onto the final “DRZ” image, with all bad pixels
removed and the sky levels matched. The final DRZ image has an exposure time
equal to the sum of the exposure times of the component FLTS. A weight file is also
created that stores the exposure time sum that contributed to each pixel of the DRZ.

The weight varies between pixels based on the effects of the distortion correction and

41

the number of frames contaminated by cosmic rays at that position. This weight is
essential to gauge how well the Shifts are calculated. In the case of poorly determined
Shifts, the locations of objects in the DRZ image are matched by lower values in the
weight image, as the objects peaks will have been flagged as cosmic rays on the frames

that have the worst shifts.

4.1.3 DETERMINING SHIFTS

Accurate interframe shifts are essential to ensure the highest quality final images. AS
bad pixels are masked by the Drizzle algorithm, based on flagging statistical outliers,
having inconsistent shifts will skew these statistics.

The first step in measuring shifts is projecting all of the data frames onto dis-
tortion corrected frames. This projection mosaics the multiple chips together, and
arranges the frames onto the same WCS. This can easily be accomplished by stopping
Multidrizzle after the creation of the single_sci images.

With the various frames drizzled to what should be a common frame, the errors
between the frame WCS and the true sky WCS manifest as shifts in object coordinates
between the corrected images. A first guess at the transformations needed to correct
each image is created by manually identifying a pair of objects on each image. The
images then have catalogs of objects created with Source Extractor (Bertin 83 Arnouts,
1996). The brightest objects (500 in this case) are extracted from these catalogs. The
initial transformation is used to find matching objects, which are then used to refine
the transformation by solving the least squares problem for each frame 3' at each

matched point z':

Xreferencez’ — XOj = Sj (COS 63'ij + Sin ajyjz’)

Yreferencei —Y0j = Sj (COS9j1/jz' *Sinajwjz')

42

As the transformation is improved, more matches will be found, so this is iterated
until the transformation for each frame converges.

With accurate shifts calculated, we next need to calculate correct sky background
levels. These values are calculated by the calibration process based on the distribution
of image pixel values. Large features (such as the galaxy itself in our data) will Skew
this distribution, and tend to overestimate the background level. The individual FLT
images are drizzled into single_sci frames again, incorporating the new shifts. Empty
regions of the images are found, and then the median in a box 100 x 100 pixels is
taken. Each image then has the median deviation from all boxes calculated, after
subtracting the minimum from all images. This median is written to the FLT image

header, to be subtracted during the final run of multidrizzle.

4.2 DATA SUMMARY

The data used for this thesis come from two sets of many orbit observations of the
galaxy M87. The data were taken as part of an initial microlensing survey with the
WFPC/ 2, and again in a followup survey with the ACS. In order to measure these
microlensing events, they must be monitored over their light curves. This requires
a time series of data taken fairly regularly. For our purposes of creating very deep
images, we can ignore any slight changes due to this microlensing, and simply combine

all the images together into very deep exposures.

4.2.1 WFPC/2

The data for the original microlensing survey were combined by Tod Lauer, using
his optimal Fourier method. Superimages for each of the four detectors on WFPC / 2
were created. These data were taken in two filters, the F606W (a V filter, with peak
wavelength at A = 606nm) and F814W (an I filter, peak wavelength A = 814nm),

with total exposure time W = 116003 and t I = 30160s.

43

These data were used to design and test the methods for detection on such deep
images. A very deep luminosity function was created for the globular clusters (Waters
et al., 2006). However, difficulties in constraining the PSF prevented reliable structure
fitting. Conveniently, the followup ACS data were taken at this time, and the project

switched to use the new data.

4.2.2 ACS

The followup data came from a 50 orbit series of observations with the ACS. The
images are of the core of the giant elliptical galaxy M87, extending out to a projected
radius of 8 kpc. The data were taken over the course of a three month search for mi-
crolensing events, which need multiple exposures to look for the changes in brightness.
This arrangement yields data that can be combined into Single very deep exposures.

The same two filters were used for this project as for the WFPC / 2 survey: F606W
and F814W. On each observing day, four exposures in F814W were taken with slight
pointing offsets to provide for full image sampling every day. This was done as F814W
was the primary filter used for the microlensing search. These exposures are matched
by a single exposure in F606W, which are dithered over the different days, providing
a full sampling of the image plane over the entire set of observations.

In all, 49 F606W and 205 F814W images were combined to yield final images with
exposure times of W = 245008 and t I = 738008, making these some of the deepest
images ever taken with HST. In addition to these exposures, 8 exposures in F606W
and. 13 in F814W were taken but excluded due to a loss of the teleSCOpe pointing. The
main images were combined to a resolution of 0(.’045 pixel-1, the nominal resolution
of ACS. These images therefore have the highest signal to noise possible. A second set
of images were combined at twice this resolution, 0’.’025 pixel‘l, for use in modeling
the clusters. These higher resolution images are useful as they provide a better view

of the cluster structure.

44

4.3 IMAGE PREPARATION

Given that the final combined exposure times make the final images among the deep-
est ever taken with HST, we have an excellent Opportunity to measure the globular
cluster luminosity function fainter than has been done before for any other galaxy. In
order to take advantage of this depth, however, we need to first prepare the images
to ensure that objects are detected with the best possible efficiency. This preparation
mainly involves the removal of the galaxy light. If the galaxy were not removed from
the image, then the photometry of the clusters would be biased by the addition of
extra Signal from the galaxy.

Subtracting the galaxy from the image is essential in another way as well. The
detection threshold is defined in terms of the image noise. This constant threshold
works fine for images with uniform noise, such as a sparse field of stars, but for
these images with their strongly varying noise distribution, the threshold needs to
be defined better. This variable noise arises from the fact that the brightness of any
given pixel is dependent on the number of stars that fall within that pixel. This is a
Poissonian distributed quantity, so the noise due to the galaxy scales as the square
root of the galaxy brightness. The center of the galaxy is thus a much noisier region
of the image than the edges, and so using a fixed detection threshold will miss many
real objects at the edges, and count too many noise spikes in the core. With an
accurate model of the galaxy flux, the detection threshold can be weighted to ensure
equal detection efliciency across the image.

The model of the galaxy light must not be biased by any objects on the image, and
must handle the steep changes near the core of the galaxy. There are two common
methods used for removing galaxies, unsharp masking and isophote fitting. Unfortu-
nately, both of these techniques have failings, and so a hybrid method was developed

to ensure the best quality subtraction possible.

45

4.3.1 UNSHARP MASKING

The unsharp masking method estimates the galaxy light by smoothing the image so
as to remove the contributions from small scale objects. A simple Gaussian filter will
not work to accomplish this, as any final image will still contain the flux from all of
the small objects, just smeared out onto a larger scale. Instead, a median filter is
used, as this can ignore the bulk of the light from these objects. The quality can be
increased by clipping the highest and lowest pixels from the box before the median is
taken. As the main objects found on the image are mainly globular clusters, we can
choose a box size larger than the expected sizes of these objects. This consideration
leads to the choice of a box 100 x 100 pixels in size.

Two main issues arise from this fitting method. First as each pixel in the final im-
age is calculated a different filter box, adjacent pixels may not have smoothly changing
values. Although the majority of pixels in the filter box will remain the same, there
can be sufficient changes to significantly change the output values. This effect will
increase as the image noise increases, to the point where the galaxy subtraction may
actually increase the final image noise. Secondly, bright objects can skew the median,
even if the brightest pixels are excluded. This effect will cause the unsharp masked
image to oversubtract around any bright objects on the image. It also influences how
the core of the galaxy is subtracted, where the steep galaxy profile skews the median
to lower values. This prevents the core from being correctly modeled, and it will show
up as undersubtracted. Figure 4.5 shows the results of unsharp masking on the final
combined F814W image. These defects in the method can be clearly seen, with the
most obvious oversubtraction occuring around the spiral galaxy at the bottom of the

frame.

46

4.3.2 ISOPHOTE FITTING

The ELLIPSE routine in IRAF is commonly used for modeling the light distribution
of galaxies. It fits the galaxy as a series of concentric elliptical sophotes, allowing the
ellipticity, orientation, and center to vary with radius. This method has the benefit
then that all adjacent pixels represent similar points on a smoothly changing model,
which should reduce the noise in the fit. As it fits the values based on only pixels
which match an isophote, it is able to generate a model that accurately accounts for
the steep gradient in the galaxy core. This method can be biased by objects on the
image, but as an iSOphote is generally fit using many pixels, only objects that take
up a significant fraction of the isophote will have much influence. Such objects can
be masked out, which helps remove this problem.

Unfortunately, the isophote fitting can create odd artifacts at radii where EL-
LIPSE has determined the ellipticity has changed. This leads to small ‘fivaves” where
the galaxy light is only partially subtracted, often with matching areas of oversub-
traction on the opposite side of the center. Finally, the ELLIPSE algorithm requires
the isophotes have full angular sampling. This is not true at the largest radii, as the
square image prevents such sampling. As the algorithm stops when this sampling is
not possible, the corners of the image have no fitting done, and hence no subtraction.
Figure 4.6 shows the results of ELLIPSE on the F814W image, illustrating these

problems.

4.3.3 FINAL METHOD

The solution to the difficulties presented is to use a hybrid of these methods. The
image is first scanned by a large box (100 x 100 pixels) and the pixel statistics calcu-
lated within the box. All pixels 40 higher or lower than the box median are flagged,
and the statistics recalculated without them until no more pixels are being flagged.

The box is shifted by half its width, and the process repeated to create a mask that

47

contains all real objects on the image. Since the data region of the final DRZ image
is smaller than the total image, the empty regions created by distortion correction
are also flagged on the mask file.

With a reliable mask created, ELLIPSE is used to find a model for the galaxy.
The fitting algorithm used by ELLIPSE only allows this to work out to a radius
where a full isophote can be constructed. To fill the regions that are not modeled, we
assume that the ellipticity and position angle are fixed at the values of the last fitted
isophote. Incomplete isophotes are then constructed by taking the median values in
annuli of increasing semi-major axis. This extends the fit to cover the entire image.

Although the galaxy model created by ELLIPSE does an excellent job of removing
the majority of the galaxy signal, especially in the sharply peaked core, it still leaves
the small “waves” in the final image. These waves are removed by generating a new
mask from the ELLIPSE subtracted image, and then sampling the image with a wide
median filter and constructing a bicubic spline model between these sampled points.
This spline model is then evaluated across the image, and that difference removed.
With the waves removed, we are left with an image that is fully cleaned of the galaxy.
The final step of the image preparation is the removal of any remaining constant
background offset. This is done by subtracting the mode of the image intensity
histograms from the image. The final mode subtracted histograms are shown in
figure 4.7. The pixel distribution is asymmetric in the galaxy subtracted images,
which is not entirely surprising, as the real objects on the frame should increase the
number of pixels with positive values. The final galaxy subtracted F814W image is
shown in figure 4.8, which shows the clear improvement over the other methods.

The output of this filtering is saved as the “data” image, which has had the galaxy
light subtracted off, and is used for all subsequent analysis. The “background” image
is the final galaxy model, and the “noise” image is the galaxy model plus the fixed

read noise for our data. This image is used as the weight for the detection routine.

48

 

 

Figure 4.1 Original F606W raw FLT frame before reduction. Cosmic rays can be
seen, such as the feature between the bright GC and the companion galaxy on the
right edge.

49

 

 

Figure 4.2 Distortion corrected F606W frame. The rhombus shape of the ACS detec-
tor footprint can be clearly seen as a result of distortion correction.

50

 

Figure 4.3 Final F606W image after combining with Multidrizzle. The cosmic rays
have been eliminated from this image, and the interchip gap has been filled with data
from other frames.

51

 

Figure 4.4 Final F814W image after combining with Multidrizzle.

52

 

Figure 4.5 Subtraction of galaxy light using unsharp masking, with a 100 x 100
pixel box. Although it creates a generally smooth image, the fringing effects can be
seen at the eges of the data, as well as in the oversubtraction that surrounds bright
objects. The sharp rise toward the galaxy core skews the median filter, which prevents
pixels within a box that contains the very center of the galaxy from being correctly
subtracted

53

 

 

 

 

Figure 4.6 Subtraction of the galaxy light using the IRAF ELLIPSE routine. Although
the core of the galaxy shows much improvement (illustrated by the fact that the
central dust lanes can now be seen), the two main failures are also visible. The
corners of the image show where ELLIPSE has stopped fitting, and the unevenness
around the core shows where the isophotes have failed to correctly match the galaxy’s
true ellipticity.

54

 

 

5 I 5 | I l l I I l |
4— - 4
1
“53— —%3
H H
x _x
1 i
z?— ‘22
1— — 1

 

 

 

 

 

 

91000 -500 0 500 1000 1500 0 —3000 -1000 0 1000 3000
counts counts

Figure 4.7 Image histograms for the F606W and F814W images. The centering about
zero indicates the background subtraction has correctly removed the galaxy light. The

asymmetry points out that these images contain real objects that are not subtracted
by the galaxy removal.

55

 

Figure 4.8 Final galaxy cleaned image. Although some artifacts of the galaxy sub-
traction process can be found (such as around the companion galaxies), a clear im-
provement is evident over the simpler methods. This image also shows that the noise
remaining after galaxy subtraction increases towards the center of the galaxy.

56

 

CHAPTER 5:
GLOBULAR CLUSTER LUMINOSITY
FUNCTION

 

5.1 OBJECT DETECTION

Source Extractor (Bertin 8; Arnouts, 1996) was used to generate a database of objects.
The galaxy subtracted “data” images were used as the search images, and the galaxy
model was used to weight the detection process. By default, Source Extractor looks
for objects a specified number of standard deviations above the background, taking
the statistics across the entire image. As the noise in this data is highly dependent on
the galaxy flux, this method does not work uniformly across the image. By supplying
a weight image, however, Source Extractor scales the detection threshold across the
image to reflect how the image noise truly behaves.

For this project, a detection threshold of 30 was used, with a minimum area of
2 pixels. This area criterion requires that an object must satisfy the threshold on
at least two adjacent pixels, and helps to keep noise spikes from being detected as
real objects. We also set the requirement that an object must be detected at this
threshold on both the F606W and F814W frames. By requiring all objects be found
in both filters, we reject any unusual features that appear in only one. Table 5.1 lists

the quantities directly measured by Source Extractor for each detected object.

57

Table 5.1: Quantities measured for each object in both filters by Source Extractor.

 

 

 

Description Units
(X,Y) position on data image pixels
((1,6) position in the sky (2000 epoch) degrees
m,- apparent magnitude within a fixed radius mag
Am,- error in the apparent magnitudes mag
threshold local detection threshold counts
background local background level counts
Fmax maximum flux value counts
Fiso total flux above detection threshold counts
Aiso area above detection threshold pixels
#threshold surface brightness at the threshold detection level mag / arcsec2
,umax surface brightness at object peak mag / arcsec2
A,- area at isophotal levels: I,- = threshold-(51%;“é‘1—dyy8 pixels
flags any internal flags about measurement problems
FWHM full width at half maximum pixels
stellarity classification ranging from one (star) to zero (galaxy)
r1/2 radius containing half of Pisa pixels
A E major axis length pixels
BE minor axis length pixels
63 position angle of major axis degrees
elongation fig
ellipticity 1 — §§
dmerge distance between V and I image centroids arcsec
Rgal projected distance to galaxy center arcsec

 

58

5.1 . 1 INSTRUMENTAL MAGNITUDES

Source Extractor automatically converts the raw fluxes to instrumental magnitudes,
using the photometric zeropoints for the individual filters supplied. This zeropoint is
also modified by the exposure time of the image, so the instrumental magnitude is
defined as

m = —2.510g10 (F) + 2.510g10 (temp) + Zeropoint (5.1)

The zeropoint is listed in the calibration manuals for the individual detectors, and is
reasonably well calibrated.

Unfortunately, this instrumental magnitude only measures the light within the
apertures defined. A single aperture clearly does not fully account for all the cluster

light, requiring a more complete study of the photometry.

5.2 COMPLETENESS CORRECTION

Even though we have set the detection threshold fairly liberally, and have extraordi-
narily deep data, we still expect that we are not likely to be 100% efficient at detecting
objects. The best way to quantify this detection efficiency is to add simulated clusters
to the images, and then check to see how many of these are detected by searching
with Source Extractor and the same detection limits.

Since the detectability of any given object is related to the surface brightness,
and not just the total object flux, we must also incorporate the sizes of the objects.
Simulated globular clusters with a fixed central potential of W0 = 5 (King 0 = 1.03)
were generated for a grid of apparent instrumental magnitudes and tidal radii. At each
grid point, 200 simulated clusters were randomly added to the background subtracted
images. The detection and measurement was repeated as was done for the real data
with each simulated cluster stored along with the values of the parameters calculated

by Source Extractor.

59

Once this is finished, all objects at a given grid point are analyzed, and the ratio
of the number found to the number input is calculated. Due to some significant
objects that remain on the image after galaxy subtraction (most notably the central
jet and the companion galaxy), some regions of the image are manually excluded,
which changes the input number between grid points.

Figure 5.1 shows the completeness surfaces for the two filters. Over the range
of tidal radii that we expect to find globular clusters, the size dependence of the
completeness is fairly weak, so for all further analysis, we define the completeness
solely as a function of apparent instrumental magnitude. Due to the radial dependence
of the galaxy noise, the detection efficiency changes with distance from the core. To
minimize the scatter in the completeness, we divide the data into two radial bins,
breaking at the median cluster distance 68’.’ 95. The final completeness for each filter
was then calculated independently for each bin. This allows the completeness in the

outer bin to extend to slightly fainter levels.

5.3 PHOTOMETRIC CALIBRATION

5.3.1 APERTURE CORRECTION

The completeness values measured from the simulated clusters require the instru-
mental magnitudes of an observation to be calculated. The magnitudes measured by
Source Extractor are taken at fixed radii. As the sizes of real clusters are not all
the same, taking a single aperture magnitude for all objects will not equally measure
the total light. This fact suggests that we need a way to correct the fixed aperture
photometry to account for the variable cluster size. This is done with the aperture

correction. This correction is parameterized by an estimate of the size of the object,

60

 

 

I I

ll. / o [0 9/1 9 I‘ p o ’9 o o 1
o o I o ‘9 9‘0 0 o lo 9] o o o o
l
\o o o ‘6 o ‘0 o ‘0 o“ .l o o o .—

\ 9 k 0 \o o \o o \o ‘9 0‘ 9 o 0
cl 0 \ 9 3 o X o b r o ‘0 o ‘o

7 i L I. '. '
0‘ o 1 o I o ‘ o ‘ o [o OI 9'9
0‘ o \ o W o x o 1 o ‘9 o ‘0 o .. %
o\ o 9‘ o o\ o o\ o 0‘ o ‘9 o ‘0 o

oo‘ov‘oooooc-

I I
o¢¢¢o
‘\

ol

o oi. 9'0 0|. 0

o 0 0‘9 9‘4- 0‘9 9‘9 9‘9 9‘. oi-n
o

 

\ \ ‘
o \b o b o >
0“. \ \o o\¢ \ \
o o\o o\ o\ o X o o\o o\o 0M»

 

 

 

 

 

 

020422 l 24 ' 26 ' 28 ' 30 018 I 20 ' 22 A 24T26 . 28
mF606W mF8l4W

Figure 5.1 Completeness levels as a function of instrumental magnitude and input
cluster size. The lines denote constant completeness, from 100% to 10%, with the
completeness of the line equaling the data completeness. Note that the radius depen-
dence is fairly weak, suggesting that the instrumental magnitude is the main factor
in the completeness. The crosses show the locations of the simulated clusters used to
estimate the completeness.

based on the measured magnitudes with different aperture radii:

R = m4 pxl — m2 pxl (5'2)

The logic behind this parameter is that a small object (with a radius smaller than
2 pixels) will have approximately the same magnitude in both apertures, yielding a
value of ’R. ~ 0. As an object increases in size, more light is measured in the large
aperture compared to the smaller, which pushes ”R. away from zero (and to more
negative values due to the definition of magnitudes).

Given this size parameter, and one of the aperture magnitudes (taken as m4px1

61

for convenience), it should be possible to construct an aperture correction, such that

minstrumental : m4 pxl + “MRI m4 pxl) (5-3)

Conveniently, the data used for the completeness correction samples this function
A(’R., m4 pxl): as the simulated clusters were created with a known instrumental mag-
nitude, and ”R, and m4px1 are measured by Source Extractor. At each simulated
cluster grid point, all detected clusters are measured, and the median values of ”R

and m4px1 are stored at that point, with the value of the aperture correction:
A (median(R), median(m4 pxl» = minput instrumental - median(m4 pxl) (5:4)

This grid of points has the upper end fixed, such that A(0,m4pxl) = 0 to anchor
the small object end of the aperture correction (which is not well sampled by the
completeness data). This grid of values (72, m4”), A) is irregularly sampled as the
values of ’R do not linearly match the input sizes. For each detected object, the
aperture correction is calculated using thin plate splines, which can interpolate such
data. Figure 5.2 shows the surfaces of the aperture correction for the two filters.

In addition to this aperture correction, another correction of 0.1 magnitudes is
applied to the final instrumental magnitude. This extra aperture correction compen—
sates for light that is scattered to large angles by the optics of HST. This is a standard
calibration step, and must be applied even though an aperture correction has already
been applied. The simulated clusters used do not incorporate this scatter, as the PSF

used ignores very wide angle scattering.

5.3.2 COLOR CORRECTION AND EXTINCTION

Most previously published results on globular clusters present the magnitudes using

the standard Johnson Cousins BVRI system. Converting the instrumental magni-

62

 

 

 

 

 

 

 

 

 

K
I I I I I I | 'O-U I I I I I I I I
h .
1". o o o _,-¢‘ {-“nngz ; ...-f4". ..... 2 ...._°'__ F ....... o .3 o 3 3 : 3 : .x . .’ e "9. 3., ,
' * o “om-9"": 0 3 ..3 :‘-._’_..~.o*. 9, o 3.. o. o- o l) v ’3 3 '1' .1 2 . ” v e . "4
.»o' ,1 .9. o. f I D o I ‘ t OI 8", O- 3:“ .. -
’~ ‘ ‘1 .‘. I ': , . I . *: L ». 2 ~:."".. .
" “C. '2»: _ ‘ ' ' . 7' .‘ I- "rm-'A‘i; .,-' '.-' '- ‘ .'
J (I;.>:."{Z::"“If.-"ZI _ ~91 ,‘ '4 0 O o o o I, o 1*. 3 _‘;,"I,:_'~'-‘:f" k._ ‘0‘... ‘0’: j'..‘_ .0 . .
' -~ .I .. . ' e "..I ' " ~ ‘. ’ A n . -' '1'1-— '« V L. I ‘
”.... 1' __ '.:‘ 1L _. . . . g , V... =.' J N...” ,r_,. v.
-. «,1» .. «... '1' ~' _. . , - _ 5nd, gf'»...__ ‘ , ,
g1!) £";“:' If. - '__‘ @100” 3 ' ° ‘ 9| 8 3» 3 2 6,}? .’ 1‘”: 1‘ F" .’ 0‘ o 0"
.x- * ~ . ..... . z z :m . '3’ ~. .
¢ ..... * . z. :. r : 3 .
3' :.-9 :0. o o o O .9 Q 3' gr}; ° ’0 o ‘9 o o o
.. ‘ ..
ire»: ' 9* °’ < __°__°_ :::: o o t. 3.4—3-: * o ... . . 4
:g o 1o ., o 3"— > o .-.__.__.-z z z , v 9 ,'o_ 9 3"3‘
3 ’01 o o '5 1 . o o o o o 9 ‘0' 3 z o 6 o o o o o
.9 n 0.; o . . x o o o o . . g : ‘ ‘ o o o .0 . .
t h. 0 o ’9 9 q 0 O o . . ‘ 3 g. 0 O 6 .29 ,
0 o .......... o o o o 4'0 o
1"."2"--‘-:—: L—————.—..__——"b—"-:;-‘ ~'—v‘fi ——————— .
. .
.0 0 o o ‘ 9
o
5-....9r-‘3 ........ _+ .— ....... : ’. o _....’—'.‘ ...... _‘
‘ _________ o:
g 0. 3 03. o o. o + o o o o;-rr ._‘£,.: 03 0 o o a. o o -I
or ------- fi. uuuuu d P --------- “z .......... . ...............
t 4 _.t-“'*--...*.....9.,..t....t. ~° * ‘ r..:....?.,.t. 4...: ° ‘. ......... f.....‘.....‘:e, 3..
01”.. O o o *3 o "“0 o o o o ’0 in... ...... ‘ 0 9 o ‘9 h"
K J l + t: l l 4 l I l l l
- .1 - .1

Figure 5.2 Aperture correction levels as a function of instrumental magnitude and
input cluster size. The weak dependence on the measured magnitude indicates that
the fraction of light lost is generally a function only of the object size. The crosses

Show the where the simulated clusters used for evaluation were placed, and the small
dots show where the real globular cluster data falls.

tudes mF606W and mF814W to V and I requires a color correction. The correction

parameters are taken from Sirianni et a1. (2005) and reproduced in table 5.2

Table 5.2: Color correction parameters

 

 

al

a2

 

Filter (V — I )break a0
F606W < 0.4 26.394 :1: 0.005
> 0.4 26.331 :1: 0.008
F814W <0.1 25.489 :t 0.013
> 0.1 25.496 i 0.010
F775W < 1.2 25.241 :l: 0.005
> 1.2

25.292 i 0.033

0.153 d: 0.018
0.340 :I: 0.008
0.041 d: 0.211
-0.014 :I: 0.013
-0.061 :1: 0.021
-0.105 :I: 0.026

0.096 :I: 0.085
-0.038 :1: 0.002
-0.093 :1: 0.803
0.015 :I: 0.003
0.002 :I: 0.021
0.007 :t 0.004

 

63

The corrections are defined as:

V = mF606w+00+al(V—I)+02(V-1)2

I = mF814W + b0 +bl(V—1)+b2(V — I)2

and are defined piecewise in color on two intervals separated at (V — I )break- The
standard application of this color correction uses an iterative process, where the color
V — I is updated at each step, and V and I are re—evaluated with that new color.
Unfortunately, the piecewise definition of the color correction creates problems in the

evaluation around the breaks. Instead, an algebraic solution was used, noting:

V— I = (mF606W — mF814W) + (00 — b0) + (GI - bl)(V— 1) + (a2 -— b2)(V — [)2 (5-5)

which can be solved for a corrected color in terms of the instrumental color. Formally,
this yields two solutions, but it can be seen by checking a sample of colors that only
the negative solution gives realistic magnitudes. Although this method seems to have
the same problem at the color correction breaks, we can check that the calculated
color falls within the necessary range for the coefficients used. Figure 5.3 shows a plot
of the calculated V — I color as a function of the mF606W — mF814W instrumental
color. As this function is continuous even around the color breaks, we can be sure
that the color correction has been applied correctly.

The final photometric step is the application of an extinction term. This term
accounts for the absorption of light as it passes through dust and gas in space. Using
the Schlegel et a1. (1998) value for M87’s position on the sky, we can estimate this

extinction as:

AV = —0.074 (5.6)

A; = —0.043 (5.7)

64

 

 

 

 

 

0 I 1 0.5 1.0 I 1.5 ‘ 2.0
mF606W — mF814W

Figure 5.3 Plot of calculated V — I color as a function of instrumental mF606W —
mF814W color. The continuity ensures that the correction has been applied correctly,
even at the color breaks, which are marked with lines.

5.4 CONTAMINATION

Even with the various cuts placed on the data, it is still possible that some non-
cluster objects may remain in the sample. Although these objects masquerade as real
clusters, we can statistically remove them by estimating how many of these objects
should appear on our image. To do this, we need to find images that only contain

the contaminating objects.

5.4.1 NOISE OBJECTS

One method to look for noise spikes on the image is to run the detection code on an
inverse image, generated by multiplying the data image by —1. The real objects on
this inverted image have peaks that are then negative, and so fall below the detection

threshold. Instead, only the dark areas of the data image will have peaks that allow

65

them to be detected. If the image histogram was perfectly symmetric, than the
number of objects detected on this inverse image would directly give the number of
noise spikes in the data image. The fact that the histograms Shown in Figure 4.7
are not symmetric suggests that our data is deep enough to measure a bias in how
the stars in M87 are distributed. Despite the fact that we cannot use these objects
directly to correct for the contamination, we can still use this as a lower limit, and

correct for the points that are detected.

5.4.2 BACKGROUND GALAXIES

The contamination by background galaxies is obviously a major contribution to the
number of objects found. A quick look at the galaxy subtracted data image shows a
large number of background spiral galaxies that are easy to identify. The Hubble Ultra
Deep Field (Beckwith et al., 2006) is a very deep data image taken of an “empty” field
with a goal of looking at the distribution of distant galaxies. This gives an obvious
way to estimate the contribution from these galaxies.

Unfortunately , simply searching for UDF galaxies will overestimate this contri-
bution. Instead, we need to ensure that the noise characteristics of the UDF match
that of the data. This first requires a rebinning of the UDF data from a pixel scale of
(If/030 to the pixel scale of the data (0’.’045). Next, the UDF is overlaid on the image
footprint of the data, and trimmed to match the field of view. With the geometry
matched, the noise characteristics can be matched to our data. The UDF is given in
counts per second, so must first be scaled to match the exposure time of our data.
Finally a noise image is created from the data galaxy model. This noise model as-
sumes each pixel is drawn from a Gaussian distribution with mean zero and standard
deviation equal to the square root of the galaxy level. This noise model is then scaled
such that when added to the UDF, the final image variance equals that of the data

image. This method makes it certain that objects in the UDF frame that fall “near

66

the center of the galaxy” have more noise than objects that fall at the edge. Figures
5.4 and 5.5 show the original UDF image, as well as the final cropped and noise added
final version. The apparent loss of objects from the center of the image is a result of
the increase in the image noise in that region.

The reduction is identical for objects detected in the UDF as for the data objects,
with the only exception being that the UDF filters are F606W and F77 5W. However,
we can convert these filters to the same standard V and I system, minimizing the

changes.

5.5 CLUSTER CANDIDATES

The cluster detection code used so far has been designed to detect as many objects
as possible. Unfortunately, this also means that some fraction of the objects detected
are not true globular clusters. The two main sources of contaminating objects are
background galaxies and noise fluctuations from the galaxy light (including random
local overdensities in the stars that make up the galaxy). To eliminate as many
of these contaminating objects as possible, we make a series of cuts in the many
dimensional space of parameters.

The first cut is to remove objects that have dmerge > 0’.’ 015. This removes things
that are likely to be random matches of noise between the two images, that do not
represent real objects. Next, we can apply a color cut where only objects with
0.5 < V — I < 1.7 are considered to be globular clusters. These limits are based
on observations of other globular cluster systems and cover the range of colors where
clusters are known to exist. Finally, we only want to consider objects that have a
final completeness value greater than 0.5. This limit is chosen because if we are only
detecting half of these objects, then any count of them is influenced more by the
errors in the completeness than in the true number of objects.

These cuts so far are all based on obvious ways to limit the sample. We can

67

also make the assumption that globular clusters are circular, and set limits on their
ellipticity such that g VI < 1.5. We also assume that clusters are peaked, such
that the central pixel is bi‘ighter than the rest of the cluster. By comparing the peak
surface bright to a calculated average surface brightness, we can exclude objects where

the average surface brightness is brighter than the peak surface brightness:

000455 2
Hpeak < V + 2-510g10 (AISO ( 1px] ) ) . (5'8)

 

Such an object is likely a very extended background galaxy, or a very blurry noise
object.

The numbers of clusters detected with each of the various cuts are given in table
5.3. This shows that the various restrictions on the data remove the majority of the
contaminating objects, and yields a final number of clusters consistent with what is
expected based on earlier surveys. The one object that is manually excluded from
the sample is the HST-1 knot in the jet, which successfully eludes these cuts. This
is likely due to the fact that it is dominated by a small flaring emission source, that
makes it appear reasonably point like (Perlman et al., 2003; Waters & Zepf, 2005).
Figure 5.6 shows the color magnitude diagrams of both the full sample of detected
objects, as well as the full sample of objects on the contaminating frames. It is clear

that at the brightest magnitudes, the globular clusters are easy to identify.

68

 

Figure 5.4 F606W UDF image rescaled to match the pixel scale of the M87 data.
Nearly all objects on this image are galaxies, making it an excellent field to investigate
the expected contamination from background galaxies.

69

 

Figure 5.5 F606W UDF image cropped to the footprint of the M87 data, and with
added noise to simulate the effects of the galaxy subtraction process. This added
noise can be seen by the apparent preferential loss of objects from the center of the
image, where the noise added has a stronger variance.

70

Table 5.3: Impact of Selection Cuts on Sample Size

 

 

 

 

 

Cut Nclusters Ninverse NUDF
Raw V 9454 721 790
Raw I 12495 249 255
Joint V+I 6624 194 249
dmerge < 01015 5392 194 249
Completeness > 0.5 4045 47 231
V — I < 1.7 3113 26 176
V — I > 0.5 2919 16 114
§|V < 1.5 2434 6 60
§| I < 1.5 2191 4 34
#peak < (p) 2092 2 33
Final Sample 2091 2 33

 

5.5.1 PHOTOMETRIC ACCURACY

As the ACS pointing contains the majority of the clusters imaged by Kundu et a1.
(1999) and Waters et a1. (2006), we can test the photometric accuracy of the survey
by comparing the measured magnitudes to those previously done. For the Kundu
et al. (1999) sample, there are 989 clusters matched between the two samples, with

magnitude differences

V — VKundu = —0.028 :t 0.100
I — [Kundu = 0.033 a: 0.107

(V — I) — (V — I)Kundu = —0.06 :1: 0.110

71

20 . ..

 

 

 

 

 

 

 

 

 

)- _
el 1
24- .
>_

26- -

28- -

30-4‘-‘2'é‘é‘430:4‘J2‘6‘é‘2
V—I V—I

Figure 5.6 Color magnitude diagrams of detected objects. The left panel shows all
objects detected on the data image. The globular clusters are clearly visible as the
vertical stripe around V — I ~ 1.0. The right panel Shows the objects detected in the
inverse (shown as asterisks) and UDF (shown as boxes) images. The objects detected
in the data image, but excluded from the sample due to their properties are also
shown as small crosses.

For the Waters et a1. (2006) sample, there are 814 clusters in common, and the

differences are

V — VW,tars = 0.035 :t 0.102
I — 1W,ters = 0.073 :t 0.086

(V — I) - (V — I)Waters = ——0.035 :t 0.127

These deviations are consistent with the expected photometric errors.

72

5.6 LUMINOSITY FUNCTION

The luminosity function for the globular clusters shows how many are detected at a

given magnitude. Each detection is weighted by the a factor

__ i1
_ Completeness,

 

z' (5.9)
where the positive case is for objects from the “cluster” sample, and the negative case
for the “false” and “UDF” samples. These samples are included to correct for the
expected contaminating objects that still exist in the sample. The error in each bin

is then

Ebin = (II/VII) V Nbin (5-10)

Figure 5.8 shows the final GCLF, including both the completeness weighting and
the contamination rejection. Overplotted with the GCLF is the smoothed complete—
ness as a function of magnitude. This shows that our sample is 50% complete to
nearly V ~ 27. As this is the product of the completenesses in both the F606W and
F814W filters, the final completeness is color dependent. For the range of colors we
consider for our globular cluster sample, the 50% completeness limit generally follows
the line V50% ~ 25.9 + 0.8 (V — I).

GCLFS have in general been modeled by Gaussians, due to their peaked natures.
Recently, more physically motivated models have been used, based on the expected
mass lost from the clusters. Such models are discussed further in Chapter 6. The
best fitting Gaussian to the binned GCLF gives a turnover of (V) = 23.62 with a
Width of 0V = 1.40. Repeating this for the I data yields (I) = 22.58 with a] = 1.35.

5.6.1 COLOR DEPENDENCE

The M87 globular clusters are bimodal in color, as are the GCs in most galaxies. This

bimodality is obvious from the color magnitude diagram, in figure 5.7, with a blue

73

peak around V - I ~ 0.9 and a red peak at V — I ~ 1.1. If we divide the clusters
between two samples in color separated at V — I ~ 1.03, we can create luminosity
functions for each. This divisions places 959 clusters in the blue sample, and 1133 in
the red sample. Figure 5.9 shows the GCLFS for these two samples, along with their
associated completeness functions.

It is clear from this figure that the faint end of the blue sample falls off faster
than the red sample. This drop is somewhat expected based on the color dependent
completeness function, but as shown, the fallofl' occurs at a level where the sample
should still be roughly 90% complete. Therefore, this drop must be a real effect, such
that the faintest clusters tend to be from the red sample. On the bright end of the
GCLF, the blue sample also drops off at brighter magnitudes. Fitting Gaussians to
these two samples confirms these eflects, with the blue sample having a higher mean
than the red sample, with <V>blue = 23.31 and <V>red = 23.89. In addition to the
differences in the means, the widths of these Gaussian fits are also different, with the
blue sample significantly narrower than the red sample: Ublue = 1.21 and 0red = 1.40.

We can check that these fits are not just an effect of the sampling by comparing
them with the t— and F— tests. These tests indicate that the means of the two
samples are significantly different (t = 10.042, p-value = 1.6 x 10’”), as are the

variances (F = 1.34, p-value = 2.5 x 10—6).

5.6.2 RADIAL DEPENDENCE

Although we can only measure the projected distance from the clusters to the center
of M87, any radial dependence of the cluster luminosities should still be evident.
To investigate this, the clusters were divided into five bins in projected radius with
widths of 2 kpc. Table 5.4 shows the Gaussian fits to the GCLFS for these bins, along
with the number of clusters in each bin. The fall in the number of clusters beyond

RGC ~ 6 kpc is due to the shape of our field, and not from a sudden drop in the

74

cluster density. Figure 5.10 plots the binned GCLFS as simple lines, along with the
completeness functions for each bin. We can see that the peak of the GCLF moves
fainter with increasing distance from the center of the galaxy, along with an increase

in the GCLF dispersions.

Table 5.4: Radial Dependence of the GCLF

 

 

Bin NGC (V) 0’

 

0 < R < 2 186 23.30 0.985
2 < R < 4 481 23.46 1.246
4 < R < 6 569 23.56 1.280
6 < R < 8 568 23.73 1.437
8 < R < 10 244 23.75 1.488

 

We can again use the t— and F— tests to confirm that this radial gradient is a
real property of the GCLF. Comparing each bin against the next larger indicates that
only the outer two bins have consistent means (p ~ 0.43). This is likely due to the
fact that the outermost bin preferentially samples the smaller radii due to the square
shaped frame, so it is sampling clusters most like its neighboring bin. The variances
also match for the two outer bins (p ~ 0.53) and also for the R2_,4 and R4_,6 bins
(12 ~ 0.53). Again, the differences in completeness make an attractive explanation for
this, but as before, the data is complete beyond the point where the GCLF begins to
fall.

5.6.3 COMPARISON WITH PUBLISHED RESULTS

A selection of published Gaussian GCLF fits are presented in Table 5.5. Figures 5.11
and 5.12 show these fits overplotted withithis new ACS GCLF data for both filters.

75

The first GCLF to compare our data to is based on the Harris (1996) catalog of
Milky Way globular clusters, using the most recent (2003) revision. The histogram
of this data is shown in figure 5.11, illustrating the statistical advantage found in our
sample. The peak for the Milky Way is only marginally fainter than we observe in
M87, with a t—test p value that suggests they are not significantly different. The vari-
ances, however, are very much different, with the Milky Way GCLF much narrower.
This variance is in fact consistent with the M87 blue sample, which is reasonable as
the Milky Way GC system is composed predominantly of blue clusters.

For our V GCLF, we compare against the two WFPC/ 2 surveys (Kundu et al.,
1999; Waters et al., 2006) as well as one ground based survey (McLaughlin et al.,
1994). The WFPC / 2 samples match well, with reasonable t— and F— test p values.
Incorporating the photometric offsets calculated in Section 5.5.1 shifts the turnover
values even further into agreement. The McLaughlin et a1. (1994) sample has much
brighter completeness limits, such that they are only sensitive to the bright half of
the GCLF. This leads to a significant problem in fitting a Gaussian, as is seen by the
fainter peak and wider dispersion.

The agreement with the I GCLF is not quite as good, likely due in part to the
worse completeness effects for the I data. The Kundu et a1. (1999) and Waters et
a1. (2006) samples are again reasonably consistent, with the p values improving again
when the offsets between the samples are applied. The comparison with the ACSVCS
sample (Jordan et al., 2007) indicates Significant disagreement. This is likely biased
due to the fact that the ACSVCS used the F850LP (z) filter, which extends further
into the infrared that our F814W. To compare the samples, we’ve used a simple
conversion of I — z = 0.3. This ignores any more complicated color dependence of
the correction, which may account for the narrower dispersion in the ACSVCS fits.
However, it seems far more likely that the ACSVCS dispersion is an underestimate,

based on the lack of the very faintest clusters in their sample.

76

 

 

 

 

19

 

V'—I

Figure 5.7 Color magnitude diagram of the final sample of clusters with color his-

togram. The blue and red clusters can be seen as grouping around V —— I ~ 0.9 and

V — I ~ 1.1 respectively.

77

100

 

 

75

01
O
ssauaqudquQ

25

 

 

27

 

Figure 5.8 Final globular cluster luminosity function along with the completeness
function. This shows that the sample is at least 50% complete down to 27th mag-
nitude, Where the cluster distribution is seen to drop to roughly zero. The dashed
curve represents the best fitting Gaussian, with u = 23.623 and a = 1.40.

78

 

 

100.....fifi,.,

 

 

 

 

Figure 5.9 Comparison of the red and blue cluster GCLFS. The Gaussian fits are
plotted, along with the smoothed completeness curves for both samples.

79

 

 

nn
UV

250-

25-

 

 

 

 

Figure 5.10 Comparison of the GCLF between five bins in projected distance from
the center of M87 along with their completeness curves. The lines are for R042 (red),
R2_.4 (blue), R445 (green), R643 (magenta), and R8_,10 (cyan).

80

 

 

 

 

$2 BER em 2 38 H Bass 25. $2

83 Ream a 2 38 > ease... was. 5:

TS x 3 $2 8-2 x was $3 $3 23m 8 2 BE 8 A38 .5 B 52% $2
5.3 was So $3 92 RR R 2 22 H 883 .B s 8883 52
25.: £3 :3 83. was Ram R 2 was 2 88¢ a 8 35m 52
ex; 83 S3 was a: 8% ea 2 22 > 883 a 5 5.583 $2
2.3 was $3 $3 a: $2 em 2 as > 38: a 5 Beam $2
88 m5 2:. 82: m: 8.8 R 2 $3 > :83 a s £55352 5:
TS x we 23 82C 33 w: E? R 2 ...2 > 8%: 25m .55 .552
a SEQ a ease e a sea 82 48E 8888a 5380

 

 

BOO 2: 3 BE 5330 as 2.3.

81

5.7 BLUE TILT

Recent deep observations of extragalactic globular cluster systems have suggested the
existence of a “blue tilt,” where fainter clusters from the blue metal poor sample are
bluer (and hence, more metal poor) than the brighter blue clusters (Strader et al.,
2006; Harris et al., 2006). This tilt is explained as a result of a metallicity-luminosity
relation of the form Z 0: L055.

We can look for such a trend in our data by dividing the sample into bins of width
one magnitude, and running a KMM (Ashman et al., 1994) test on these samples.
This test looks for bimodality in the sample, and finds the best fitting peaks and

dispersions for the two groups. Table 5.6 shows the results of this test, and the

results are plotted over the globular cluster color-magnitude diagram in figure 5.13.

Table 5.6: Bimodality with Cluster Luminosity

 

 

Bin NGC %blue I-‘blue ablue %red ”red cred

 

V < 20 9 0.22 0.59 0.11 0.78 1.12 0.19
20 < V < 21 63 0.53 0.85 0.09 0.47 1.10 0.08
21 < V < 22 229 0.47 0.84 0.08 0.53 1.11 0.11
22 < V < 23 421 0.54 0.87 0.10 0.46 1.14 0.06
23 < V < 24 613 0.45 0.85 0.09 0.55 1.15 0.08
24 < V < 25 459 0.35 0.84 0.10 0.67 1.19 0.13
25 < V < 26 224 0.49 0.95 0.15 0.51 1.33 0.13
26 < V < 27 74 0.48 0.98 0.16 0.52 1.42 0.14

 

The brightest bin contains too few objects to give reliable results, and the color
dependent completeness can be seen to deplete the edge of the CMD below V ~ 25,

making the faintest two bins also unreliable for investigating the color bimodality.

82

The remaining bins have almost identical values for ”blue: indicating that our sample
shows no evidence for a blue tilt. AS this sample is much deeper (generally by at least
a magnitude) and more complete than those that appear to have a blue tilt, it seems

likely that the tilt is not a real effect, but rather some systematic photometric error.

83

nnn
QUU ' I ' l ' I ' I ' I ' I ' I

 

 

 

 

Figure 5.11 Comparison of the final GCLF to the M87 GCLF of Waters et al. (2006)
(red histogram) and the Milky Way GCLF based on the data of Harris (1996) (black
histogram). This comparison shows the number advantage this new sample presents.
In addition, the truncation of the Waters et al. (2006) sample at V ~ 26 can be seen
to limit the faint end of the GCLF. Overplotted are the best fitting ACS Gaussian
fit (dashed line), the Milky Way fit (dotted line), and the ground based McLaughlin
et a1. (1994) fit (dot dashed), rescaled to match the number of ACS clusters.

84

 

200 ' I ' I ' I ' I ' I ' I I I

150

2100

50

 

 

 

018

Figure 5.12 Comparison of the final I band GCLF to the M87 GCLF of Waters et
a1. (2006). The Gaussian fits are for the ACS (solid), WFPC/ 2 (dashed), and the
ACSVCS z filter (dot dashed).

85

19

 

_I_.

21

>23

I

25

 

 

Figure 5.13 Bimodality of the cluster sample as a function of

which we expect to see completeness issues influence the fitting.

86

magnitude. As the
peaks of the blue cluster distribution stay constant from V ~ 21 to V ~ 25, we see
no evidence for a blue tilt. The black line illustrates the 90% completion limit, below

 

 

CHAPTER 6:
GLOBULAR CLUSTER MASS FUNCTION

 

The mass spectrum of globular clusters is defined as 2/2(M)dM, which gives the number

of clusters with mass between M and M + dM. We then define the mass function as

= M WM )
1081009)

 

‘1’ (10810(M)) (6.1)

to directly relate to observational counts of clusters. In the literature, both 29 and \II
are interchangeably called the “mass function.”

For an initial mass function Ibo E «MM, 0) = 212(M), we can see that the evolved
form of this must be 7,!)(M, t) 0: 2120(M0) where M0 is the original mass of the cluster
at t = 0. Considering the mass loss mechanisms of section 2.3, we can see that for

both stellar evolution and gravitational shocks, the mass loss has the form

M (t) _,,t
— = 6.2
M0 6 ( )
whereas two body relaxation yields
M(t) = M0 — at (6.3)

Since the form of both stellar evolution and gravitational shocks simply scales the
initial mass function, we can neglect the effects of these mechanisms, and instead,

simply allow the normalization to change over time. This defines the evolution with

87

respect only to the number of clusters that survive. Evaporation does not follow this

form, and so must be treated by actually evolving the cluster mass spectrum.

6.1 INITIAL MASS FUNCTION OF GLOBULAR CLUSTERS

For this study of the M87 globular cluster system, we want to consider two forms for

the initial mass function. The first IMF we consider is a power law, where
(MM) oc M—a (6.4)

The standard choice of a is a = 2.0, as this is what is observed to be the luminosity
function of young star clusters in merging galaxies (Zhang & Fall, 1999). Since these
merging galaxies are thought to be in the process of creating new globular clusters,
this observed function is likely to be a reasonable estimate of the mass function of
all systems of newly formed clusters. This IMF is also reasonable from a theoretical
standpoint, as it assumes that the collapse of the giant molecular clouds that form
GCs is effectively scale free.

Unfortunately, such power law mass functions overpredict the number of high mass
globular clusters in many evolved systems. A solution is to use a Schechter (1976)
like function, in which above some cutoff mass, M0, the number of clusters decays
exponentially:

’I/J(M) oc M’arfM/MC (6.5)

This form retains a power law shape at low masses, but falls off to more closely match
observations at higher mass. Burkert & Smith (2000) present this form of IMF as the
result of the expected mass function for the sizes of giant molecular clouds that form
from the merging of smaller clouds. They give a = 1.5 and also provide a fit to value
for the cutoff in M87, M0 = 5 x IOGMQ. We use this form for the fitting, although

we allow the cutoff value to be fit to match the data.

88

6.2 MASS LOSS RATES

If we define the mass loss rate as

dM
—_ = _ —7
dt kM (6.6)

we would like to find a form for the total mass function as a function of time, given
the initial mass function 1120(M). As the mass loss equation is so simple, we can solve

it by direct integration:

MMM = —kdt
M(t) t
M’YdM = —k/ dt
M0 0
; 1+7_ 1+7 _ _
1+7(M(I) M0 ) _ kt

So, the cluster mass at any time is just
1+
M(t)1+7 = M0 ’7 — kt(1 + 7) (6.7)

Since we observe the evolved mass, we can invert this to find the starting mass for
a given cluster, which gives the relative weights for that cluster from the IMF. This

gives the mass Spectrum for evolved clusters as

wevolved(M) = VIMF (M0) (6.8)

= my((M(t>1+l+kt<1+v))1“””) (6.9)

This evolved mass function depends on the product of the mass loss rate and the
current time of evolution. As we only know the current ages of the clusters indirectly,
it makes sense to consider this as a Single parameter, 7' = k - t. This represents the

total mass lost by the cluster over its entire life. This also provides the peak value

89

of the GCMF, as clusters at the peak have lost roughly half their mass, and so have

current masses equal to 7'.

6.2.1 FALL & ZHANG

Following Fall & Zhang (2001), we can consider the mass loss of a globular cluster
to be due mainly to the effects of two-body relaxation. As stellar evolution happens
quickly, we can assume that it happens immediately after the clusters’ formation, and
as such, can be treated as a simple shift of the IMF. Similarly, the mass loss due to
shocks affects the mass function as a whole, and can also be ignored as a simple Shift.

Using the general form above, we can write the mass loss as

dM

d—t = —kM0'0 (6.10)

where is k is the #60 of evaporation defined in section 2.3.4.

6.2.2 BAUMGARDT 82 MAKINO

We would also like to consider the model presented by Baumgardt & Makino (2003),
based on N-body calculations of multimass globular clusters. These simulations in-
cluded all of the mechanisms presented in section 2.3. Baumgardt 82 Makino (2003)

parameterize the dissolution timescale of clusters in their simulation as

_ :1: 1—36
tdiSS _ ktrhtcrossing (6-11)
1—3:
1 2 3/2 ‘5 3/2
N / rh rh

 

 

Gl/2m1/21nA G1/2N1/2m1/2

k l x52
lnA VG

90

Noting that N = 1%? we can see that this gives
tdiss 0C M33 (6.12)

We can convert this to the general form by noting

0 tdiss
/ —kM7dM = / dt (6.13)
M0 0
1 +1

so, 7 = :1: — 1. The simulations they performed suggested a range of x from 0.75 to

0.85, with 0.75 being their preferred value. This gives a general mass loss:

M
dEI— = —kM“0°25 (6.15)

Some authors have suggested that the correct way to convert the Baumgardt &
Makino (2003) dissolution time to a mass loss is to set the mass loss rate It o< Mg , and
leave that fixed for the entire lifetime of the cluster. For an evolved system like M87,
this will give a result very similar to the constant mass loss rate M = kMO°0. All
massive clusters (M > 7') will have lost such a small fraction of their mass that they
Show no change. The least massive clusters (M << 7') will also have effectively the
same mass loss rate, as they all come from progenitor clusters with nearly identical
masses. Only around the turnover (M ~ 7') will there be much change, where the
progenitor masses are sufficiently different to create changes in the mass loss rate.
This will have the effect of creating a slightly more peaked mass function, although
the scale of this effect is still very small.

A final check that setting the mass loss rate as a function of the current cluster
mass is correct can be seen by overplotting the evolution curve of this model against

the evolution curves plotted by Baumgardt & Makino (2003) of their N-body simu-

91

lations. Figure 6.1 does this, illustrating that both this model and the N-body data
drop faster than linear early in the cluster’s life, and then slow down to a slower than

linear loss in the later stages of evolution.

 

1.0 .'<,'~ f I I l I I I I I I

\ \
\ '4
x .‘
\ \ ..
0 8 _ \ X _.
o \ ‘ .'
\ “ 5‘
\ ‘.
x ‘3

M/Mo

0.4

I
l

I
/

‘. ~, ‘
~, '. \ \
.‘~ ._~ . ~. ‘
7.” ‘-_ ~. \
‘~_ ‘-, x \
‘-. ‘-_ ‘ \ -
‘. “ ' \
. ‘._ ‘._ ‘ \ \
‘. '-‘ ‘ '\
‘. ', ‘ \
-lV ‘._ ' \
'._ -.‘ \

 

 

 

0‘80 ' 0.2 ' 0.4 ' 0.6 l 0.8 ‘ 1.0 ”1.2
t/tdiss

Figure 6.1 Comparison of M or M ‘0'25 mass loss model to various Baumgardt &
Makino (2003) N-body simulations. The deviation around t/tdiss = 1.0 arises from
the differences when the number of stars in the cluster becomes small.

6.2.3 LAMERS ET AL.

The final published model we consider is that given by Lamers et a1. (2005). This
model builds upon the results of Baumgardt & Makino (2003), by noting that as In A
also depends on the number of particles, it must slowly vary over the life of the cluster
as the number of stars drops. They use a polynomial fit to this function, adopting

the form

N 0.75 N 0.75

This polynomial is equivalent around N ~ 3 x 105, which corresponds to a mass

92

Mequiv ~ 1.2585, with a slight overestimate of the dissolution time for lower masses.
Lamers et al. (2005) also notes that for an object moving in a circular orbit,

222 _.
7m = V - F = 4770p (6.17)
r
and uses this to substitute the local galaxy density for the second term in the Baum-

gardt & Makino (2003) solution, giving their final dissolution time as

M 0.62 —1/2

where Cam, is a parameterization of the dissolution time. They quote a value for Cam)
of 810, with a range stretching down to 300.

They compare this relation to the star cluster populations of the Milky Way,
M51, M33, and the Small Magellanic Cloud. Although these samples contain large
numbers of young clusters that are unlikely to be bound globular clusters, they insist
the results are valid for all types of clusters. Using the general form again, this has a
mass loss

E = —]€M_0'38 (6.19)

following the conversion from dissolution time to mass loss used above for Baumgardt

& Makino (2003).

6.3 MASS FUNCTION

The observed clusters are used to make the mass function for M87, by first converting

the measured magnitudes to masses, via the relations:

_L__ = 10<V+5—5loele<d)—Ve>/(—2.5>
LO

M L

_ = T—

MO LC

93

where we take the mass to light ratio to be equal to T = 3.0. A histogram of these
masses is created in boxes of logarithmic M, weighted by their completion factors (in
the same way as for the GCLF), yielding the mass function \I!. This mass function is
plotted in figure 6.2, along with the best fitting mass loss models. This figure shows
that the Burkert 82 Smith (2000) IMF fits the high mass clusters best.

The fits to the mass loss models were created by minimizing the x2 deviations for
each bin. The models were allowed two free parameters: A, an arbitrary normalization
defined to be equal to the value of loglo (\II(6)), and T, the total mass loss experienced
by the clusters. For the Burkert 82 Smith (2000) IMF, the value of M0 was also
allowed to vary to ensure the best possible model. This makes only a small difference
in the value, as MC tends to stay reasonably close to the initial value. The best
fitting values for these fits for each mass loss model are given in table 6.1 along with

the reduced X2 for each fit.

Table 6.1: Best fitting mass loss fits to ACS data.

 

 

 

 

 

 

IMF 'y 7' A x2
Powerlaw a = 1.8 0.00 150980 2.275 2.9933
0.25 9346.2 2.277 3.5235
0.38 2273.8 2.275 3.8841
Powerlaw a = 2.0 0.00 202558 2.281 2.4658
0.25 12175 2.279 3.0689
0.38 2927.4 2.275 3.4850
Burkert & Smith (2000) log10(Mc) = 6.4 0.00 178462 2.375 1.3648
0.25 11950 2.382 1.6782
0.38 3238.6 2.390 1.8368
Burkert 85 Smith (2000) log10(Mc) = 6.376 0.00 186644 2.382 1.3512
10g10(MC) = 6.343 0.25 13481 2.399 1.6053
log10(MC) = 6.340 0.38 3772.2 2.408 1.7451

 

94

 

 

 

10310 (‘1’)
10810 (‘1’)

 

 

 

 

 

 

 

l l l l I 1 l 1 l L
0 4 5 6 7 0 4 5 6 7
10810 (M / MO) 10810 (M/MO)

Figure 6.2 M87 GCMF with best fitting mass loss models. The left panel shows
models with powerlaw IMFS, and the right panel shows models with Burkert 82 Smith
(2000) IMFs. Fit parameters are given in table 6.1. The mass dependence for the
models are MO'0 (red solid line), M ‘0'25 (blue dot dashed line), and M “0'38 (green

dashed line).

6.3.1 COMPARISON TO MASS LOSS MODELS

We can compare these best fitting values to the expected values of 7' from the mass

loss coefficients in each of the models. For the Fall 82 Zhang (2001) case, this is easy,

as 7' is simply equal to pev - tfifetime, or
TFZ = 26966 (GPII/2 m In Atlifetime (6-20)

Taking the lifetime to be 12 Gyr, {e = 0.045, m = 0.7MQ, lnA = 12, and using the

value of p ~ 3.7%? gives a value of TF Z = 157200, which is close to the the best fit

95

value.

For the other two models, 7' does not directly relate to the mass lost, but to some
power of this quantity. Given the general form of the mass loss (equation 6.6), we can
note that for a model with mass dependence M '7, 7' will have units of 5%. Regardless,
the expected coefficient of such a model can be calculated from the dissolution time,
by noting that all of the cluster’s mass is gone by that time, such that tdiss = Ail.

The dissolution time for the Baumgardt 8c Makino (2003) model is given by equa-

tion 6.11, so

7 = L (mlnA)0'75 E —1 __3’_ 1 M29. (6 21)
BM 1.91 lkpc 220km/s 1 — e 1 x 106yr '

Using the same quantities as above, and assuming a circular orbit where e = 0 and

 

G M 1/2 , , M075 . .
v = (fim) ~ 488km/s, this gives TB M = 12778——®—yr , also conSIStent w1th the
best fitting value.
Repeating this same procedure for the Lamers et al. (2005) model dissolution time

yields

_ 0.62 p 1/2 t. .
TLamers = 0871122 (104MO) (W) (fifi) (6°22)

where we take Cam, = 810. This yields a value of TLameTs = 8605.9 that is much
greater than the best fit values. To match the best fit value, we must change Gem, ~
1847.9, which greatly increases the lifetimes for these clusters compared to the Lamers
et a1. (2005) prediction, a result that is not surprising, given that Gem, is based on

samples that likely contain objects that are not fully bound.

6.3.2 MASS TO LIGHT RATIO

As the cluster evolves, mass segregation moves the most massive stars to the core, and
the lighter stars to the outskirts of the cluster. These light stars are then preferentially

stripped from the cluster, leaving the fraction of more massive stars in the cluster to

96

rise over time. The lightest stars contribute only slightly to the mass of the cluster,
but contribute even less to the cluster luminosity. Therefore, as these stars leave the
cluster, the mass to light ratio drops.

Generally, the mass to light ratio changes more dramatically due to the evolution
of the stars in the cluster. However, for old GC populations such as M87, this effect
can be assumed to be negligible, as the stars with short evolution times will have
all evolved. The N-body simulations of Baumgardt 8c Makino (2003) track the effect
the loss of individual stars has on T, as well as the changes due to stellar evolution.
They correct for the stellar evolution, giving just the changes due to dynamics. For

simplicity, we model these changes as

. t
AT ( t ) = 00’ {dissolution < 02 (623)
tdissolution 1 _ 5 t ' t > 0 2
6 6 tdissolution ’ tdissolution .

 

with the initial mass to light ratio T0 = 3.0.
As it is the best fitting model, we use the c) = 0.0, Burkert 83 Smith (2000) IMF
model to examine the effects of this changing T. This also has the benefit that the

fraction of life is easy to calculate

t M
_ = 1 _
tdissolution M + T

 

(6.24)

due to the linear nature of the mass loss. This allows the individual cluster masses
to be corrected for the changes in T. As this depends on the fit value of 7', the
correction and the fit must be iterated until a stable solution is found. The final
fit gives a smaller value for 7', but with a better value of x2. This new best fitting
7" also matches much better with TF2, suggesting that the changes in mass to light
are significant to the low mass end of the GCMF. Table 6.2 presents the best fitting

model for the data using this variable mass to light ratio, and the fit is shown in

97

figure 6.3.

 

 

 

 

 

3 I l ,a T
'
l
'l
.’
.I
.I
.I
.I
.I
I
f
I"' i
0"
2 - _
A
91
v
o
F.
b0
O
l—.
1 r _
...... 6‘
~““
“‘
n x
‘\
‘\
‘\
§
\\ .J
‘S
‘\
‘\
\
“‘
“
‘\
~‘
0 . 1 . l L 1 ---------
3 4 5 7
10810(M/Mo)

Figure 6.3 Best fitting mass loss model, incorporating the variable mass to light ratio
from Baumgardt & Makino (2003). The dashed line shows the calculated t/tdiss for
the mass, and the dot dashed line shows the mass to light ratio used at each mass.

98

Table 6.2: Best fitting mass loss fits incorporating variable mass to light ratio.

 

 

 

Dataset Burkert 85 Smith (2000) log10(MC) '7 7' A X2

ACS 10g10(Mc) = 6.419 0.00 149604 2.311 0.7587
WFPC/2 log10(MC) = 6.350 0.00 155177 1.996 0.7979
Harris (1996) 10g10(Mc) = 6.046 0.00 113297 0.936 0.9486

 

6.3.3 COMPARISON TO PREVIOUS RESULTS

Figure 6.4 shows how this new mass function compares to the one presented by Waters
et a1. (2006). The extra depth the new sample provides a tighter constraint on the
mass loss rate, as the smallest clusters have more influence on the mass loss rate. The
effects of the changing mass to light ratio also produce the effect predicted by Waters
et a1. (2006), that the premature loss of clusters was an artifact of such changes.
When this new variable mass to light ratio method is used to fit the Waters et al.
(2006) data, the results match up very well. The differences between the ACS and
WFPC / 2 mass loss rates are within 5% of each other.

We can also compare our results to those presented by Jordan et al. (2007) based
on the results of the ACS Virgo Cluster Survey. They fit “evolved Schecter functions”
to their data, similar to what is done here and in Waters et a1. (2006) with our
Burkert & Smith (2000) IMFs. Unfortunately, they do these fits in magnitude, forcing
a conversion to our mass based units. For M87, Jordan et a1. (2007) present the

following values:

('1‘) = 2.67 (6.25)
5 = —7.287:l:0.089 (6.26)
m0 = —9.850:l:0.232 (6.27)

99

where we can convert these values to our units as

5 = C — 2.510g107' (6.28)
mc = C — 2.516g10 MC (6.29)
C = V + 2.516g10 M = 5.928 (6.30)

This gives values TACSVCS = 1976612 and loglo M0 = 6.321.

This neglects the fact that the F435W filter used by the ACSVCS does not directly
match the profile of a V filter, but based on the reasonable agreement in the cutoff
mass, this appears to work well regardless. The high value for TACSVCS derived
in this manner, along with the truncation of their published luminosity functions
indicates that they suffer from the same overestimate of the decay rate as the original
Waters et al. (2006) data. In addition, as shown above, neglecting the dynamical
effects on the mass to light ratio also tends to increase the best fit values of 1'.

The Harris (1996) Milky Way sample was also fit, although the number of clusters
used is far smaller. It is interesting to note that the best fitting cutoff mass is

significantly lower in the Milky Way.

6.3.4 COLOR DEPENDENCE OF THE MASS Loss

The analysis of the red and blue GCLFS indicated that the distributions were signif-
icantly different than each other. We can investigate this further by examining the
best fit mass loss models for these samples. Due to the splitting of the sample, we
are restricted in how small a mass we can reliably fit. For both samples, the clusters
below loglo M ~ 4.5 were excluded, as the evolution of T depletes the bins in such a
way as to bias the fits to large mass loss rates. Figure 6.5 shows the GCMFs for the
two samples, along with their best fitting mass loss models.

This plot clearly shows that the blue sample has a higher cutoff mass than the

100

 

10810 (‘1’)
N

p—A

 

 

 

 

.' ,- A l \ \
g, . . ( .
. » , ~ ~ ~.
.7 ,- . 1 . \ .
.' , _, \ v, -
_. ,9 \ .
. ' ’ ,‘ L ‘ '- '
1 x \ '.
X , \ -
’ /' - \
z, ' \ l '
r/ _/ \ \ -‘
‘ / 'v . \ .' _
.' ‘ \ \-
, \ -
,‘ K \'.
,‘ q 1
/ l l
«A ‘ \
0 n l l 1 l 1 l n
' 6

5
10510 (M / Mo)

Figure 6.4 Comparison of this GCMF (red) to that presented by Waters et al. (2006)
(green) . The mass functions are scaled to provide the same normalization, and the
best fitting mass loss models are shown for both. Note how the added depth in the
new data offers better constraints on the mass loss. In addition, the best fit model
from Jordan et al. (2007) is plotted for this new ACS sample (pink dotted), and the
Milky Way data of Harris (1996) along with the best fit mass loss model (blue).
red sample (log MC blue = 6.45, log MC red = 6.32). This is a rather significant shift,
corresponding to a 30% difference in the masses for the clusters at the turnover.
Such a large difference may be explained by a simple difference in the mass to light
ratios due to the metallicity differences, as such changes can range up to around 10%
(Ashman et al., 1995).

As the metal poor clusters are believed to be older, any evolutionary effect should
be more pronounced on those clusters. Without such a signature, it seems unlikely
that some high mass evolutionary effect can explain this difference.

The mass loss rates found by the best fitting models show that the blue clusters

do seem to have lost more mass than the red (75),“, = 205113, Tred = 130304). This

101

is to be expected, as stellar evolution is linear with time, so the older blue clusters
will have lost more of their mass over time. However, directly taking these values of
7' suggests that the blue population is nearly twice as old as the red population. This
is unreasonable, and suggests that the fit values may be biased due to the lack of the
lowest mass clusters.

However, if we shift the cutoff masses to be equivalent (in other words, assuming
that the metallicity dependence of T explains the effect fully), we can determine the
new mass loss for the models. Holding loglo MC red constant, we find a new fit to the
blue clusters where 71,1119 = 152052. If we assume that the mass loss is the same for

both types of cluster, this gives a ratio for the ages of the red and blue clusters as

t . .
hfetlmered N 0.87 (631)
tlifetime blue

For the nominal age of 12 Gyr for the blue clusters, this suggests that the red clusters

formed about 1.5 Gyr later, consistent with expectations.

102

 

10310 (‘1’)

 

 

P—
N]

 

 

0‘.—

0 4 - .
10810 (M/Mo)

Figure 6.5 Comparison of the GCMFs for the red and blue samples. The blue sample
clearly has a higher cutoff mass, and as it turns over faster, a larger amount of mass

loss.

103

 

CHAPTER 7:
STRUCTURAL PARAMETER FITTING

 

The superior resolution offered by our data Opens up the possibility of looking at more
than simply the brightness of the clusters. The fact that the cluster widths are signif-
icantly wider than the local PSFS indicates that we are resolving the actual structure
of the cluster. By generating simulated clusters, and matching these templates to the
data, we can estimate the parameters of that real cluster. The program superking
(see Appendix A.2) was written to perform the generation and fitting of simulated
clusters. Briefly, superking reads in a list of extracted GC images and associated
PSFS, estimates the initial parameters from the data image, and then uses a variety
of fitting stages to calculate the set of parameters that yield the best X2 value.

As discussed in chapter 2, only three parameters are generally needed to fully
define a King model: the concentration, any of the defined radii, and a brightness
calibration. We have chosen the tidal radius and the total cluster flux to define our
models, but due to the relations between the model parameters, any equivalent choice
would be equally valid (say, the core radius and the central surface brightness). In
addition to the parameters that define the model, there are parameters related to how
the cluster projects onto the image plane. As the sampling of the model can change
significantly with small changes in the center of the cluster, we must include the peak
position of the cluster, 9:0 and yo, as parameters in the fits. This must be fit as part
of the template, as asymmetries in the PSF can skew a simple centroid away from

the “true” value. Finally, since the local background is not guaranteed to be exactly

104

zero, a constant background level is also added to the parameters.

Other parameters can be added to this list, but, in general, are less important
to the fit quality. The PSF centering can be allowed to move, to more accurately
represent the effects of the pixel response function on subpixel shifts in the centering.
This was included in the fitting, but most PSF shifts are negligible, so this parameter
tends not to change much during fitting. The cluster ellipticity can also be included
by laying the model onto a template image as a function of semi-major axis instead of
a simple radius. However, this is not included in the King model formalism, and the
definition of the tidal radius is not well constrained for such an objects. As we have
required all cluster candidates to have at most B? = 1.5, we should have few highly
elliptical objects in the fitting sample, so this parameter can be safely excluded.

The data images are constructed by extracting a box of 128 x 128 pixels centered
on the cluster peak. These were drawn from the F814W image that had been drizzled
to twice the nominal ACS resolution, such that 1 pixel = 0(.’025. The choice was made
to use this instead of the single resolution data as the increase in resolution provides
more pixels on a given cluster, so even though the signal decreases, the quality of a
given fit is better constrained. The image size was chosen to correspond to a projected
radius of ~ 128 pc, which is larger than the tidal radii of most Milky Way globular

clusters. Therefore, all extracted images should fully contain the G0.

7.1 ERROR ANALYSIS

Even before running the fitting, it is clear that by fitting six main parameters by
minimizing the X2, there will be some degeneracy in the fits. This can cause coupled
errors between the best fitting values. The first source of error is between :20, yo, and
re. This arises from the fact that miscentering the cluster forces the fit to move light
away from the model peak to accommodate the peak in the data, leading to a model

with a higher core radius. A similar effect can be caused from errors in the PSF. If

105

the PSF used for fitting is wider than the true PSF, then the best fitting model will
have a smaller core radius to ensure that the final convolved simulated cluster has
the correct effective width.

The largest problem in the fitting couples 0, 7,3, cluster magnitude, and the back-
ground level together. This coupling is easiest to visualize by assuming an error in the
background level in which BFit > BTme. In this case, the best fitting tidal radius will
be smaller than the true value, as any wide tails will be swallowed by the increased
background level. Since the cluster is then smaller, fewer pixels are considered to
be part of the cluster, so the best fit flux is also smaller. As the size of the core is
changed little by this error, but the tidal radius has decreased, the fit concentration
will be lower as well. The opposite case works similarly, with a smaller background
creating clusters with larger rt, larger F, and higher c.

Given the set of expected errors, it is reasonable to wonder how well any of these
parameters can be fit. To test this, artificial globular clusters were generated for the
full range of expected parameters. These fake clusters have Poissonian noise for the
cluster light, as well as Gaussian noise with the value of a chosen from the range
observed in the data images. Looking at the quality of these fits and the variance in
the parameters provides an estimate of how well we can believe the fits to the real
data.

This estimate of the error does ignore two major sources of uncertainty. First, as
the template images are created using the same PSF as they are fit with, they are
certain to have better fits than any real data. The PSFS used are the best estimate
of the true image PSFS, but are not likely to be perfect matches. This uncertainty
will make the fits to the real data noisier in the core. The second problem is that the
templates are made with Gaussian distributed noise over a fixed fiat background. The
real images often contain other globular clusters, and may not have such an idealized

noise and background. This will cause issues with the background estimate, and as

106

discussed above, that can seriously alter the fit.

7 .1.1 POSITION ERRORS

Figure 7.1 shows the histograms of the position deviations from the simulated cluster
fitting, along with the best fitting Gaussians to the distribution. It is clear that the
simple position errors have wide wings that are poorly fit by the Gaussian model.
However, incorporating the PSF offsets largely eliminates these wings. This points
to the PSF offset changing the modeling only Slightly, and working more as a second
component of the positioning. If we compare the “effective centers,” defined as 930 —
aspsp, between the input and output values, we can see that there are no offsets
between the two (Ax = —0.0005, Ay = —0.0007) and very small scatter (ox =
0y = 0.016). This confirms that the centers found during the fitting process are

exceptionally accurate.

7.1.2 MAGNITUDE AND BACKGROUND ERRORS

The error in the measured magnitudes is in general small (0M = 0.004). For the
faintest clusters, however (V > 25), there is an extended tail where the measured
magnitude is even fainter. This suggests that the fitting routine fails to correctly
identify and fit these clusters, and provides the limit for the faintest cluster that can
be reliably fit.

The background is generally well measured, although the distribution is wider
than a simple Gaussian. Figure 7.2 shows the scatter in the background as a function
of the input magnitude, with the largest scatter coming from the brightest clusters.
This is to be expected, as these bright clusters will tend to skew the statistics used
for calculating the background. In addition, the objects in this sample with the
largest scatter also are the clusters with the largest input tidal radii, which would

also naturally lead to errors in the background, as these large clusters will serve to

107

400 . l r 1 T l . 100 T l r r . l

 

 

300

l
l

750 ‘ 4

T

Z200 2500

100* * 250* *
0 MMW. ___I _l__ J .

'- .0 —0.5 0.0 0.5 1.0 0R).

I
l

l
l

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(5 -0. 25 0. 25 0. 5

5'30 IN — 5'30 OUT $01N * $0 OUT)0- ($13 IN - 1‘13 OUT)
400 . I . I . I . 100C 1 l l
300* ~ 750* l? d

r)

2200* * 2500* *
100* * 250* _
0._Ql *‘J 4‘ ‘ i ‘ "‘ ml ‘ “ ‘

.0 -05 0.0 05‘“ 1.0 050.5 -025 0.0 0. 25 0.5
yOIN ‘ 1’10 OUT (90 IN * 90 OUT) * (yPIN- yPOUT)

Figure 7.1 Expected errors in position from fitting simulated clusters.

raise the image average.

7 .1.3 STRUCTURE ERRORS

Given that the background shows a problem in fitting the largest clusters, it is rea-
sonable to investigate how well the structure of such large clusters is fit. The central
potential fits for these clusters are very good, but show a large scatter and offset in
the best fitting Rt. A detailed check of the source of this scatter shows that the limit

of reliable Rt fits occurs at RtIN ~ §image size. This is an obvious limit, as larger

108

clusters are not fully contained on the fitting image, and so should not be expected
to be fit well.

The effects of the input radius on the central potential fits show almost the exact
opposite trend. The largest clusters have excellent fits to the central potential, due to
the large number of pixels devoted to the central core. The small clusters, however,
have a much larger scatter. This is reasonable, as such small clusters begin to all look
like a point source. If we consider a fairly diffuse cluster, with c = 1.0, but then set
Rt = 5 pixels, then the central pixel will contain 3% of the total light of the cluster.
This effect is even worse as the concentration increases, making all small clusters
effectively point like. Using this as a guide, we can set a lower limit of Rt > 5 pixels

for reliable fits to the central potential.

7.2 FITTING RESULTS

From the sample of all clusters detected, 1194 were chosen for fitting. These clusters
were selected to be brighter than V = 25, and were chosen to be the brightest clusters
within 25 pixels. This spacing requirement was added to allow PSFS to be generated
without problems from overlapping. Once the best fitting model was calculated with
superking, the various derived parameters defined in section 2.4 were calculated as
well. Given the results from the fitting of the simulated data, all clusters that have
final tidal radii less than 5 pixels or greater than 64 pixels were excluded from further
analysis. This leaves a final sample of 1096 clusters.

Figure 7.4 shows one of the extracted cluster images, along with the residual
for the best fitting model. The deviations in the core are generally less than 5%,
consistent with the expected error in the PSF. Based on this fact, it seems unlikely
that any fit can be made with a better residual, without some new understanding of

the PSF.

109

800 - l I I I I I 600 I I - I I I

 

 

 

 

 

 

 

  

 

 

 

 

 

 

 

 

90.25 -0.125 0 0.125 0.25 90.25 -0.125 0 0125 0.25
mIN * mOUT mIN * mOUT
500 I r I I #71 r 10 I + I - I I I
- j )- i +
400— — f +
5 — I i —
I-I 3 i
_ D E i 2 i
O 'l E 3;
“i I - E 1 i I I —
— z 1 E 5 g
a: i I f l
_5 _ i + _
_ . 1T 1 I l 4 I .
0—2.5 -1.25 0 1.25 2.5 HIS 18.75 22.5 26.25 30
BIN * BOUT mIN

Figure 7.2 Expected errors in magnitude and background level based on fitting simu-
lated clusters. The upper right panel shows the improvement in the magnitude errors
when the clusters fainter than I ~ 25 are excluded.

110

 

 

 

 

 

 

 

 

 

T I
0.375 - 0.375 —
:6: :5
O I O
20.25 .- 620.25- —
\ -' \ ;,
Z L ES Z r.
0.125 5;; — 0.125 5'; -
1| . .
' g .' '.
l ;'. - f '.
I '-'- 3 E
n A l m.
0‘-‘2 5 —112 0 1 25 ' 2 5 0‘0-5 —2 5 0' 2 5 5
tr3V0 IN - 9V0 OUT RtIN * 9from

Figure 7.3 Expected errors in the structural parameters based on fitting simulated
clusters. The red solid line shows the results for R¢IN < 5, and the green dashed
line for Rth > 50. The blue dot dashed line is the result for the remaining clusters,
showing the small scatter for such objects.

 

Figure 7.4 Original extracted cluster (x,y) = (4562,6772) and residual for best fitting
model. The companion cluster near the top of the frame does not seem to skew the
fitting at all.

111

 

CHAPTER 8:
GLOBULAR CLUSTER STRUCTURE
RESULTS

 

Given this extensive sample of structural parameters, we can compare against the
Milky Way results of Djorgovski & Meylan (1994). The Milky Way is obviously the
best observed, and has a some advantages over any extragalactic sample. First, the
true three dimensional distance from the center of the galaxy can be easily found,
whereas only the projected distance is observed for M87. In addition, the small
distances to the Milky Way clusters allows Spectra to be taken in a reasonable amount
of time, providing a direct measure of the central velocity dispersion of the stars.
Finally, even the faintest clusters of the Milky Way have values for at least some of
their structural parameters, and as discussed in Chapter 7, the fitting in our data is
limited to V ~ 25.

Despite these issues, the M87 sample does have one main advantage over the Milky
Way. The shear number of globular clusters ensures that for the range of parameters
that can be fit, the results should be more statistically sound. Over the range in
luminosities that we consider, Djorgovski & Meylan (1994) used 116 clusters (out
of a total sample size of 143 clusters), giving this new sample nearly ten times the

number of objects.

112

8.1 EFFECT OF CLUSTER LUMINOSITY ON STRUCTURE

Figure 8.1 shows the relations between several structural parameters with the absolute
V magnitude. These are plotted long with the Milky Way data taken from Harris
(1996) , which is effectively the same set of measurements used by Djorgovski & Meylan
(1994). This figure shows that the two samples line up rather well, and illustrates
that brighter clusters tend to be more concentrated, and have smaller cores than lower
luminosity clusters.

These trends become more clear in figure 8.2, in which the M87 data is median
binned in magnitude, with the error bars calculated using the non-parametric method

presented by Djorgovski & Meylan (1994), where

0 = 0-7415(Q75 - Q25) (8-1)

where Q75 and Q25 are the 75th and 25th percentiles in the bin. This gives the same
value as the standard deviation if the bin is populated with Gaussian distributed
data, but in the case of data with many outliers, provides an estimate that ignores
such points.

The binned data shows a clear trend of increasing concentration with brighter
clusters. This trend continues up to My ~ —10, beyond which the concentration
flattens out, with a constant concentration level for the brighter objects. This flat-
tening may represent a failing in the quality of the fits for clusters with the smallest
core radii. However, the binned results of Djorgovski & Meylan (1994) also show such
a flattening above My ~ —9. This suggests that the effect may be real, and as the
M87 data extends about two magnitudes brighter, we are simply seeing more of this
flattening.

The enhancement in concentration also shows as a drop in the binned core radius

over the same luminosity range. Djorgovski & Meylan (1994) fit this drop as re ~

113

L;0.5i0.25, or in terms of magnitudes

log10 Tc oc 0.2MV (8.2)

Fitting this trend over the full range of the M87 data gives loglo rc oc 0.04MV, which
does not match well. However, limiting the sample to only those clusters fainter than

My ~ —9 provides a fit of
loglo rc = 0.15MV + 0.943 (8.3)

well within the error range given by Djorgovski & Meylan (1994)
As the central density increases with both increasing concentration and luminosity,
it is not surprising that the trend in this parameter is fairly strong. The Djorgovski

& Meylan (1994) fit gives p0 ~ Lgfl, or
loglo p0 oc —0.8MV (8.4)

Restricting the fitting again to the range of magnitudes considered for the Milky Way,

we find a trend of

16g10 p0 = —0.69MV — 0.939 (8.5)

for the M87 data, which again matches.

The lack of a trend between the half light radius with luminosity is a well known
result for globular clusters, due to the fact that it contradicts what is observed for
elliptical galaxies, which have larger eflective radii with increasing mass. This con-
stant half light radius has been observed in young cluster systems (Zepf et al., 1999;

Larsen, 2004), which suggests that the half light radius is set during formation. The

114

M87 data provided a confirmation of this lack of trend, with

1ch10 n. = 0.013MV + 0.604 (8.6)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 8.1 Parameter correlations with absolute magnitude. The red crosses are the
new M87 data, and the green asterisks are the Harris (1996) data. All distances are
in parsecs, and the central density is given in solar masses per cubic parsec.

115

8.2 EFFECT OF DISTANCE ON STRUCTURE

The structure of globular clusters is expected to change with the distance from the
center of the galaxies. The clusters that are close to the core will shrink due to the
higher potential felt by the cluster, which in turn raises the concentration and central
density. However, due to the projection eflects, any such trend is washed out in this
sample.

Figure 8.3 shows the trends with projected distance, although no attempt was
made to try to fit trends to the data. Qualitatively, the clusters with the lowest
concentration do seem to be at the largest projected distance. Similarly, the half
light radii increases with distance. Unfortunately, without a good deprojection, no

firm trends can be established.

8.3 RELATIONS BETWEEN CORE PARAMETERS

Djorgovski & Meylan (1994) also Show the correlations between various core structural
parameters, plotted here in figure 8.4. The fairly tight correlations are believed to
be the result of clusters evolving toward core collapse. The central density of such a
cluster is expected to follow the relation p0 ~ rc'2'23 (Binney & Tremaine, 1987), the
same result found by Lynden-Bell 88 Eggleton (1980) by modeling core collapse in a

thermally conducting gas. A fit to our data yields
p0 o< 7‘6—2‘12 (8.7)

remarkably close to the theoretical expectation.

116

8.4 THE FUNDAMENTAL PLANE OF GLOBULAR CLUS-

TERS

In addition to simple correlations between parameters, Djorgovski & Meylan (1994)
noted that the parameters were related by a manifold of only three dimensions: a
size, a shape, and a brightness. This is reasonable to expect, as globular clusters are
reasonably well fit by King models, which have exactly the same parameters. Beyond
this, Djorgovski (1995) determined that the parameters can be fit to “Fundamental
Plane.”

This plane is given in terms of the cluster surface brightness as

140 = (—4.9 :l: 0.2) (loglo a — 0.45l0g10 rc) + (20.45 :I: 0.2) (8.8)

(11);, = (—4.1 :l: 0.2) (loglo a — 0.7llog10 rc) + (19.8 :I: 0.1) (8.9)

for the core and half light radii. Djorgovski (1995) notes that the core result can be
transformed to match the virial theorem (1°C ~ 0'2I0-1), and that the half light radius
result can be transformed to match the fundamental plane result for elliptical galaxies
(Th N 01410—08).

Unfortunately, this derivation for a fundamental plane requires information that

we do not have for M87. However, if we accept that the clusters are fit by King

models, then we can take advantage of a “cheat.” We have defined the core radius as

902
= 8.10
Tc V 47TG'p0 ( )

This a is formally the 3-d velocity dispersion, and not the line of sight dispersion mea—

 

sured by spectra and used for the fundamental plane relations. However, McLaughlin
(2000) shows that for concentrations greater than 0 ~ 1, these two values deviate by

less than 10% from each other. Therefore, we can trade the central density to get a

117

velocity dispersion. Doing this, and then fitting the best plane to the M87 data yields

110 = —5.0 (loglo a — 0.53log10 Q) + 21.4 (8.11)

04)), = —5.0 (loglo a — 0.6010g10 TC) + 21.16 (8.12)

which are very close to the cited fundamental plane relations.

A closer look at the cheat indicates that our calculation of po is based directly on
values of 110 and re (equation 2.53). Therefore, if simply rearranging the structural
parameters can provide a “fundamental plane,” for this data set that closely matches
those published, then it suggests that the fundamental plane is not a surprising effect,
but merely the consequence of the fact that globular clusters are well defined by King

models.

8.5 LMXB PROBABILITY

Low mass X-ray binaries (LMXBS) are binary star systems in which one star overfills
its Roche lobe, dumping mass onto the companion star. This companion is assumed
to be a compact object, like a neutron star or black hole, as the X-ray luminosity is
so great. The source of these X—rays is believed to be gas that is heated as it spirals
into the companion gravity well. These LMXBS have been found with surprisingly
high numbers in globular clusters, leading to the theory that these objects are formed
when two previously unrelated stars form a close binary due to stellar interactions.
Such an event would be vastly more likely to occur in the dense stellar environment
of a globular cluster than in the field.

To investigate these objects in the M87 sample, we used the catalog of X-ray
sources in Jordan et al. (2004) and matched the objects against our catalog of glob-
ular clusters. Figure 8.5 shows the structural parameters and CMD for the objects

matched, along with the general cluster population. As has been noted by many au-

118

thors before (Jordan et al., 2004; Sivakoff et al., 2007; Kundu et al., 2007 ), the bright
clusters and the more metal rich (or redder) clusters are more likely to host LMXBS.

Furthermore, the structural parameter data suggests that for the GCs that host
LMXBS that are not among the brightest clusters, those clusters with high concen-
trations (or equivalently, small core radii and high central densities) are the most
likely. This is consistent with the LMXB probability being a function of the stellar

encounter rate

2 3 2 3
r r
I‘ oc 8% ~ fl00_0_ ~ p657"? (8.13)

As shown above (section 8.4), both of these parameters tend to scale with the cluster
luminosity, so this effect is not surprising. However, by looking at the plot of I‘,
we can see that at a given magnitude, clusters with LMXBS tend to have the most

extreme values of I‘.

119

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 8.2 Parameter trends in averaged bins of luminosity. The solid green lines
show the best fit trends from Djorgovski & Meylan (1994), and the dashed blue lines
the new M87 fits.

120

 

 

 

 

1.3 —0.3 ‘ 0.02 ‘ 0.34 ' 0.66 ' 0.98 '

 

 

1081036 (kpc)
. , . . .

   

 

 

 

 

 

 

.98 ' 1.3 "30.3 ' 0.02 ' 0.34 . 0.66
10810R0(kPC)

A l
-0.3 0.02

' 0.34 ' 0.66 0
logioRc (kPC)

Figure 8.3 Parameter correlations with projected galactocentric distance.

121

‘ 0.98 ‘

 

 

 

 

 

 

10g 10 trc
00
l

K)
I

 

 

 

 

 

 

 

‘0l1‘224.20l16112
logloRc “V(0)

Figure 8.4 Correlations between core structural parameters. The blue line shows the
trend of central density with core radius.

122

 

 

 

 

 

 

 

 

I l I I I I 1 m I
-6 —8 —10 -12 -6 -8 -10 -12

 

 

 

 

 

 

 

 

 

- I ! 1T1; I '- 7’ I I 1 I I
%.5 0.75 1 1.25 1.5 3-6 -8 -10 -12

Figure 8.5 Plot of the locations of LMXBS with various structural parameters. The
small dots show the general cluster population, with the asterisks denoting clusters
that host LMXBs.

123

 

CHAPTER 9:
CONCLUSIONS

 

The data used for this thesis have shown that very deep images of galaxies can provide
a wonderful Opportunity to investigate the evolution and structure of globular clusters.
Although there are likely few datasets that have the same depth as the one used for
this project, other galaxies that are closer can easily be adapted to use the same
methods. However, it is clear that space based observing provides a tremendous
advantage over ground based imaging, due to the enhanced resolution.

The globular cluster luminosity function observed for M87 is consistent with pre-
vious space based surveys. The dispersion measured for this new sample does tend to
be larger than those found previously, likely due to the truncation of the GCLF in that
shallower and less complete data. Although the bimodality of cluster colors has been
observed before, this sample shows no evidence for a mass-metallicity relation that
some recent studies have found. Since the data presented here is so much deeper and
better sampled, it seems likely that and the appearance of such a relation is simply
the effect of photometric errors. By separating the cluster luminosity function into
radial bins, it is clear that the bins closer to the core of the galaxy have brighter mean
values and smaller dispersions. This result is entirely consistent with the expectation
of enhanced mass loss from clusters in those regions.

The effect mass loss has on the shape of the cluster luminosity function becomes
even more clear after considering the cluster mass function. The observed mass

function is fit very well by considering mass loss only due to the effects of evaporation

124

from two-body relaxation. By accounting for the changes in the cluster mass to
light ratio from the preferential loss of the lowest mass stars, the observed mass loss
rate is consistent with theoretical predictions. The mass function also explains the
differences that are seen between the GCLFS for the red and blue cluster samples.
The shift in the blue peak to brighter luminosities is due to a lower mass to light ratio
for these clusters compared to the red sample. After accounting for this, the peak
masses are effectively the same between both samples. Assuming that both types of
cluster lose mass via the same mechanism, then the smaller observed dispersion in
the blue GCLF is likely due to those clusters being older and having lost more of
their total mass due to that difl'erence in age.

The high resolution images used for this project show that reliable King model
structural parameters can be measured, even at the distance of the Virgo cluster. The
trends between these parameters appear to be universal, as the relations between the
cluster structure and luminosity observed for the Milky Way match well with those
found for M87. However, care must be taken when finding relations between the
various parameters. As globular clusters are well fit by King models, the observed
parameters are coupled together by the definition of such a model. Specifically, the
proposed “fundamental plane” for globular clusters appears to simply be an artifact
of this parameter coupling.

As has been observed before, low mass X-ray binaries are more likely to be found
in clusters that are brighter and more metal-rich. As these objects are believed to
form due to the capture of a star to a compact object to create a close binary, they
probability of finding an LMXB should correlate with the rate of stellar interactions,
I‘. This interaction correlates strongly with luminosity, obscuring which parameter is
most important. However, for the fainter clusters, only those clusters with larger than
normal values of I‘ are found to host an LMXB, suggesting that stellar interactions

are the important factor.

125

APPENDICES

126

 

APPENDIX A:
DETAILED SOURCE CODE

 

A.1 LIBKINGMODEL

LIBKINGMODEL is a C library written as part of this thesis to simplify the cre-
ation and evaluation of King (1966) models. In addition to creating the one di-
mensional volume and surface densities, it can also project this onto a two dimen-
sional image, and perform the convolution with an instrumental PSF. The library is
linked against the Gnu Scientific Library (GSL, http://www.gnu.org/software/gsl/),
the FFTW Fast Fourier Transform library (http://www.fftw.org/), and CFITSIO
(http: / / heasarc. gsfc.nasa. gov / docs / software / fitsio / ) to handle many numerical meth-
ods, and to allow standard fits images to be used.

In general LIBKINGMODEL is very fast, able to generate a complete template image
from scratch in about a second on a 1 GHz processor. Individual model generation
is faster than this, although the dominant component of the execution time remains
in solving the King model differential equation.

There are two main structures used in the library. The kingmodel, which contains

the 1D model, and the kingmodel_template, which holds the 2D image information.

 

The kingmodel internally calculates the main defining parameters of the model: the
central potential (W0), the core (Tc) and tidal radii (rt), and the King concentration
(c). All of these are in an undefined “internal unit” system, as are the volume and
surface densities. Since these parameters are generally normalized after the fact, this

is not major problem.

127

A.1.1 KINGMODEL.C

The functions defined in this file handle the allocation, generation, and evaluation of

the one dimensional King model.

kingmodel It km_model_alloc(int size); This function allocates a kingmodel struc-

 

ture pointer for evaluation. The size parameter denotes how many radial sam-

ples are used for the model.

void km_model_free (kingmodel *M) ; This destroys a kingmodel pointer, and frees all

 

memory associated with it.

void 1nn_model_solve(kingmodel *M, double W0); This solves the King model with

 

specified central potential @, and stores the calculated volume density. The

differential equation for the model density comes from Poisson’s equation, where

d 2dW _ 2 W 4 2
dr (1" d7 ) — 47707 p1(e erfi/W 77W (1 + 3W)) (A.1)

as the density can be written in terms of the reduced potential as

p(W) = p1 (8W erf\/W — :W (1 + gW)) (A.2)

This second order diflerential equation can be solved by breaking it into two

first order difl'erential equations

dW

217 = Y (A.3)
dY 2
fl — ~47er(W) — ;Y (A.4)

This system of differential equations is solved using a very accurate 8th or-

der Runge-Kutta Prince-Dormand solver from GSL. This method requires the

128

Jacobian of the system, which is easy enough to find:

0 1

J = —47rG (6W erf (x/W) — %W) -%

(A.5)

' . . _ 9
The first step to solve the model IS to calculate the core radlus, rc — , “W'
The next step is to find the King tidal radius, where W becomes zero. Unfor-

tunately, there is no analytic solution for this. Instead, it is found by using a
binary search algorithm to bisect the range of values considered until the range
becomes very small. The final range is then the uncertainty in the tidal radius,
and is restricted to be smaller than 1 x 108-re. Although this is a fairly stringent
constraint, the speed of the binary search algorithm allows this to be found with

only a dozen or so iterations.

Once the tidal radius is found, the radius is divided into samples. The first five

points are fixed to have the values

r1: {00,1 x 10’1°.1 x 10‘8, 1 X 10'6I1 X 1041} (A6)
t

to ensure that the very core of the cluster is reasonably sampled. The remaining

points are distributed logarithmically, such that

10810 (Ti) = 10810 (Tt) X i/_Size (A-7)

With the radius array filled, the volume density is then evaluated at each sample

radius, and the values stored as well.

Once both arrays have been filled, GSL is used to construct a cubic spline
through all the data points. This ensures that the evaluation can be performed

quickly at all radii, and that no further calls to solve the differential equation

129

are needed. This Spline is also stored in the kingmodel. Finally, the tidal radius
is solved once more, to determine the cutoff for the spline. This stage usually
only corrects the tidal radius slightly, and is done to check the limit beyond
which the spline cannot be evaluated anymore. This new tidal radius is then

used to calculate the model concentration.

void km_model_project (kingmodel *M); This function projects the volume density

 

into a surface density. GSL is again used to perform this integration. The

standard integral for projection is

 

E (R) = 2 / rt p(7‘)—-—T——d7° (A.8)

but, by noting that we can write f = p(r) and gfi = fl, we can use

integration by parts to remove the Singularity from the integral:

2 (R) = 2 ((p(I~)I/II2 — R2)” — f” “p—(7‘)I/I~2 — R265) (A.9)

R R (17'

which then allows a faster and simpler integration method to be used that does
not need to work around the point where r = R. This can be simplified further,
as we can notice that the first term is zero at both endpoints. This leaves the
projection integral in terms of the volume density derivative, which is easy to
calculate with the spline implementation. This integral is evaluated on the same

set of radial points as the volume density, and a Spline is created for it as well.

king-model It km_model_simple(double W0); For quick cases in which only a single

 

model is needed, this function allows a kingmodel to be allocated, and the

model with central potential fig solved and projected in one step.

double km_model_eva1_rho(kingmodel *M, double r);

 

double km_model-eva1_drho(kingmodel *M, double r);

 

130

double km_model_eva1_surf(kingmodel *M, double r);

 

double km-mode1_eva1_dsurf(kingmodel #11, double r); These four functions evalu-

 

ate the King model at a specified radius, and return the value of the p(r), (1%;2,
2(R), or d—figfl. The radius is again in the internal radius units, and so needs

to be converted using the stored king-model tidal or core radius.

A.1.2 KINGIMAGE.C

The functions in this file handle the creation of template images that contain a two
dimensional representation for a kingmodel. The layout parameters and image data

are stored in the kingmodeLtemplate structure.

 

kingmodeLtemplate *km_temp1ate_a11oc(int size) ; This allocates space for the im-

 

age data on a square image of length size.

void km_temp1ate-init (kingmodel_template *T,

 

double x0,double y0,doub1e px1_tidal_radius); This function sets the center

 

and tidal radius scale for use in the filling stage.

void km_temp1ate_fi11(kingmodel *M,kingmodel_temp1ate *T,char *method);
This fills the template image with the specified kingmodel surface density data.
The method can be either Sampled or Elliptical, although only the Sampled

method has been thoroughly tested.

The Sampled method functions by first assuming that the model has circular

symmetry. It then evaluates the projected radius as
2 . 2 . 2
R = (a: —1nt(x(_))) + (y —1nt(fl)) (A.10)

This ensures that the peak of the cluster is located on a pixel corner. The phase

values for the template, phi_x = 39 — int(x_0_) and phi_y = fl — int(y0), are

131

then stored to be applied following convolution. This is a necessary step, as
simply laying the pixels out directly will tend to underestimate the core value

for highly concentrated models after convolution.

The ’Elliptical’ method works in a similar fashion, using a radius equation

 

2 2
I = (6-5-4.466New» I

((—(2 — x9) adage) + (y — fl) cos(t_he_’g_a_)) /§) (AM)

where theta is the position angle of the major axis, and g is the ellipticity. This
method is not recommended, as it does not ensure the core is properly sampled,
nor does it correctly model the cluster as an elliptical object, but merely skews

the standard spherical form.

void km_temp1ate-normalize (kingmodeLtemplate *T,doub1e flux); This function

 

normalizes the image in a kingmodeLtemplate to have total counts flux.

 

void km_template_free(kingmodel_temp1ate *T);

 

This destroys the kingmodeLtemplate, and frees all associated memory.

 

A.1.3 KINGCONVOLVEC

As this library is designed largely to help fit King models to observed data, the gen-
erated images must be convolved with the PSF of the detector. This file contains the
functions needed to perform this convolution using the FFTW fast Fourier transform

library.

void km_template_plan(kingmodel_temp1ate *T);

 

To perform the Fourier transforms, FFTW requires “plans” to be generated to
Speed up the calculation. This function generates the plans necessary for a given

template. This only needs to be done once for a given template, as it depends

132

only on the template size.

void km-temp1ate_convolve(kingmodel-template *T, fitsimage *PSF);

 

This performs the actual convolution with Q. The PSF data is first padded
to the size of the template image, and then Shifted so that the peak value is
on pixel (0,0). The PSF data is wrapped around the edges, which ensures that

after the convolution, no unwanted shifts are introduced.

Following this, both the template and the PSF are transformed into the Fourier
domain, and convolved by multiplying the complex values pixel by pixel. This
product image is then transformed back to the real domain, yielding the con-
volved image. This image is then normalized to unit flux. If the values of
w and M are not zero, the convolved image is Shifted back to the correct

centering.

void km-template_shift(kingmodel_temp1ate *T); This function performs the shift

 

to remove the efl'ects of phi-x and phi_y using a bicubic interpolator.

A.1.4 KINGUTILS.C

This file contains miscellaneous routines, largely for converting between the various
model parameters. Since these parameters are functions of a given King model, they

are included here.

double km_concentration_to_wo(kingmodel #14, double c); The concentration for a

 

given W0 can be found by simply evaluating the model. This function works

the other way, searching for the value of W0 that has a concentration of 9.

double hn_integrate-f1ux(kingmodel *M, double R) ; This function finds the total

 

flux (in the internal units) enclosed within a given projected radius.

133

double km_find_rha1f (kingmodel *M) ; This finds the radius that contains half of the

 

total flux of the cluster.

void km_scale_radii(kingmodel *M, double rt, double *rc, double *rh);

 

This calculates the core and half light radii that correspond to the given

kingmodel, and have tidal radius g.

void km_calculate_surface_brightness(kingmodel *M, double mag,

 

double pixel_sca1e, double *muO, double *meanmu-half);

 

This calculates the central and half light surface brightnesses for a cluster of

given mag and pixel_sca1e in arcseconds.

 

void km_ca1cu1ate_centra1_density(kingmodel *M, double mag,

 

double pixel_sca1e,double distance, double *rhoO);

 

Calculates the central mass density for the model in solar masses per cubic

parsec.

A.2 SUPERKING

Superking is a program to determine the best fitting single mass King (1966) model
to a given globular cluster image. It automatically finds the best central potential,
tidal radius, center, cluster magnitude, and background level, with few required in-
puts.

The default way to run superking is to simply specify an extracted GC image and
a corresponding PSF. Alternatively, a list of GC and PSF images can be specified in a
file, and superking will proceed to fit all the images. The options supported by the

program are as follows:
—help, -h Display a list of the options available, and a usage summary.

—list, -L Read a list of GC and PSF image pairs from a file.

134

 

—gain, -g Specify the image gain for the GC images. The fitting is performed in

units of electrons, so the gain is required for the conversion.
—readnoise, -r Give the detector read noise to be included in the noise calculation.

—zeropoint, -z The zeropoint used for converting the measured flux into an instru-

mental magnitude.

—background, -K Fix the background level to be the Specified constant value for

the entire fitting.

—quick, -q Use only a single pass of the Levenberg-Marquardt solver to find the best

fitting model.

-slow, -s Run a more thorough search method that adds a second pass grid search

in W0 and RI.

—randomize, -R Randomize the initial model values by searching for a new maxi-

mum around the first guess.

I—brute, -B Change how many samples are taken during the grid search. The search

is square, so this value determines the number of samples on a side.
—iterations, -I Change the number of iterations used at a given sampling step.
—range, -J Multiplier for all the ranges considered when doing a search.
-N Set the radial resolution of the King model.

-Z Set the Size of the King model template image.

A.2.1 INITIALIZATION

Superking loads the GC and PSF images, and calculates a first set of estimates

for the fit parameters directly from the data. The uncertainty in the background

135

 

is found by calculating the standard deviation of the image, using an iterative 3o
clipping process. The cluster is assumed to be in the center of the image, so that
is set as the first guess. The image minimum is set to the background, although
that is nearly always an underestimate. As the central potential and tidal radius are
not easy to estimate directly from the image, they are set to 5 and 15 respectively.
These values are chosen simply because they are good “middle” values. The initial
magnitude is found by summing the image flux and converting it.

An error map is calculated and stored along with the GC image. This contains
the expected noise on each pixel, and is used in calculating the x2 value. The exact

form of the error map can be altered by changing options in the superking.h header

 

file. The default error calculation incorporates all of the major sources of error
ECU, y) = Gain ' IGC($’ 31“ + Ugeadnoise + Gain ° Ufiackground (A.12)

Once all the initial values are calculated, and the necessary images loaded, a logfile

is opened to store the results of the calculation for each step of the fitting process.

A.2.2 FIT EVALUATION

The quality of the fits is determined by calculating the x2 value for the model. This

is calculated as

 

 

24:? (Gain GCC’EI y)-(1'“111x-'1"(1”I%’)’13‘3“’1‘gr°‘1n‘i))2 (A.13)

=Nx'1Nyx E($)y)

For each new set of fit values, superking determines what has changed from the
previous run. As the various stages of the evaluation of a new model run at different
rates, it is worthwhile to attempt to prevent any unneeded calculations. Changes in
the photometric parameters, the magnitude and the background, require no changes

to the model, and are implemented as just a simple reevaluation of the x2 value.

136

Changes in how the cluster is laid out, either in the cluster centering or in the scaling
of the tidal radius, do require a new template image to be generated and convolved
with the PSF. However, the King model used previously can continue to be used in
all cases, except that of a change in W0.

Because of the different speeds involved in the changes of different parameters,
the speed of the fitting process can be improved by holding the slow parameters fixed
while the faster parameters are varied. This can help alleviate the degeneracy of the
fitting, as the best fitting magnitude and background can be quickly found for a given
set of W0 and Rt.

A.2.3 FITTING PROCEDURE

If the “—randomize” option is specified, the program first randomly samples magnitude
and background values in an attempt to find the best fitting initial values. The ranges
for these parameters are Am = :l:1, AB = 21:1000 counts. This is repeated for the
position, with a range in both 2:0 and yo of :l:1 pixel. This uses a randomized hill
climbing algorithm, which checks if the proposed value is better than the current
value, and if so, sets the new center to that location. This ensures that a good value
can be found quickly.

A grid search is the next stage of the fitting, with a range of W0 from 3.5 to
13, and a range of Rt from 0 pixels to half the image size. At each grid point, the
magnitude and background are searched using the same hill climb algorithm, over
ranges (Am, AB) = (:l:0.1, :l:150). This helps ensure that these parameters are the
best values possible for the given (W0, Rt). If this were not done, the degeneracy of
these four parameters would likely prevent equally good fits across the grid.

With the optional randomized refinement of the parameters finished, a Levenberg-
Marquardt solver (using the implementation provided by GSL) is used to attempt

to find a good set of fit parameters. This fit is generally reasonable, but the odd

137

 

degeneracy of the parameters can cause the Levenberg-Marquardt solver to often find
a local minimum instead of the final global best fit.

To rectify this, unless the “—quick” Option is given, superking repeats the hill
climb and grid searches again, with an optional second pass if the “—slow” Option is
also used. The ranges searched for this are given in table A.1. As is clear from these
ranges, the second pass is designed to refine the fit around the final best fitting model.

With the fitting completed, the best fitting model is calculated, and the results
saved to a series of images: the best fitting template, the difference image between
the data and the template, the x2 map, and the calculated error map. If the fitting
is being run from a list, a new file with the “.best” extension is created that saves the

best fitting values for each set of input images.

138

 

 

 

mam” modmu 3mm“ flour mdnm mdnm QNH mdfi m mmmm
one 8.3“ SSH made 39 is $3 I 5% E - m H mam
can.” :1“ 83H 93” QB“ cam” mummmwdfim m I o 2 I ......m qoBeEEovaam 33:

m condom EEO E Agnew BOG m E on Pa 9m ox: 88:8me

 

 

mafia 8m moms—Rm 50.8% JAN 2an

139

 

APPENDIX B:
GLOBULAR CLUSTER CATALOG

 

In order to allow easy access to the globular cluster data, a SQLITE database was
created to store the various measured parameters. Given that there are thousands of
objects, with up to a hundred parameters for each one, no other method of organizing
the data is practical.

Each Object stores all of the values measured by Source Extractor for both data
images. After matching the filters, the final RA and DEC are chosen, and the dis-
tance from the center of the galaxy is calulated. The next set of parameters are
the photometric data calculated during the calibration phases, with instrumental and
final magnitudes, and the aperture corrections and completeness values for both fil-
ters. These are combined to create the Object color and final completeness. Finally,
a quality flag is calculated to note how likely an object is to be a cluster. A quality
of zero is reserved for the good objects, with higher values being a bitmask to denote
all the reasons why an object is rejected.

For clusters that have matches in either the Kundu et al. (1999) or Waters et
al. (2006) samples, that data is stored as well, with a “merge” value that gives the
angular separation of Objects between catalogs. The Jordan et al. (2004) Xray data
is also matched in the same way, and stored for the clusters that match. Finally, the
raw output of superking is stored, along with the calculated structural parameters
defined in section 2.4.

This same database is also used to store the results from the false and UDF

140

frames, along with all the calibration data from the simulated cluster searches. These

components are included to allow for any reanalysis at a later date.

141

BIBLIOGRAPHY

142

REFERENCES

Anderson, J., 88 King, I. R. 2003, PASP, 115, 113

Anderson, J ., King, I. R. 2006, ACS Instrument Sci. Rep. 2006—01 (Baltimore: STSci)
Ashman, K. M., Bird, C. M., 88 Zepf, S. E. 1994, AJ, 108, 2348

Ashman, K. M., Conti, A., 88 Zepf, S. E. 1995, AJ, 110, 1164

Ashman, K. M., 88 Zepf, S. E. 1998, Globular Cluster Systems, (Cambridge: Cam-
bridge Univ. Press)

Baltz, E. A., Lauer, T. R., Zurek, D. R., Gondolo, P., Shara, M. M., Silk, J ., 88 Zepf,
S. E. 2004, ApJ, 610, 691

Barmby, P., 88 Huchra, J. P. 2001, AJ, 122, 2458
Baumgardt, H., 88 Makino, J. 2003, MNRAS, 340, 227
Beckwith, S. V. W., et al. 2006, AJ, 132, 1729

Bertin, E., 88 Arnouts, S. 1996, A88AS, 117, 393

Binney, J ., 88 Tremaine, S. 1987, Galactic Dynamics, (Princeton, NJ: Princeton Univ.
Press)

Brosche, P., Odenkirchen, M., 88 Geffert, M. 1999, New Astronomy, 4, 133
Burkert, A., 88 Smith, G. H. 2000, ApJ, 542, L95
Chemofl, D. F., 88 Weinberg, M. D. 1990, ApJ, 351, 121
COté, P., et al. 2004, ApJS, 153, 223

Da Costa, G. S., 88 Freeman, K. C. 1976, ApJ, 206, 128
Djorgovski, S., 88 Meylan, G. 1994, AJ, 108, 1292
Djorgovski, S. 1995, ApJ, 438, L29

Fall, S. M., 88 Zhang, Q. 2001, ApJ, 561, 751

Fruchter, A. S., 88 Hook, R. N. 2002, PASP, 114, 144
Gunn, J. E., 88 Griffin, R. F. 1979, AJ, 84, 752

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

143

Harris, W. E., Whitmore, B. C., Karakla, D., Okon, W., Baum, W. A., Hanes, D. A.,
88 Kavelaars, J. J. 2006, ApJ, 636, 90

Hénon, M. 1961, Annales d’AstrOphysique, 24, 369

Herschel, W. 1789, Philosophical Transactions Series I, 79, 212
Innanen, K. A., Harris, W. E., 88 Webbink, R. F. 1983, AJ, 88, 338
Jordan, A., et al. 2004, ApJ, 613, 279

Jordan, A., et al. 2005, ApJ, 634, 1002

Jordan, A., et al. 2007, ApJS, 171, 101

King, I. 1962, AJ, 67, 471

King, I. R. 1966, AJ, 71, 64

Koekemoer, A. M., Fruchter, A. 8., Hook, R., Hack, W., 2002, HST Calibration
Workshop, 337

Krist, J. 1993, Astronomical Data Analysis Software and Systems II, 52, 536
Kroupa, P. 2001, MNRAS, 322, 231

Kundu, A., Whitmore, B. C., Sparks, W. B., Macchetto, F. D., Zepf, S. E., 88 Ashman,
K. M. 1999, ApJ, 513, 733

Kundu, A., Maccarone, T. J., 88 Zepf, S. E. 2007, ApJ, 662, 525

Lamers, H. J. G. L. M., Anders, P., 88 de Grijs, R. 2006, A88A, 452, 131
Lamers, H. J. G. L. M., Gieles, M., 88 Portegies Zwart, S. F. 2005, A88A, 429, 173
Larsen, S. S. 2004, A88A, 416, 537

Lauer, T. R. 1999, PASP, 111, 1434

Lauer, T. R. 1999, PASP, 111, 227

Lynden-Bell, D., 88 Eggleton, P. P. 1980, MNRAS, 191, 483

Macri, L. M., et al. 1999, ApJ, 521, 155

McLaughlin, D. E., Harris, W. E., 88 Hanes, D. A. 1994, ApJ, 422, 486
McLaughlin, D. E. 2000, ApJ, 539, 618

Meurer. G. R., et al. 2002, HST Calibration Workshop (Baltimore: STSci)
Meylan, G., 88 Heggie, D. C. 1997, A88A Rev., 8, 1

144

 

Ostriker, J. P., Spitzer, L. J., 88 Chevalier, R. A. 1972, ApJ, 176, L51

Perlman, E. S., Harris, D. E., Biretta, J. A., Sparks, W. B., 88 Macchetto, F. D. 2003,
ApJ, 599, L65

Plummer, H. C. 1911, MNRAS, 71, 460

Schechter, P. 1976, ApJ, 203, 297

Schlegel, D. J., Finkbeiner, D. P., 88 Davis, M. 1998, ApJ, 500, 525
Schultz, H. 1866, Astronomische Nachrichten, 67, 1

Shapley, H., 1921, Bull. Nat. Res. Coun., 2, 171

Sirianni, M., et al. 2005, PASP, 117, 1049

Sivakoff, G. R., et al. 2007 , ApJ, 660, 1246

Spitzer, L. 1987, Dynamical Evolution of Globular Clusters, (Princeton, NJ: Princeton
Univ. Press)

Stoughton, C., et al. 2002, AJ, 123, 485

Strader, J., Brodie, J. P., Spitler, L., 88 Beasley, M. A. 2006, AJ, 132, 2333
Vesperini, E., Zepf, S. E., Kundu, A., 88 Ashman, K. M. 2003, ApJ, 593, 760
Waters, C. Z., 88 Zepf, S. E. 2005, ApJ, 624, 656

Waters, C. Z., Zepf, S. E., Lauer, T. R., Baltz, E. A., 88 Silk, J. 2006, ApJ, 650, 885

Whitmore, B. C., Sparks, W. B., Lucas, R. A., Macchetto, F. D., 88 Biretta, J. A.
1995, ApJ, 454, L73

Zepf, S. E., Ashman, K. M., English, J., Freeman, K. C., 88 Sharples, R. M. 1999,
AJ, 118, 752

Zhang, 6)., 88 Fall, S. M. 1999, ApJ, 527, L81

145

 

IIIIIIIIIIII)