Shaping the Globular Cluster Mass Function by Stellar-Dynamical
  Evaporation
THE ASTROPHYSICAL JOURNAL, 679:1272–1287, 2008 JUNE 1
Preprint typeset using LATEX style emulateapj v. 08/22/09
SHAPING THE GLOBULAR CLUSTER MASS FUNCTION BY STELLAR-DYNAMICAL EVAPORATION
DEAN E. MCLAUGHLIN1,2 AND S. MICHAEL FALL3,4
The Astrophysical Journal, 679:1272–1287, 2008 June 1
ABSTRACT
We show that the globular cluster mass function (GCMF) in the Milky Way depends on cluster half-mass
density, ρh, in the sense that the turnover mass MTO increases with ρh while the width of the GCMF decreases.
We argue that this is the expected signature of the slow erosion of a mass function that initially rose towards low
masses, predominantly through cluster evaporation driven by internal two-body relaxation. We find excellent
agreement between the observed GCMF—including its dependence on internal density rhoh, central concen-
tration c, and Galactocentric distance rgc—and a simple model in which the relaxation-driven mass-loss rates
of clusters are approximated by −dM/dt = µev ∝ ρ
h . In particular, we recover the well-known insensitivity
of MTO to rgc. This feature does not derive from a literal “universality” of the GCMF turnover mass, but rather
from a significant variation of MTO with ρh—the expected outcome of relaxation-driven cluster disruption—
plus significant scatter in ρh as a function of rgc. Our conclusions are the same if the evaporation rates are
assumed to depend instead on the mean volume or surface densities of clusters inside their tidal radii, as
µev ∝ ρ
t or µev ∝ Σ
t —alternative prescriptions that are physically motivated but involve cluster properties
(ρt and Σt) that are not as well defined or as readily observable as ρh. In all cases, the normalization of µev
required to fit the GCMF implies cluster lifetimes that are within the range of standard values (although falling
towards the low end of this range). Our analysis does not depend on any assumptions or information about
velocity anisotropy in the globular cluster system.
Subject headings: galaxies: star clusters—globular clusters: general
1. INTRODUCTION
The mass functions of star cluster systems provide an im-
portant point of reference for attempts to understand the con-
nection between old globular clusters (GCs) and the young
massive clusters that form in local starbursts and galaxy merg-
ers. When expressed as the number per unit logarithmic mass,
dN/d log M, the GC mass function (GCMF) is character-
ized by a peak, or turnover, at a mass MTO ≈ 1–2× 10
that is empirically very similar in most galaxies. By con-
trast, the mass functions of young clusters show no such fea-
ture but instead rise monotonically towards low masses over
the full observed range (106 M⊙ & M & 10
4 M⊙ in the best-
studied cases), in a way that is well described by a power law,
dN/d log M ∝ M1−β with β ≃ 2 (e.g., Zhang & Fall 1999).
At the same time, for high M > MTO, old GCMFs
closely resemble the mass functions of young clusters, and
of molecular clouds in the Milky Way and other galaxies
(Harris & Pudritz 1994; Elmegreen & Efremov 1997); and
it is well known that a number of dynamical processes
cause star clusters to lose mass and can lead to their com-
plete destruction as they orbit for a Hubble time in the
potential wells of their parent galaxies (e.g., Fall & Rees
1977; Caputo & Castellani 1984; Aguilar, Hut, & Ostriker
1988; Chernoff & Weinberg 1990; Gnedin & Ostriker 1997;
Murali & Weinberg 1997). It is therefore natural to ask
whether the peaks in GCMFs can be explained by the deple-
1 Dept. of Physics and Astronomy, University of Leicester, University
Road, Leicester, UK LE1 7RH
2 Permanent address: Astrophysics Group, Lennard-Jones Lab-
oratories, Keele University, Keele, Staffordshire, UK ST5 5BG;
dem@astro.keele.ac.uk
3 Institute for Advanced Study, Einstein Drive, Princeton, NJ 08450
4 Permanent address: Space Telescope Science Institute, 3700 San Martin
Drive, Baltimore, MD 21218; fall@stsci.edu
tion over many Gyr of globulars from initial mass distribu-
tions that were similar to those of young clusters below MTO
as well as above.
Our chief purpose in this paper is to establish and inter-
pret an aspect of the Galactic GCMF that appears fundamen-
tal but has gone largely unnoticed to date: dN/d log M has
a strong and systematic dependence on GC half-mass den-
sity, ρh ≡ 3M/8πr
h (rh being the cluster half-mass radius), in
the sense that the turnover mass MTO increases and the width
of the distribution decreases with increasing ρh. As observed
facts, these must be explained by any theory of the GCMF. We
argue here that they are an expected signature of slow dynam-
ical evolution from a mass function that initially increased
towards M < MTO, if the long-term mass loss from surviv-
ing GCs has been dominated by stellar escape due to internal,
two-body relaxation (which we refer to from now on as either
relaxation-driven evaporation or simply evaporation).
Fall & Zhang (2001; hereafter FZ01) explain in detail why
cluster evaporation dominates the long-term evolution of the
low-mass shape of observable GCMFs. Briefly, stellar evolu-
tion removes (through supernovae and winds) the same frac-
tion of mass from all clusters of a given age, and so cannot
change the shape of dN/d log M (unless special initial con-
ditions are invoked; cf. Vesperini & Zepf 2003). Meanwhile,
for GCs like those that have survived for a Hubble time in
the Milky Way, the mass loss from gravitational shocks dur-
ing disk crossings and bulge passages is generally less than
that due to evaporation for M < MTO (FZ01; Prieto & Gnedin
2006).5
As we discuss further in §2 below, the evaporation of tidally
5 It is possible that there existed a past population of GCs with low den-
sities or concentrations, or perhaps on extreme orbits, that were destroyed in
less than a Hubble time by shocks or stellar evolution. Our discussion does
not cover such clusters.
http://arxiv.org/abs/0704.0080v4
2 McLAUGHLIN & FALL
limited clusters proceeds at a rate, µev ≡ −dM/dt, that is ap-
proximately constant in time and primarily determined by
cluster density. FZ01 show that a constant mass-loss rate
leads to a power-law scaling dN/d log M ∝ M1−β with β → 0
(corresponding to a flat distribution of clusters per unit lin-
ear mass) at sufficiently low M < µevt in the evolved mass
function of coeval GCs that began with any nontrivial ini-
tial dN/d log M0.
6 To accommodate this when dN/d log M0
originally increased towards low masses as a power law, a
time-dependent peak must develop in the GCMF at a mass
of order MTO(t) ∼ µevt (FZ01). But then, since µev depends
fundamentally on cluster density, so too must MTO.
A β ≃ 0 power-law scaling below the turnover mass has
been confirmed directly in the GCMFs of the Milky Way
(FZ01) and the giant elliptical M87 (Waters et al. 2006), while
Jordán et al. (2007) show it to be consistent with dN/d log M
data for 89 Virgo Cluster galaxies, and it is apparent in deep
observations of some other GCMFs (e.g., in the Sombrero
galaxy, M104; Spitler et al. 2006). As regards the peak itself,
old GCs are observed (e.g., Jordán et al. 2005) to have rather
similar densities on average—and, therefore, similar typical
µev—in galaxies with widely different total luminosities and
Hubble types. (Inasmuch as cluster densities are set by tides,
this is probably related to the mild variation of mean galaxy
density with total luminosity; see FZ01, and also Jordán et al.
2007.) Thus, an evaporation-dominated evolutionary origin
for a turnover in the GCMF appears to be consistent with
the well-known fact that the mass scale MTO generally dif-
fers very little among galaxies (e.g., Harris 2001; Jordán et al.
2006).
If this picture is basically correct, it implies that, even
though MTO may appear nearly universal when considering
the global mass functions of entire GC systems, in fact the
GCMFs of subsamples of clusters with similar ages but dif-
ferent densities should have different turnovers. In §2, we
show—working for definiteness and relatively easy observ-
ability with the half-mass density, ρh—that this is the case for
globulars in the Milky Way. We fit the observed dN/d log M
for GCs in bins of different ρh with models assuming that
(1) the initial distribution increased as a β = 2 power law at
low masses and (2) the mass-loss rates of individual clusters
can be estimated from their half-mass densities by the rule
µev ∝ ρ
h . In §3 we discuss the validity of this prescription
for µev, which is certainly approximate but captures the main
physical dependence of relaxation-driven mass loss. In par-
ticular, we show that the alternative mass-loss laws µev ∝ ρ
and µev ∝ Σ
t —where ρt and Σt are the mean volume and
surface densities inside cluster tidal radii—lead to models for
the GCMF that are essentially indistinguishable from those
based on µev ∝ ρ
h . The normalization of µev required to fit
the observed GCMF implies cluster lifetimes that are within
a factor of ≈ 2 (perhaps slightly on the low side, if the ini-
tial power-law exponent at low masses was β = 2) of typical
values in theories and simulations of two-body relaxation in
tidally limited GCs.
We also show in §2 that when the observed densities of in-
dividual clusters are used in our models to predict GCMFs
in different bins of Galactocentric radius (rgc), they fit the
6 Throughout this paper, we use “initial” to mean at a relatively early time
in the development of long-lived clusters, after they have dispersed any rem-
nants of their natal gas clouds, survived the bulk of stellar-evolution mass
loss, and come into virial equilibrium in the tidal field of a galaxy.
much weaker variation of dN/d log M and MTO as functions
of rgc, which is well-known in the Milky Way and other large
galaxies (see Harris 2001; Harris, Harris, & McLaughlin
1998; Barmby, Huchra, & Brodie 2001; Vesperini et al. 2003;
Jordán et al. 2007). Similarly, applying our models to the
GCs in two bins of central concentration, with only the mea-
sured ρh of the clusters in each subsample as input, suffices
to account for previously noted differences between the mass
functions of low- and high-concentration Galactic globulars
(Smith & Burkert 2002). The most fundamental feature of the
GCMF therefore appears to be its dependence on cluster den-
sity, which can be understood at least qualitatively (and even
quantitatively, to within a factor of 2) in terms of evaporation-
dominated cluster disruption.
There is a widespread perception that if the GCMF evolved
slowly from a rising power law at low masses, then a weak
or null variation of MTO with rgc can be achieved only in GC
systems with strongly radially anisotropic velocity distribu-
tions, which are not observed (see especially Vesperini et al.
2003). This apparent inconsistency has been cited to bol-
ster some recent attempts to identify a mechanism by which
a “universal” peak at MTO ∼ 10
5 M⊙ might have been im-
printed on the GCMF at the time of cluster formation, or
very shortly afterwards, and little affected by the subsequent
destruction of lower-mass GCs (e.g., Vesperini & Zepf 2003;
Parmentier & Gilmore 2007). However, given the real suc-
cesses of an evaporation-dominated evolutionary scenario for
the origin of MTO, as summarized above and added to below,
it would be premature to reject the idea in favor of requiring a
near-formation origin, solely on the basis of difficulties with
GC kinematics. (And, in any event, formation-oriented mod-
els must now be reconsidered in light of the non-universality
of MTO as a function of cluster density.)
We are not concerned in this paper with velocity anisotropy
in GC systems, because we only predict an evaporation-
evolved dN/d log M as a function of cluster density (and age)
and take the observed distribution of ρh versus rgc in the Milky
Way as a given, to show consistency with the observed be-
havior of MTO as a function of rgc. Most other models (FZ01;
Vesperini et al. 2003; and references therein) predict dynam-
ically evolved GCMFs directly in terms of rgc, and in doing
so are forced also to derive theoretical dependences of cluster
density on rgc. It is only at this stage that GC orbital distri-
butions enter the problem, and then only in conjunction with
several other assumptions and simplifications. As we discuss
further in §3 below, the radially biased GC velocity distribu-
tions that appear in such models could well be consequences
of one or more of these other assumptions, rather than of the
main hypothesis about evaporation-dominated GCMF evolu-
tion.
2. THE GALACTIC GCMF AS A FUNCTION OF
CLUSTER DENSITY
In this section we define and model the dependence of
the Galactic GCMF on cluster density. First, we describe
the dependence that is expected to arise from evaporation-
dominated evolution.
Two-body relaxation in a tidally limited GC leads to a
roughly steady rate of mass loss, µev ≡ −dM/dt ≃ constant
in time. Thus, the total cluster mass decreases approxi-
mately linearly, as M(t) ≃ M0 −µevt. This behavior is exact
in some classic models of GC evolution (Hénon 1961) and
is found to be a good approximation in most other calcula-
tions (e.g., Lee & Ostriker 1987; Chernoff & Weinberg 1990;
SHAPING THE GCMF BY EVAPORATION 3
Vesperini & Heggie 1997; Gnedin, Lee, & Ostriker 1999;
Baumgardt 2001; Giersz 2001; Baumgardt & Makino 2003;
Trenti, Heggie, & Hut 2007). The result comes from a variety
of computational methods (semi-analytical, Fokker-Planck,
Monte Carlo, and N-body simulation) applied to clusters with
different initial conditions (densities and concentrations) on
different kinds of orbits (circular and eccentric; with and with-
out external gravitational shocks) and with different internal
processes and ingredients (with or without stellar mass spec-
tra, binaries, and central black holes). To be sure, deviations
from perfect linearity in M(t) do occur, but these are generally
small—especially away from the endpoints of the evolution,
i.e., for 0.9&M(t)/M0 & 0.1—and neglecting them to assume
an approximately constant dM/dt is entirely appropriate for
our purposes.
When gravitational shocks are subdominant to relaxation-
driven evaporation, as they generally appear to be for
extant GCs, they work to boost the mass-loss rate µev
slightly without altering the basic linearity of M(t) (e.g.,
Vesperini & Heggie 1997; Gnedin, Lee, & Ostriker 1999; see
also Figure 1 of FZ01). A time-dependent mass scale ∆ ≡
µevt is then associated naturally with any system of coeval
clusters having a common mass-loss rate: all those with initial
M0 ≤ ∆ are disrupted by time t, and replaced with the rem-
nants of objects that began with M0 >∆. As we mentioned
in §1, if the initial GCMF increased towards low masses as a
power law, then ∆ is closely related to a peak in the evolved
distribution, which eventually decreases towards low M <∆
as dN/d log M ∝ M1−β with β = 0 (FZ01).
In standard theory (e.g., Spitzer 1987; Binney & Tremaine
1987, Section 8.3), the lifetime of a cluster against evapora-
tion is a multiple of its two-body relaxation time, trlx. For
a total mass M of stars within a radius r, this scales to first
order (ignoring a weak mass dependence in the Coulomb log-
arithm) as trlx(r) ∝ (Mr
3)1/2 ∝ M/ρ1/2, where ρ∝ M/r3. In a
concentrated cluster with an internal density gradient, trlx(r)
of course varies throughout the cluster, and the global re-
laxation timescale is an average of the local values (see the
early discussion by King 1958). This can still be written as
trlx ∝ M/ρ
1/2, with M the total cluster mass and ρ an appro-
priate reference density. We then have for the instantaneous
mass-loss rate, µev ≡ −dM/dt ∝ M/trlx ∝ ρ1/2. Insofar as this
is approximately constant in time, a GCMF evolving from an
initial β > 1 power law at low masses should therefore de-
velop a peak at a mass that depends on cluster density and age
through the parameter ∆∝ ρ1/2t.
It remains to identify the best measure of ρ in this context.
A standard choice in the literature, and the one that we even-
tually make to derive our main results in this paper, is the
half-mass density ρh = 3M/8πr
h. However, in a steady tidal
field, the mean density ρt inside the tidal radius of a cluster is
constant by definition, and thus choosing ρ = ρt instead is the
simplest way to ensure that µev ∝ ρ
1/2 and µev ≃ constant in
time are mutually consistent. In fact, King (1966) found from
direct calculations of the escape rate at each radius within his
standard (lowered Maxwellian) models, that the coefficient
in µev ∝ ρ
t is only a weak function of the internal density
structure (concentration) of the models, and thus only a weak
function of time for a cluster evolving quasistatically through
a series of such models.
The rule µev ∝ ρ
t is routinely used to set the GC mass-
loss rates in models for the dynamical evolution of the GCMF,
although such studies normally express µev immediately in
terms of orbital pericenters, rp, most often by assuming ρt ∝
r−2p as for GCs in galaxies whose total mass distributions fol-
low a singular isothermal sphere (e.g., Vesperini 1997, 1998,
2000, 2001; Vesperini et al. 2003; Baumgardt 1998; FZ01).
This bypasses any explicit examination of the GCMF as a
function of cluster density, which is our main goal in this pa-
per. But it is done in part because tidal radii are the most
poorly constrained of all structural parameters for GCs in the
Milky Way (their theoretical definition is imprecise and their
empirical estimation is highly model-dependent and sensitive
to low-surface brightness data), and they are exceedingly dif-
ficult if not impossible to measure in distant galaxies. We
deal with this here by focusing on the GCMF as a function of
cluster density ρh inside the less ambiguous, empirically bet-
ter determined and more robust half-mass radius, asking how
simple models with µev ∝ ρ
h fare against the data.
Taking µev ∝ ρ
h in place of µev ∝ ρ
t , which we do to
construct evaporation-evolved model GCMFs in §2.2, is most
appropriate if the ratio ρt/ρh is the same for all clusters and
constant in time. This is the case in Hénon’s (1961) model
of GC evolution, and in this limit (adopted by FZ01 in their
models for the Galactic GCMF) our analysis is rigorously jus-
tified. However, real clusters are not homologous (ρt/ρh dif-
fers among clusters) and they do not evolve self-similarly (ρh
may vary in time even if ρt does not). The key assumption in
our models is that µev is approximately independent of time
for any GC, which is well-founded in any case. By using cur-
rent ρh values to estimate µev, we do not suppose that the half-
mass densities are also constant, but we in effect use a single
number for all GCs to represent a range of (ρt/ρh)
1/2. Equiv-
alently, we ignore a dependence on cluster concentration in
the normalization of µev ∝ ρ
h . As we discuss further in §3,
it is reasonable to neglect this complication in a first approx-
imation because (ρt/ρh)
1/2 varies much less among Galactic
globulars than ρt and ρh do separately. We demonstrate this
explicitly by repeating our analysis with ρh replaced by ρt and
recover essentially the same results for the GCMF.
In §3 we also discuss some recent results, which indicate
that the timescale for relaxation-driven evaporation depends
on a slightly less-than-linear power of trlx (Baumgardt 2001;
Baumgardt & Makino 2003). We point out that this implies
that µev may increase as a modest power of the average sur-
face density of a cluster as well as (or, in an important special
case, instead of) the usual volume density. However, we show
in detail that making the appropriate changes throughout the
rest of the present section to reflect this possibility does not
change any of our conclusions.
2.1. Data
Figure 1 shows the distribution of mass against half-mass
density and against Galactocentric radius for 146 Milky Way
GCs in the catalogue of Harris (1996),7 along with the dis-
tribution of ρh versus rgc linking the two mass plots. The
Harris catalogue actually records the absolute V magni-
tudes of the GCs. We obtain masses from these by apply-
ing the population-synthesis model mass-to-light ratios ΥV
computed by McLaughlin & van der Marel (2005) for indi-
vidual clusters based on their metallicities and an assumed
age of 13 Gyr. However, we first multiplied all of the
7 Feb. 2003 version; see http://physwww.mcmaster.ca/∼harris/mwgc.dat .
http://physwww.mcmaster.ca/~harris/mwgc.dat
4 McLAUGHLIN & FALL
FIG. 1.— Left: Mass versus three-dimensional half-mass density, ρh ≡ 3M/8πr
h , and versus Galactocentric radius, rgc, for 146 Milky Way GCs in the
catalogue of Harris (1996). The dashed line in the first panel is M ∝ ρ
h , a locus of approximately constant lifetime against evaporation. Right: Half-mass
density versus rgc for the same clusters.
McLaughlin & van der Marel ΥV values by a factor of 0.8 so
as to obtain a median Υ̂V ≃ 1.5M⊙L
⊙ in the end,
8 consistent
with direct dynamical estimates (see McLaughlin 2000 and
McLaughlin & van der Marel 2005; also Barmby et al. 2007).
By assigning mass-to-light ratios to GCs in this way, we
allow for expected differences between clusters with differ-
ent metallicities. Our application of a corrective factor to
the population-synthesis values, ΥpopV , is motivated empiri-
cally by the fact that their distribution among Galactic GCs
is strongly peaked around a median Υ̂popV ≃ 1.9 M⊙ L
⊙ , while
the observed (dynamical) ΥdynV lie in a fairly narrow range
around Υ̂dynV ≃ 1.5 M⊙ L
⊙ (McLaughlin & van der Marel
2005). However, it is worth noting that the size of this differ-
ence is similar to what is found in some numerical simulations
of two-body relaxation over a Hubble time in clusters with a
spectrum of stellar masses (e.g., Baumgardt & Makino 2003).
In such simulations, ΥdynV falls below Υ
V due to the prefer-
ential escape of low-mass stars with high individual M∗/L∗
(population-synthesis models do not incorporate this or any
other stellar-dynamical effect). Thus, a median Υ̂dynV < Υ̂
may itself be a signature of cluster evaporation. We might
then also expect that more dynamically evolved clusters—that
is, those with shorter relaxation times—could have systemat-
ically lower ratios of ΥdynV /Υ
V . However, this is a relatively
small effect, which is not well quantified theoretically and
is not clearly evident in real data (the numbers published by
McLaughlin & van der Marel 2005 show no significant corre-
lation between Υ
V and trh for Galactic globulars). We
therefore proceed, as stated, with a single ΥdynV /Υ
V = 0.8
assumed for all GCs.
Harris (1996) gives the projected half-light radius Rh for
141 of the clusters with a mass estimated in this way, and for
these we obtain the three-dimensional half-mass radius from
the general rule rh = (4/3)Rh (Spitzer 1987), which assumes
no internal mass segregation. The remaining five objects have
mass estimates but no size measurements. To each of these
clusters, we assign an rh equal to the median value for those
of the other 141 GCs having masses within a factor two of the
one with unknown rh. In all cases, the half-mass density is
8 Throughout this paper, we use bx to denote the median of any quantity x.
ρh ≡ 3M/8πr
The leftmost panel in Figure 1 shows immediately that
the cluster mass distribution has a strong dependence on
half-mass density: the median M̂ increases with ρh while
the scatter in log M—that is, the width of the GCMF—
decreases. The first of these points is related to the fact that
rh correlates poorly with M (e.g., Djorgovski & Meylan 1994;
McLaughlin 2000). The second point, that the dispersion of
dN/d log M decreases with increasing ρh, is behind the find-
ing (Kavelaars & Hanes 1997; Gnedin 1997) that the GCMF
is broader at very large Galactocentric radii. We return to this
in §2.2.
A natural concern, when plotting M against ρh as we have
done here, is that any apparent correlation might only be a
trivial reflection of the definition ρh ∝ M/r
h. This may seem
particularly worrisome because, as we just mentioned, it is
known that size does not correlate especially well with mass
for GCs in the Milky Way (or, indeed, in other galaxies).
However, the lack of a tight M–rh correlation does not imply
that all GCs have the same rh, even within the unavoidable
measurement errors. The root-mean-square (rms) scatter of
log rh about its average value is ±0.3 for Galactic GCs, and
the 68-percentile spread in log rh is slightly greater than 0.5,
or more than a factor of 3 in linear terms (from the data in
Harris 1996; see, e.g., Figure 8 of McLaughlin 2000). This
compares to an rms random measurement error (from formal,
χ2 fitting uncertainties) of δ(log rh) ≈ 0.05, or about 10% rel-
ative error; and an rms systematic measurement error (i.e.,
differences in the rh inferred from fitting different structural
models to a single cluster) of perhaps δ(log rh) . 0.03; see
McLaughlin & van der Marel (2005). Most of the scatter in
plots of observed half-light radius versus mass is therefore
real and contains physical information. The left-hand panel
of Figure 1 displays this information in a form that highlights
clear, nontrivial overall trends requiring physical explanation.
The dashed line in the plot of mass against density traces
the proportionality M ∝ ρ
h , or Mr
h = constant. Insofar as
the half-mass relaxation time scales as trh ∝ (Mr
1/2, and
to the extent that µev ∝ M/trh ∝ ρ
h approximates the av-
erage rate of relaxation-driven mass loss, this line is one of
equal evaporation time. That such a locus nicely bounds
the lower envelope of the observed cluster distribution is
SHAPING THE GCMF BY EVAPORATION 5
TABLE 1
MILKY WAY GC PROPERTIES IN BINS OF DENSITY AND GALACTOCENTRIC RADIUS
Bin N bρh
brgc a Mmin Mmax bM
a MTO
[M⊙ pc−3] [kpc] [M⊙] [M⊙] [M⊙] [M⊙]
ρh bins
0.034 ≤ ρh ≤ 76.5 M⊙ pc
−3 48 8.48 12.9 5.63× 102 8.84× 105 4.12× 104 3.98× 104
78.8 ≤ ρh ≤ 526 M⊙ pc
−3 49 232 5.6 8.37× 103 1.67× 106 1.22× 105 1.58× 105
579 ≤ ρh ≤ 5.65× 10
4 M⊙ pc−3 49 973 3.2 1.93× 104 1.30× 106 2.82× 105 2.88× 105
rgc bins
0.6 ≤ rgc ≤ 3.2 kpc 47 597 1.9 4.47× 103 1.02× 106 1.15× 105 2.14× 105
3.3 ≤ rgc ≤ 9.4 kpc 50 261 5.2 2.02× 103 1.67× 106 1.27× 105 1.66× 105
9.6 ≤ rgc ≤ 123 kpc 49 18.4 18.3 5.63× 102 1.30× 106 7.42× 104 8.71× 104
a The notation bx represents the median of quantity x.
b MTO is the peak mass of the model GCMFs traced by the solid curves in each panel of Figure 2, which
are given by equation (3) of the text with β = 2, Mc = 10
6 M⊙, and individual ∆ given by the observed ρh
of each cluster through equation (4).
itself a strong hint that relaxation-driven cluster disruption
has significantly modified the GCMF at low masses (re-
call that Mr3h = constant defines one side of the GC “sur-
vival triangle” when the M–ρh plot is recast as rh versus M:
Fall & Rees 1977; Okazaki & Tosa 1995; Ostriker & Gnedin
1997; Gnedin & Ostriker 1997). It is also further evidence
that the weak correlation of observed rh with M is due to sig-
nificant and real differences in cluster radii, since if rh were
intrinsically the same for all GCs, then we would see M ∝ ρh
instead.
The middle panel of Figure 1 shows the well-known result
that the typical GC mass depends weakly if at all on Galacto-
centric radius, at least until large rgc & 30–40 kpc, where there
are too few clusters to discern any trend. The right-hand panel
of the figure shows why this is true even though the GCMF
depends significantly on cluster density: although there is a
correlation between half-mass density and Galactocentric po-
sition, the large scatter about it is such that convolving the
observed M versus ρh with the observed ρh versus rgc results
in an almost null dependence of M on rgc.
We now divide the GC sample in Figure 1 roughly into
thirds, in two different ways: first on the basis of half-mass
density, and second by Galactocentric radius. These ρh and
rgc bins are defined in Table 1, which also gives a few sum-
mary statistics for the globulars in each subsample. We count
the clusters in every subsample in about 10 equal-width bins
of log M to obtain histogram representations of dN/d log M,
first as a function of ρh and then as a function of rgc. These
GCMFs are shown by the points in Figure 2, with errorbars
indicating standard Poisson uncertainties. The curves in the
figure trace model GCMFs, which we describe in §2.2. For
the moment, it is important to note that the dashed curve is the
same in every panel, apart from minor differences in normal-
ization, and is proportional to the GCMF for the whole sample
of 146 GCs. (In the middle-left panel of Figure 2, which per-
tains to clusters distributed tightly around the median ρh of
the entire GC system, the dashed curve is coincident with the
solid curve running through the data.)
The left-hand panels of Figure 2 show directly that the
GCMF is peaked for clusters at any density, and that the mass
of the peak increases systematically with ρh (see also the last
column of Table 1, but note that the turnover masses there
refer to the model GCMFs that we develop below). The sta-
tistical significance of this is very high, and qualitatively it
is the behavior expected if MTO owes its existence to cluster
disruption at a rate that increases with ρh, as is the case with
relaxation-driven evaporation.
The right-hand panels of Figure 2 confirm once again that
the GCMF peak mass is a very weak function of Galactocen-
tric position. In fact, the observed distributions in the two rgc
bins inside ≃10 kpc are statistically indistinguishable in their
entirety, and the main difference at larger rgc & 10 kpc is a
slightly higher proportion of low-mass clusters rather than a
large change in MTO. All of this is consistent with the pri-
mary dependence of the GCMF being that on ρh, since Fig-
ure 1 shows that the GC density distribution is not sensitive to
Galactocentric position for rgc . 10–20 kpc but has a substan-
tial low-density tail at larger radii (with a broader associated
GCMF, as seen in the upper-left panel of Figure 2).
2.2. Simple Models
We now assess more quantitatively whether these results
are consistent with evaporation-dominated evolution of the
GCMF from an initial distribution like that observed for
young clusters in the local universe. We model the time-
evolution of the distribution of M versus ρh in Figure 1 but do
not attempt this for the distribution of ρh over rgc—the details
of which likely depend on a complicated interplay between
the tidal field of the Galaxy, the present and past orbital pa-
rameters of clusters, and the structural nonhomology of GCs.
To compare our models to the current GCMF as a function of
rgc, we simply calculate them using the observed ρh of indi-
vidual clusters in different ranges of Galactocentric radius.
We assume that the initial GCMF was independent of clus-
ter density, and that all globulars surviving to the present day
have been losing mass for the past Hubble time at constant
rates. We use the current half-mass density of each cluster
to estimate µev ∝ ρ
h . As we discussed earlier, an approxi-
mately time-independentµev is indicated by most calculations
of two-body relaxation in tidally limited GCs. We give a more
detailed, a posteriori justification in §3 for using ρh, rather
than other plausible measures of cluster density, to estimate
Consider first a group of coeval GCs with an initial mass
function dN/d log M0 and a single, time-independent mass-
loss rate µev. The mass of every cluster decreases linearly as
M(t) = M0 −µevt, and at any later time each has lost the same
amount ∆ ≡ M0 − M(t) = µevt. FZ01 show rigorously that in
6 McLAUGHLIN & FALL
FIG. 2.— GCMF as a function of half-mass density, ρh ≡ 3M/8πr
h (left panels), and as a function of Galactocentric radius, rgc (right panels), for 146 Milky
Way GCs in the catalogue of Harris (1996). The dashed curve in all cases is an evolved Schechter function for the entire GC system (Jordán et al. 2007): equation
(3) with β = 2, Mc = 106 M⊙, and ∆≡ 2.3×105 M⊙ for all clusters (from equation [4] and a median bρh = 246 M⊙ pc
−3), giving a peak at MTO = 1.6×10
5 M⊙ .
Solid curves are the GCMFs predicted by equation (3) with β = 2 and Mc = 106 M⊙ but individual ∆ given by the observed ρh of each cluster (equation [4]) in
the different subsamples.
this case, the evolved and initial GCMFs are related by
d log M
d log M0
(M +∆)
d log (M +∆)
. (1)
This is the basis for the claim that the mass function scales
generically as dN/d log M ∝ M+1 (a β = 0 power law) at low
enough M(t)<∆—that is, for the surviving remnants of clus-
ters with M0 ≈∆—just so long as the initial distribution was
not a delta function.
We follow FZ01 (see also Jordán et al. 2007) in adopting a
Schechter (1976) function for the initial GCMF:
dN/d log M0 ∝ M
0 exp
−M0/Mc
. (2)
With β ≃ 2, this distribution describes the power-law mass
functions of young massive clusters in systems like the Anten-
nae galaxies (e.g., Zhang & Fall 1999). An exponential cut-
off at Mc ∼ 10
6 M⊙ is generally consistent with such data,
even if not always demanded by them; here we require it
mainly to match the curvature observed at high masses in old
GCMFs (e.g., Burkert & Smith 2000; Jordán et al. 2007).
Combining equations (1) and (2) gives the probability den-
sity that a single GC with known evaporation rate and age has
an instantaneous mass M. The time-dependent GCMF of a
system of N GCs with a range of µev (or ages, or both) is
then just the sum of all such individual probability densities:
d log M
[M +∆i]
M +∆i
. (3)
Here the total mass losses ∆i = (µevt)i may differ from clus-
ter to cluster (ti being the age of a single GC) but both β and
Mc are assumed to be constants, independent of ρh in particu-
lar.9 Given each ∆i, the normalizations Ai in equation (3) are
defined so that the integral over d log M of each term in the
summation is unity.
Jordán et al. (2007) have introduced a specialization of
equation (3) in which all clusters have the same ∆. They refer
to this as an evolved Schechter function and describe its prop-
erties in detail (including giving a formula for the turnover
mass MTO as a function of ∆ and Mc) for the case β = 2. Here
we note only that, at very young cluster ages or for slow mass-
loss rates, such that ∆ ≪ Mc and only the low-mass, power-
law part of the initial GCMF is significantly eroded, any one
evolved Schechter function has a peak at MTO ≃ ∆/(β − 1)
(for β > 1). As ∆ increases relative to Mc, the turnover at
first increases proportionately and the width of the distribu-
tion decreases (since the high-mass end at M & MTO is largely
unchanged). For large ∆≫Mc, however, the peak is bounded
above by MTO →Mc and the width approaches a lower limit.
Thus, the dependence of MTO on ∆ is weaker than linear when
Mc is finite in the initial GCMF of equation (2). Any peak in
the full equation (3) for a system of GCs with individual ∆
values is an average of N different turnovers and must be cal-
culated numerically.
In their modeling of the Milky Way GC system, FZ01 ef-
9 Note that Mc appears to take on different values in the GCMFs of other
galaxies, varying systematically with the total luminosity Lgal (Jordán et al.
2007). The reasons for this are unclear, as is the origin of this mass scale in
the first place.
10 The increase of MTO and the decrease of the full width of dN/d log M
for increasing ∆ eventually saturate when the mass loss per GC is so high
that it affects clusters in the exponential part of the initial Schechter-function
GCMF. This is because dN/d log M ∝ M+1 exp(−M/Mc) is a self-similar
solution to equation (1).
SHAPING THE GCMF BY EVAPORATION 7
fectively compute mass functions of the type (3)—based on
the same initial conditions and dynamical evolution—with a
distribution of ∆ values determined by the orbital parame-
ters of clusters in an idealized, spherical and static logarith-
mic Galaxy potential (used both to fix µev in terms of clus-
ter tidal densities and to estimate additional mass loss due
to gravitational shocks). Jordán et al. (2007) fit GCMF data
in the Milky Way and scores of Virgo Cluster galaxies with
their version of equation (3) in which all GCs have the same
∆. They thus estimate the dynamical mass loss from typi-
cal clusters in these systems. Here, we construct models for
the Milky Way GCMF using ∆ values given directly by the
observed half-mass densities of individual GCs.
We adopt β = 2 for the initial low-mass power-law index
in equation (2), which carries over into equation (3) for the
evolved dN/d log M. Jordán et al. (2007) have fitted the full
Galactic GCMF with an evolved Schechter function assuming
β = 2 and a single ∆ ≡ ∆̂ for all surviving globulars. They
find Mc ≃ 10
6 M⊙ and ∆̂ = 2.3× 10
5 M⊙. We use this value
of Mc in equation (3) and we associate ∆̂ with the mass loss
from clusters at the median half-mass density of the entire GC
system, which is ρ̂h = 246 M⊙ pc
−3 from the data in Figure 1.
Since we are assuming that ∆ = µevt ∝ ρ
h t for coeval GCs,
we therefore stipulate
∆ = 1.45× 104 M⊙
ρh/M⊙ pc
−3)1/2 (4)
for globulars with arbitrary ρh. Assuming a typical GC age of
t = 13 Gyr, this corresponds to a mass-loss rate of
µev ≃ 1100 M⊙ Gyr
−1 (ρh/M⊙ pc−3
. (5)
In §3 we discuss the cluster lifetimes implied by this value
of µev. We emphasize here that the scaling of µev and ∆
with ρ
h follows rather generically from our hypothesis of
evaporation-dominated cluster evolution, while the numerical
coefficients in equations (4) and (5) are specific to the assump-
tion of β = 2 for the power-law index at low masses in the
initial GCMF.
The dashed curve shown in every panel of Figure 2 is
the evolved Schechter function fitted to the entire GCMF of
the Milky Way by Jordán et al. (2007). This has a peak at
MTO ≃ 1.6× 10
5 M⊙ (magnitude MV ≃ −7.4 for a typical V -
band mass-to-light ratio of 1.5 in solar units) and gives a very
good description of the observed dN/d log M in the middle
density bin, 79 . ρh . 530 M⊙ pc
−3, and in the two inner ra-
dius bins, rgc ≤ 9.4 kpc. This is expected, since the median
half-mass density in each of these cluster subsamples is very
close to the system-wide median ρ̂h = 246 M⊙ pc
−3 (see Ta-
ble 1). Even in the outermost rgc bin, a Kolmogorov-Smirnov
(KS) test only marginally rejects the dashed-line model (at
the ≃95% level), because this subsample still includes many
GCs at or near the global median ρ̂h (see Figure 1). By con-
trast, the average GCMF is strongly rejected as a model for
the lowest- and highest-density GCs on the left-hand side of
Figure 2: the KS probabilities that these data are drawn from
the dashed distribution are <10−4 in both cases. This is also
expected since, by construction, these bins only contain clus-
ters with densities well away from the median of the full GC
system, for which the total mass lost by evaporation should be
significantly different from the typical ∆̂ = ∆(ρ̂h).
The solid curves in Figure 2, which are different in ev-
ery panel, are the superpositions of many different evolved
Schechter functions, as in equation (3), with distinct ∆ values
given by equation (4) using the observed ρh of each cluster in
the corresponding subsample. These models provide excel-
lent matches to the observed dN/d log M in every ρh and rgc
bin, with χ2 < 1.3 per degree of freedom in all cases. This is
the main result of this paper.
The last column of Table 1 gives the mass MTO at which
each of the solid model GCMFs in Figure 2 peaks. We note
that these turnovers increase roughly as MTO ∼ ρ̂
0.3−0.4
h for
our specific binnings in ρh and rgc, somewhat shallower than
the ρ
h scaling of the cluster mass-loss rate that defines the
models. This is partly because of the averaging over indi-
vidual turnovers implied by the summation of many evolved
Schechter functions in each GC bin, and partly because—as
we discussed just after equation (3)—the turnover mass of any
one evolved Schechter function cannot increase indefinitely in
direct proportion to ∆ ∝ ρ
h , but has a strict upper limit of
MTO ≤ Mc.
Our models are naturally consistent with the fact that the
GCMF is narrower for clusters with higher densities. This
is obvious in the left-hand panels of Figure 2; in the dis-
cussion immediately after equation (3), we described how it
follows from the increase of MTO with ∆ ∝ ρ
h for a sin-
gle evolved Schechter function. In addition, the superposi-
tion of many such functions with separate, density-dependent
turnovers and widths results in wider GCMFs for cluster sub-
samples spanning larger ranges of ρh. This accounts in partic-
ular for the breadth of the mass function at rgc ≥ 9.4 kpc. The
globulars at these radii have 0.034 ≤ ρh ≤ 4.1×10
3 M⊙ pc
corresponding to individual evolved Schechter functions with
turnovers at 2.7× 103 . MTO . 4.0× 10
5 M⊙. The compos-
ite GCMF in the lower-right panel of Figure 2 is therefore
extremely broad and shows a very flat peak, such that an over-
all MTO cannot be established precisely from the data alone.
This explains the findings of Kavelaars & Hanes (1997), who
pointed out that the GCMF of the outermost third of the Milky
Way cluster system has a turnover that is statistically consis-
tent with the full-Galaxy average, but a larger dispersion (see
also Gnedin 1997).
Finally, if the GCMF evolved dynamically from initial con-
ditions similar to those we have adopted, then the data and
models in the left-hand panels of Figure 2 argue against
the notion that external gravitational shocks, rather than in-
ternal two-body relaxation, were primarily responsible for
shaping the present-day GCMF. This is because the mass-
loss rate caused by shocks alone, −dM/dt = µsh ∝ M/ρh,
differs significantly from that caused by evaporation alone,
−dM/dt = µev ∝ ρ
h . The direct dependence of µsh on M en-
sures that shocks become progressively less important com-
pared to evaporation as clusters lose mass (at a given ρh), and
consequently shocks are not likely to have had much effect on
the observed GCMF for M < MTO. Furthermore, the inverse
dependence of µsh on ρh is contrary to the direct dependence
of MTO on ρh shown in Figure 2. The different roles played
by shocks and evaporation in shaping the observed GCMF are
discussed more fully by FZ01. We note here that gravitational
shocks may have been important in destroying very massive
or very low-density clusters early in the history of our Galaxy.
2.3. Other Cluster Properties
If the current shape of the GCMF is fundamentally the re-
sult of long-term cluster disruption according to a mass-loss
8 McLAUGHLIN & FALL
rule like µev ∝ ρ
h , then it should be possible to reproduce
the distribution as a function of any other cluster attribute by
using the observed ρh of individual GCs in equations (3) and
(4) to build model dN/d log M for subsamples of the Galactic
cluster system defined by that attribute—as we did for the rgc
binning of §2.2. Here we explore one example in which dif-
ferences in the GCMFs of two groups of globulars can be seen
in this way to follow from differences in their ρh distributions.
Smith & Burkert (2002) have shown that the mass function
of Galactic globulars with King (1966) model concentrations
c < 0.99 has a less massive peak than that for c ≥ 0.99. [Here
c ≡ log(rt/r0), where rt is the fitted tidal radius and r0 a core
scale.] They further find that a power-law fit to the low-
c GCMF just below its peak returns dN/d log M ∝ M+0.5—
shallower than the M+1 expected generically for a mass-loss
rate that is constant in time—but they confirm that the latter
slope applies for the GCMF at c ≥ 0.99. They discuss various
options to explain these results, including a suggestion that, if
the mass functions of both low- and high-concentration clus-
ters evolved slowly from the same, young-cluster–like initial
distribution, then the mass-loss law for low-c GCs may have
differed from that for high-c clusters. However, they give no
physical explanation for such a difference, and we can show
now that none is required.
The upper panel of Figure 3 plots concentration against
half-mass density for the same 146 GCs from Figure 1; the
filled circles distinguish 24 clusters with c < 0.99. There
is a correlation of sorts between c and ρh, which either de-
rives from or causes the better-known correlation between c
and M (e.g., Djorgovski & Meylan 1994; McLaughlin 2000).
The important point here is that the ρh distribution is off-
set to lower values and has a higher dispersion at c < 0.99.
Following the discussion in §2.2, we therefore expect the
low-concentration GCMF to have a smaller MTO, a flatter
shape around the peak, and a larger full width than the high-
concentration GCMF.
The lower panel of Figure 3 shows the GCMFs for c < 0.99
(filled circles) and c ≥ 0.99 (open circles). The curves are
again given by equation (3) with β = 2, Mc = 10
6 M⊙, and in-
dividual ∆ calculated from the observed cluster ρh through
equation (4). These models peak at MTO ≃ 4.3× 10
for the c < 0.99 subsample but at MTO ≃ 1.8× 10
5 M⊙ for
c ≥ 0.99, entirely as a result of the different ρh involved.
The larger width of dN/d log M and its shallower slope at
any M . 105 M⊙ for the low-concentration GCs are also
clear, in the model curves as well as the data. It is further
evident that there are no low-c Galactic globulars observed
with M & 2× 105 M⊙, above the nominal turnover of the full
GCMF (as Smith & Burkert 2002 noted). But this is not sur-
prising, given that there are so few low-concentration clusters
in total and they are expected to be dominated by low-mass
objects because of their generally low densities. Thus, the
solid curve in Figure 3 predicts perhaps ≃ 3 high-mass clus-
ters with c < 0.99, where none is found.
The apparent variation of the Milky Way GCMF with in-
ternal concentration is therefore consistent with the same
density-based model for evaporation-dominated dynamical
evolution that we compared to dN/d log M as a function of
ρh and rgc in §2.2. To show this, we have made use of the
densities ρh exactly as observed within the two concentration
bins indicated in Figure 3—just as we also took ρh directly
from the data for GCs in different ranges of rgc to construct
models for comparison with the observed dN/d log M in the
FIG. 3.— Top: Concentration parameter as a function of half-mass density
for 146 Galactic GCs. The line of points at c ≡ 2.5 comes from the practice of
assigning this value to core-collapsed clusters in the Harris (1996) catalogue
and its sources. Bottom: GCMF data and models (eqs. [3] and [4]) for 24
clusters with c < 0.99 (filled circles and solid curve) and 122 clusters with
c ≥ 0.99 (open circles and dashed curve).
right-hand panels of Figure 2. Of course, this is not the same
as explaining the distribution of ρh versus rgc or c. Doing so
would certainly be of interest in its own right, but it is beyond
the scope of our work here.
3. DISCUSSION
In this section, we first show that the mass-loss rate in equa-
tion (5) above implies cluster lifetimes that compare favorably
with those expected from relaxation-driven evaporation. Then
we discuss why it is reasonable to approximate µev ∝ ρ
the first place. Finally, we address the issue of possible con-
flict, in some other models for evaporation-dominated GCMF
evolution, between the near-constancy of MTO as a function
of rgc and the observed kinematics of GC systems.
3.1. Cluster Lifetimes
The disruption time of a GC with mass M and a steady
mass-loss rate µev is just tdis = M/µev. It is convenient,
for purposes of comparison with evaporation times in the
literature, to normalize tdis to the relaxation time of a
cluster at its half-mass radius. In general, this is trh =
0.138M1/2r
G1/2m∗ ln (γM/m∗)
, where m∗ is the mean
SHAPING THE GCMF BY EVAPORATION 9
stellar mass. For clusters of stars with a single mass,
m∗ ≃ 0.7M⊙ and γ ≃ 0.4 are appropriate (Spitzer 1987;
Binney & Tremaine 1987, equation [8-72]), in which case
equation (5) for µev from our GCMF modeling implies
µevtrh
0.57M/M⊙
0.57× 105
. (6)
Clusters with realistic stellar mass spectra will have slightly
different values of m∗ and a smaller γ in the calculation of
the relaxation time (Giersz & Heggie 1996), which changes
the numerical value of tdis/trh somewhat but does not alter any
scalings.
We obtained the normalization of µev ∝ ρ
h in §2.2 by
fitting to observed GCMFs constructed by applying a spe-
cific mass-to-light ratio ΥV to every cluster, with models as-
suming a specific form for the initial dN/d log M0. Thus,
the result in equation (6) depends both on the median Υ̂V
and on the power-law index β at low masses in the original
Schechter-function GCMF. The net scaling, for either single-
or multiple-mass clusters, is
tdis/trh ∝ Υ̂
V (β − 1)
−1 . (7)
To see the dependence of this dimensionless lifetime on
Υ̂V , note that we require µev ∝∆∝ΥV to fit the mass losses
of clusters with a given distribution of luminosities (the di-
rect observables), whereas M/trh is proportional to ρ
V (L/r
1/2. Therefore, tdis/trh ∝ (M/trh)/µev ∝Υ
V . The
mass-to-light ratios adopted in this paper, with a median value
Υ̂V ≃ 1.5 M⊙ L
⊙ , are tied directly to dynamical determina-
tions (§2.1).
To understand the dependence on β in equation (7), recall
first that the coefficients in our expressions for ∆ and µev
as functions of ρh (eqs. [4] and [5]) followed from choos-
ing β = 2 for the power-law exponent at low masses in the
initial GCMF (equation [2]). As we mentioned just after
equation (3), the turnover mass of an evolved Schechter func-
tion with any β > 1 is MTO ≃ ∆/(β − 1) in the limit of low
∆ ∝ ρ
h , and MTO → Mc for very high ∆. In this sense,
the strongest observational constraints on the normalizations
of ∆ and µev come from the low-density clusters. All other
things being equal, their GCMF can be reproduced with β 6= 2
if ∆ and µev are multiplied by (β − 1) at fixed ρh. Therefore,
tdis ∝ 1/µev ∝ 1/(β− 1). Observations of young massive clus-
ters (e.g., Zhang & Fall 1999) indicate that β is near 2; but if
it were slightly shallower, then the cluster lifetimes we infer
from the old GCMF would increase accordingly. Even a rela-
tively minor change to β = 1.5 would double tdis/trh from ≈10
to ≈20.
In the model of Hénon (1961) for single-mass clusters
evolving self-similarly (fixed ratio ρt/ρh of mean densities in-
side the tidal and half-mass radii) in a steady tidal field (ρt
constant in time), a cluster loses 4.5% of its remaining mass
every half-mass relaxation time. The time to complete disrup-
tion is therefore tdis/trh = 1/0.045≃ 22. For non-homologous
clusters in a steady tidal field, tdis/trh is a function of central
concentration and can differ from the Hénon value by factors
of about two. From one-dimensional Fokker-Planck calcula-
tions, Gnedin & Ostriker (1997) find tdis/trh ≃ 10–40 for King
(1966) model clusters with c values similar to those found in
real GCs and with gravitational shocks suppressed (see their
Figure 6). Thus, even though the evaporation time in equation
(6) may be slightly shorter than is typically found in theoreti-
cal calculations, it is certainly within the range of such calcu-
lations. Moreover, the assumptions of a steady tidal field and
a single stellar mass in Hénon (1961) and Gnedin & Ostriker
(1997) are important. Part of the difference between the typ-
ical lifetimes in these particular theoretical treatments and
our estimate of tdis/trh from the GCMF is that the former do
not include gravitational shocks, which may have accelerated
somewhat the evolution of real clusters (although we stress
again that shocks do not appear in general to have dominated
the evolution of extant Galactic GCs and are not expected
to affect the basic time-independence of the net mass-loss
rate; see Vesperini & Heggie 1997, Gnedin, Lee, & Ostriker
1999, FZ01, and Prieto & Gnedin 2006). A spectrum of stel-
lar masses in the clusters may also have contributed to an in-
crease in evaporation rate over the single-mass values (e.g.,
Johnstone 1993; Lee & Goodman 1995).
Estimates of evaporation times from other numerical
methods and for models of multimass clusters can be rather
sensitive to the detailed computational techniques and input
assumptions and approximations, and differences at roughly
the factor-of-two level in tdis/trh between different analyses
are not uncommon; see, e.g., Vesperini & Heggie (1997),
Takahashi & Portegies Zwart (1998, 2000), Baumgardt
(2001), Joshi, Nave, & Rasio (2001), Giersz (2001), and
Baumgardt & Makino (2003). Thus, although the lifetimes
in these studies tend to be broadly comparable to those in
Hénon (1961) and Gnedin & Ostriker (1997), noticeably
shorter values do occur in some models. In any case, we
are encouraged by consistency to within factors of two or
three between estimates of tdis or µev by such vastly different
methods—one purely observational, based on the mass
functions of cluster systems; the other purely theoretical,
based on idealized models for the evolution of individual
clusters—particularly since each method involves several
uncertain inputs and parameters.
3.2. Approximating µev ∝ ρ
3.2.1. Half-mass versus Tidal Density
The dimensionless disruption time in equation (6) is inde-
pendent of any cluster property other than the Coulomb log-
arithm because we have used GC half-mass densities to esti-
mate tdis = M/µev ∝ M/ρ
h , while trh also scales as M/ρ
However, as we mentioned above, the Fokker-Planck calcu-
lations of Gnedin & Ostriker (1997) in particular show that
tdis/trh is actually a function of central concentration, c, for
King (1966) model clusters in steady tidal fields. The constant
of proportionality in µev ∝ ρ
h should therefore also depend
on c, a detail that we have neglected to this point. We show
now that this has not biased any of our analysis or affected our
conclusions.
The dotted curve in Figure 4 illustrates the dependence of
tdis/trh on c for single-mass King models, as given by equation
(30) of Gnedin & Ostriker (1997). The solid curve is propor-
tional to (ρh/ρt)
1/2 = (r3t /2r
1/2, which we have calculated
as a function of c for these models and multiplied by a con-
stant to compare directly with tdis/trh. Evidently, there is an
approximate equality tdis/trh ≈ 2.15(ρh/ρt)
1/2, which holds to
within <15% over the range of concentrations shown in Fig-
ure 4 (note that all but 6 Galactic GCs have 0.7 ≤ c ≤ 2.5,
corresponding to central potentials 3 .W0 . 11). Thus, if the
10 McLAUGHLIN & FALL
FIG. 4.— Dependence of tdis/trh (dotted line; from Gnedin & Ostriker
1997) and (ρh/ρt )
1/2 (solid line; after scaling by a factor of 2.15) on cen-
tral concentration for single-mass King-model clusters. Over the range of c
shown, which includes nearly all Galactic globulars, the approximate propor-
tionality tdis/trh ∝ (ρh/ρt )
1/2 holds to within better than 15%. Thus, to this
level of accuracy the evaporation time tdis is roughly the same multiple of
t for clusters with any internal density profile.
evaporation time is written as tdis ∝ trh(ρh/ρt)
1/2 ∝ M/ρ
then the constant of proportionality in the mass-loss rate µev ∝
M/tdis ∝ ρ
t should be nearly independent of c. In fact, King
(1966) originally concluded, from quite basic arguments, that
the evaporation rate of a cluster with a lowered-Maxwellian
velocity distribution would take the form µev ∝ ρ
t with
only a weak dependence on c. An essentially concentration-
independent scaling of µev with ρ
t is also found in N-
body simulations of tidally limited, multimass clusters (e.g.,
Vesperini & Heggie 1997) and so is not an artifact of any as-
sumptions specific to the calculations of either King (1966) or
Gnedin & Ostriker (1997).
This suggests that it might have been more natural to spec-
ify cluster evaporation rates proportional to ρ
t rather than
h when developing our GCMF models in §2. For any clus-
ter in a steady tidal field, with a constant ρt , such a choice
would also have been automatically consistent with an ap-
proximately time-independent µev and the corresponding lin-
ear M(t) dependence that we have adopted throughout this
paper. As we discussed at the beginning of §2, our deci-
sion to work with ρh rather than ρt was motivated by the
fact that the half-mass density is much better defined in prin-
ciple and more accurately observed in practice. Neverthe-
less, re-writing µev ∝ ρ
t as µev ∝ (ρt/ρh)
1/2 × ρ
h makes
it clear that the validity of our models, with a fixed coefficient
in µev ∝ ρ
h , depends on the extent to which variations in
(ρt/ρh)
1/2 can safely be ignored.
Figure 4 shows that the full range of possible values for
(ρh/ρt)
1/2 in King-model clusters with c ≥ 0.7 is only a factor
of ≃ 4 between minimum and maximum. Therefore, using
a single, intermediate value of this density ratio to describe
all GCs (or a single GC evolving in time through a series of
quasi-static King models)—which we have effectively done
by using a GCMF fit to normalize ∆ and µev in equations (4)
and (5)—should never be in error by more than a factor of 2 or
so. This is a relatively small inaccuracy, given that measured
GC densities range over four to five orders of magnitude.
To confirm more directly that our models with µev ∝ ρ
are good approximations to GCMF evolution under a mass-
loss law µev ∝ ρ
t , we have repeated the analysis of §2 in full
but using the GC tidal densities ρt (derived from the values of
rt listed by Harris 1996) in place of ρh throughout. All of our
main results persist.
For example, the two panels of Figure 5, which are analo-
gous to the left- and rightmost panels of Figure 1 above, show
that (1) the GC mass distribution has a clear dependence on
ρt , with a lower envelope that is well matched by a line of
constant evaporation time, M ∝ ρ
t (the dashed line in the
plot); and (2) although the scatter in the distribution of ρt
over Galactocentric radius is smaller than the scatter in ρh ver-
sus rgc, it is still significant. Because the M–rgc distribution
can now be viewed as the convolution of the M–ρt distribu-
tion with the ρt–rgc distribution, the scatter in ρt versus rgc is
again critical in explaining the weak or null dependence of the
GCMF on Galactocentric radius. (The M–rgc distribution is,
of course, unchanged from that shown in the middle panel of
Figure 1.)11
Figure 6 shows the Milky Way GCMF for globulars in three
equally populated bins of tidal density (defined as indicated in
the left-hand panels of the plot) and in the same three bins of
Galactocentric radius that we used in §2.2 above. Our mod-
els for these distributions are based as before on equation (3)
with β = 2, but now the total mass lost from any GC is esti-
mated from its tidal density rather than its half-mass density.
Specifically, we take
∆ = 2.1× 105 M⊙
ρt/M⊙ pc
−3)1/2 . (8)
The numerical coefficient in equation (8) is such that it gives a
∆ identical to that in equation (4) for a GC with ρh/ρt = 210,
which is the median value of this density ratio for the 146 GCs
in the Harris (1996) catalogue.
As in Figure 2, the dashed curve in every panel of Figure
6 is the same, representing a fit to the average dN/d log M of
the entire Galactic GC system. Thus, it is immediately clear
that the peak mass of the GCMF increases significantly and
systematically with increasing ρt , just as it does with increas-
ing ρh. Meanwhile, the solid curves are subsample-specific
model GCMFs, obtained by using the observed tidal density
of each cluster in any ρt or rgc bin to specify individual ∆ val-
ues via equation (8) for each of the evolved Schechter func-
tions in the summation of equation (3). As expected, there is
no appreciable difference, in terms of the fits to any of the ob-
served GCMFs, between these models based on evaporation
rates µev ∝ ρ
t and our original models with µev ∝ ρ
3.2.2. Retarded Evaporation
Another potential concern comes from recent arguments
(see especially Baumgardt 2001; Baumgardt & Makino 2003)
11 As was also the case with our earlier plots involving ρh in Figure 1,
the scatter and structure in both panels of Figure 5 are real, since the rms
scatter of log rt about the best-fit lines to either of log M or log rgc is 0.3–0.35
while the rms errorbars based on formal fitting uncertainties are in the range
δ(log rt ) ≃ 0.05–0.15 for a variety of models (McLaughlin & van der Marel
2005).
SHAPING THE GCMF BY EVAPORATION 11
FIG. 5.— Scatter plots of mass M versus mean density inside the tidal radius (ρt ≡ 3M/4πr3t ) and of ρt versus Galactocentric radius rgc, for 146 Galactic GCs
from the Harris (1996) catalogue. These plots are analogous to the left- and rightmost panels of Figure 1. The dashed line in the left-hand plot traces the relation
M ∝ ρ
t , which defines a locus of constant evaporation time for µev ∝ ρ
FIG. 6.— Observed GCMF (points, with Poisson errorbars) and models (curves) as a function of mean cluster density inside the tidal radius, ρt ≡ 3M/4πr3t (left-
hand panels), and as a function of Galactocentric radius, rgc (right-hand panels). The dashed curve in every panel is an evolved Schechter function representing
the entire GC system: equation (3) with β = 2, Mc = 106 M⊙, and a single ∆, common to all clusters, evaluated from equation (8) using the median bρt of all 146
Galactic GCs. Solid curves are subsample-specific models using equation (3) with β = 2 and Mc = 106 M⊙ but a different ∆ value for every cluster (obtained
from equation [8] using individual observational estimates of ρt ) in any ρt or rgc bin.
12 McLAUGHLIN & FALL
that the total evaporation time of a tidally limited cluster
is not simply a multiple of an internal two-body relaxation
time, trlx ∝ (Mr
3)1/2, but depends on both trlx and the crossing
time tcr ∝ (M/r
3)−1/2 through the combination tdis ∝ t
with x < 1. The mass-loss rate µev ∝ M/tdis then scales as
M3/2−xr−3/2, which for x 6= 1 differs from the rates µev ∝ ρ
and µev ∝ ρ
t that we have so far adopted. However, our
GCMF models are still meaningful, because postulating tdis ∝
txrlxt
cr implies a dependence of µev on a measure of cluster
density that is, once again, well approximated by ρ
h for
Galactic GCs. Before showing this, we briefly discuss the
reasons and the evidence for a possible dependence of tdis on
both trlx and tcr.
If stars are assumed to escape a cluster as soon as they
have attained energies above some critical value as a result
of two-body relaxation, then tdis ∝ trlx is expected (and con-
firmed by N-body simulations; e.g., Baumgardt 2001). How-
ever, more complicated behavior may arise when escape not
only depends on stars satisfying such an energy criterion, but
also requires them to cross a spatial boundary. Then, al-
though the stars are still scattered to near- and above-escape
energies on the timescale trlx, they require some additional
time to actually leave the cluster. This escape timescale is
related fundamentally to tcr (but also depends on details of
the stellar orbits, the external tidal field, and the shape of
the zero-energy surface). The longer this extra time, the
higher is the probability that further encounters with bound
cluster stars may scatter any potential escapers back down
to sub-escape energies. The net result is a slow-down (“re-
tardation”) of the overall evaporation rate (Chandrasekhar
1942; King 1959; Takahashi & Portegies Zwart 1998, 2000;
Fukushige & Heggie 2000; Baumgardt 2001) and a length-
ening of the cluster lifetime tdis, by a factor that can be ex-
pected to increase with the ratio tcr/trlx. If this factor scales as
(tcr/trlx)
1−x for some x < 1, then tdis ∝ trlx (tcr/trlx)
1−x = txrlxt
While such a retardation of evaporation can be expected
to occur at some level in all clusters, there are physical sub-
tleties in the effect that are probably not captured adequately
by a simple re-parametrization of lifetimes as tdis ∝ t
In particular, it is unlikely that this expression can hold for
clusters of all masses with a single value of x < 1. Since
tcr/trlx ∝ M
−1, very massive clusters have tcr ≪ trlx, and stars
scattered to greater than escape energies by relaxation cross
the tidal boundary effectively instantaneously—implying that
the standard tdis ∝ trlx, or x→ 1, applies in the high-mass limit.
Indeed, if this were not the case, and a fixed x < 1 held for all
M, then an unphysical tdis < trlx would obtain at high enough
masses; see Baumgardt (2001) for further discussion. Unfor-
tunately, “very massive” is not well quantified in this context,
and it is not yet clear if a single value of x is accurate for the
entire GC mass regime. So far, it has been checked directly
only for initial cluster masses below the current peak of the
GCMF.
It is also worth noting that the analysis and simulations
aimed at this problem to date have dealt with clusters on cir-
cular or moderately eccentric orbits in galactic potentials that
are static and spherical. This means that any tidal perturba-
tions felt by stars within the clusters are relatively weak and/or
slow compared to their own orbital periods, leading to nearly
adiabatic or at least non-impulsive responses. In more realis-
tic situations, the galactic potential would be time-dependent
and non-spherical and there might be additional tidal pertur-
bations, including disk and bulge shocks. These perturbations
could in some cases accelerate the escape of weakly bound
stars from the clusters and thus counteract the retardation ef-
fect to some degree. Further study is therefore needed to de-
termine the regime of validity of the formula tdis ∝ t
cr and
its possible modification outside this regime.
In the meantime, Baumgardt (2001) and
Baumgardt & Makino (2003; hereafter BM03) have fit-
ted this formula to the lifetimes of a suite of N-body clusters
with initial masses M0 . 7 × 10
4 M⊙ and several different
initial concentrations and orbital eccentricities. BM03 at
first write tdis in terms of the relaxation and crossing times
of clusters at their half-mass radii, so that trlx ∝ (Mr
tcr ∝ (M/r
−1/2, and tdis ∝ M
x−1/2r
h (see their equation
[5]). However, they immediately take a factor of (rt/rh)
out from the normalization of this scaling—in effect to obtain
tdis ∝M
x−1/2r
t with a different constant of proportionality—
and then use a simple definition of the tidal radius (their
equation [1], r3t = GMr
c , which is appropriate for a
circular orbit of radius rp in a logarithmic potential with
circular speed Vc; see Innanen, Harris, & Webbink 1983)
to obtain the total lifetime of a cluster as a function of its
initial mass, perigalactic distance, and Vc (their equation
[7]). A single exponent x ≃ 0.75 and a single normalization
in this function then suffice to predict to within 10% the
lifetimes of the simulated clusters, regardless of their initial
concentrations. By implication, if trlx and tcr were fixed at rh
rather than rt , then tdis would have an additional concentration
dependence, related to the ratio (rt/rh)
3/2—very similar to
what we discussed in §3.2.1 for the case x = 1.
We now re-examine the Milky Way GCMF in terms of
this prescription for retarded evaporation (bearing in mind the
caveats mentioned above). To avoid any explicit dependences
on concentration, we also focus on the tidal radius and write
tdis ∝ M
x−1/2r
t for general x ≤ 1; but we do not substitute a
potential- and orbit-specific formula for rt in terms of rp and
galactic properties such as Vc. Instead, to keep the empha-
sis entirely on cluster densities, we re-write the scaling of the
lifetime in terms of the mean surface density inside the tidal
radius, Σt ≡ M/πr
t , and the corresponding volume density
ρt = 3M/4πr
t . This leads to tdis ∝ MΣ
−3(1−x)
−2(x−3/4)
t , which
then implies
µev ≡ −dM/dt ∝ M/tdis ∝ Σ
3(1−x)
2(x−3/4)
t . (9)
Clearly, the standard µev ∝ ρ
t , which we have already dis-
cussed, is recovered for x = 1; while for x = 0.75, we have the
equally straightforward µev ∝ Σ
BM03 find that, even with the retarded evaporation implied
by x ≃ 0.75, the masses of their simulated clusters still de-
crease approximately linearly with time after stellar-evolution
effects (which are only important for the first few 108 yr) are
separated out; see especially their Figure 6, equation (12), and
related discussion. Thus, if the GCMF initially rose towards
low masses and has been eroded by slow, relaxation-driven
cluster destruction, then in this modified description of evap-
oration we might expect the current mass function to depend
fundamentally on Σt rather than ρh or ρt . But because M(t)
still decreases nearly linearly with t, only now with µev ∝Σ
for each cluster, the shape of the evolved GCMF and its de-
pendence on Σt should resemble our earlier results for ρh and
SHAPING THE GCMF BY EVAPORATION 13
We have confirmed this expectation by repeating all of our
analyses in §2 again, now using µev ∝Σ
t to estimate cluster
mass-loss rates. As before, we calculate Σt from the data in
the Harris (1996) catalogue, although we caution once more
that the tidal radii, and thus the derived Σt , are more uncertain
than rh and ρh.
Figure 7, which should be compared to Figures 1 and 5
above, shows that the average Galactic GC mass increases
systematically with Σt; that the lower envelope of the M–Σt
distribution is described well by M ∝ Σ
t (the dashed line in
the left-hand panel of Figure 7), which is a locus of constant
lifetime against evaporation for µev ∝ Σ
t ; and that the scat-
ter in the distribution of cluster Σt versus Galactocentric ra-
dius (right-hand panel of the figure) is substantial, as required
to account for the almost non-existent correlation between M
and rgc.
The left-hand side of Figure 8 shows the mass functions
of globulars in three bins of Σt , as defined in each panel.
The right-hand side of the figure shows dN/d log M in the
same three intervals of rgc as in Figures 2 and 6 above. As
in those earlier plots, the dashed curve in all panels of Figure
8 is a model GCMF with the same parameters in every case,
representing the mass function of the entire Galactic GC sys-
tem. Once again, compared to the average MTO, the observed
turnover mass is significantly lower for clusters in the lowest
Σt bin and higher for clusters in the highest Σt bin, while the
width of dN/d log M decreases noticeably as Σt increases.
The solid curves in Figure 8 are again different in every
panel. They are the sums of evaporation-evolved Schechter
functions as in equation (3), with the usual β = 2 assumed
but with total mass losses estimated individually for each GC
in any Σt or rgc bin according to ∆ ∝ Σ
t rather than ∆ ∝
h or ∆ ∝ ρ
t . However, it turns out not to be necessary
to change the normalization of ∆ ∝ ρ
h in equation (4) to
achieve good fits to the observed GCMF as a function of either
Σt or rgc. Thus, in Figure 8 we have simply used
∆ = 1.45× 104 M⊙
Σt/M⊙ pc
−2)3/4 . (10)
The fits of these models, based on tdis ∝ t
cr with x ≃
0.75, are indistinguishable from the fits of our original mod-
els based on the standard tdis ∝ trlx, i.e., x = 1. (We have con-
firmed that adopting individual ∆ given by equation [10] also
reproduces the GCMFs of low-and high-concentration GCs
in Figure 3 as well as before.) It was somewhat unexpected
that equation (10) and equation (4) should have the same nu-
merical coefficient, but we note that this follows empirically
from the fact that the measured ρh and Σt of Galactic GCs
are consistent with the simple near-equality, ρh/M⊙ pc
(Σt/M⊙ pc
−2)1.5 in the mean. This is illustrated in Figure 9,
which also shows that there is significant scatter about the re-
lation.12 However, this scatter does not correlate with clus-
ter mass or Galactocentric radius. From a pragmatic point of
view, therefore, ρ
h and Σ
t are near enough to interchange-
able for our purposes, and there is no practical difference be-
12 Although it may be only a coincidence that the constant of proportion-
ality in ρh ∝ Σ
t is so near unity, the basic scaling itself holds because
combining the observed correlation between cluster mass and central con-
centration (Djorgovski & Meylan 1994; McLaughlin 2000) with the intrinsic
dependence of rt/rh on c in King models leads roughly to (rt/rh) ∝ M
tween GCMF models based on one or the other measure of
GC density.
One further check on this is to verify that the mass-loss rate
associated with equation (10) is roughly in keeping with that
implied by the N-body simulations pointing to x = 0.75 in the
first place. Thus, we compare the rate
µev = ∆/(13 Gyr) ≃ 1100 M⊙ Gyr
−1 (Σt/M⊙ pc
−2)3/4 (11)
to a formula implicit in BM03. Starting with their equa-
tion (7) for the lifetime tdis as a function of initial cluster
mass and perigalactic distance and circular speed in a loga-
rithmic halo potential; using their x = 0.75 and their normal-
ization of 1.91× 106 yr, multiplied as in their equation (9)
by (1 + e) to allow for eccentric orbits with apo- and peri-
galactic distances related by e ≡ (ra − rp)/(ra + rp); insert-
ing their equation (1) for rt ; taking the mean mass of clus-
ter stars to be m∗ = 0.55M⊙, as they do; using γ = 0.02 as
they do in the Coulomb logarithm, ln(γM0/m∗); and defining
Σt,0 ≡ M0/πr
t,0 (the subscript 0 denoting initial values), we
obtain
µev(BM03)≃
0.7M0
1 + e
M⊙ Gyr
0.036M0/M⊙
0.036× 105
]3/4 (
M⊙ pc−2
This is appropriate for clusters that just fill their Roche lobes
at perigalacticon, which is where Σt,0 is specified. The factor
of 0.7 in the first equality accounts for mass loss due to stellar
evolution in the BM03 simulations, which, as they discuss,
can be treated as having occurred almost immediately and in
full at the beginning of a cluster’s life.
Our GCMF-based µev is a factor of ≈ 2 faster than the N-
body value for clusters on circular orbits (with e = 0 and in
steady tidal fields) in the simulations; and our µev is still
within a factor of about three of the N-body rate for clusters on
eccentric orbits with e = 0.5 in BM03 (e ≃ 0.5–0.6 is typical
for tracers with an isotropic velocity distribution in a logarith-
mic potential; van den Bosch et al. 1999). This is very similar
to the comparison of lifetimes in §3.1 for our original models
based on µev ∝ ρ
h . Moreover, our new estimate of µev and
that in BM03 are still subject to their own, separate uncertain-
ties and reflect different idealizations and assumptions. For
example, our rate still depends on the exact power-law expo-
nent β at low masses in the initial GCMF, as discussed after
equation (7); while the rate from BM03 still neglects grav-
itational shocks from disk crossings and passages by a dis-
crete galactic bulge, and may additionally be biased low for
M0 > 10
5 M⊙ if x > 0.75 at such masses. All of this—not
to mention again the large uncertainties and possible system-
atics in the estimates of tidal radii needed to calculate Σt—
makes the near agreement between equations (11) and (12)
more striking than any apparent discrepancy.
In summary, although the relation µev ∝ ρ
h ≃ constant
in time is rigorously correct only in rather specific circum-
stances, our GCMF models based on it in §2 are good proxies,
in all respects, for models based on other plausible characteri-
zations of relaxation-driven cluster mass loss. This result will
likely be important for future studies of the mass functions of
extragalactic cluster systems, where it may well be necessary
14 McLAUGHLIN & FALL
FIG. 7.— Scatter plots of mass M versus mean surface density inside the tidal radius (Σt ≡ M/πr2t ) and of Σt versus Galactocentric radius rgc, for 146 Galactic
GCs from the Harris (1996) catalogue. These plots are analogous to the left- and rightmost panels of Figure 1, and the two panels of Figure 5. The dashed line in
the left-hand plot traces the relation M ∝ Σ
t , which defines a locus of constant evaporation time for µev ∝ Σ
FIG. 8.— Observed GCMF (points, with Poisson errorbars) and models (curves) as a function of mean surface density inside the tidal radius, Σt ≡ M/πr2t (left-
hand panels), and as a function of Galactocentric radius, rgc (right-hand panels). The dashed curve in every panel is an evolved Schechter function representing
the entire GC system: equation (3) with β = 2, Mc = 106 M⊙, and a single ∆, common to all clusters, evaluated from equation (10) using the median bΣt of all
146 Galactic GCs. Solid curves are subsample-specific models using equation (3) with β = 2 and Mc = 106 M⊙ but a different ∆ value for every cluster (obtained
from equation [10] using individual observational estimates of Σt ) in any Σt or rgc bin.
SHAPING THE GCMF BY EVAPORATION 15
FIG. 9.— Half-mass density, ρh = 3M/8πr
h , against mean surface den-
sity inside the tidal radius, Σt = M/πr2t , for 146 clusters with data in Harris
(1996). The straight line is ρh = Σ
to adopt procedures based on ρh rather than ρt or Σt because
of the difficulty or impossibility of estimating tidal radii.
3.3. MTO versus rgc, and Velocity Anisotropy in GC Systems
In this paper we have directly modeled dN/d log M as a
function only of GC density and age, and used the observed ρh
(or ρt , or Σt) of clusters in relatively narrow ranges of Galac-
tocentric position to show that such models are consistent with
the current near-constancy of the GCMF as a function of rgc.
Most other models in the literature for evaporation-dominated
GCMF evolution, in either the Milky Way or other galaxies,
instead predict the distribution explicitly as a function of rgc at
any time. They therefore need, in effect, to derive theoretical
density–position relations for clusters in galaxies alongside
their main GCMF calculations. This usually begins with the
adoption of analytical potentials to describe the parent galax-
ies of GCs. Taking these to be spherical and static for a Hub-
ble time allows the use of standard tidal-limitation formulae
to write GC densities ab initio in terms of the (fixed) peri-
centers rp of unique orbits in the adopted potentials. Cluster
relaxation times and mass-loss rates µev then follow as func-
tions of rp as well. Finally, specific initial mass, space, and
velocity (or orbital eccentricity) distributions are chosen for
entire GC systems, so that at all later times it is known what
the dynamically evolved dN/d log M is for globulars with any
single rp; how many clusters with a given rp survive; and what
the distributions of rp and all dependent cluster properties are
at any instantaneous position rgc.
In this approach, if the GCMF began with a power-law rise
towards low masses and its current peak is due entirely to clus-
ter disruption, then a dependence of MTO on rp is expected in
general, because the densities of tidally limited GCs decrease
with increasing rp. Thus, models along these lines that as-
sume the orbit distribution of a GC system to be the same
at all radii in a galaxy (i.e., that the time average of the ra-
tio rgc/rp is independent of position) have typically had diffi-
culty in accounting for the observed weak or non-correlation
between MTO and present rgc in large galaxies. This is partic-
ularly a problem if it is assumed that the initial GCMF was a
pure power law, with the same index at arbitrarily high masses
as low (e.g., Baumgardt 1998; Vesperini 2001). It is poten-
tially less of a concern if dN/d log M started as a Schechter
function with an exponential cut-off at masses M > Mc, as
we have assumed, since then the existence of a strict upper
bound MTO ≤ Mc (§2.2) means that the dependence of an
evaporation-evolved MTO on rp and rgc must saturate for small
enough galactocentric radii (high enough GC densities). Even
so, the “scale-free” models of FZ01, in which Mc ≃ 10
and all GCs in a Milky Way-like galaxy potential have the
same time-averaged rgc/rp, still predict a gradient in MTO ver-
sus rgc that is stronger than observed.
FZ01 showed that, if they left all of their other assumptions
unchanged, then a dependence of GCMF peak mass on rgc
could be effectively erased by an appropriately varying radial
velocity anisotropy in the initial GC system. Thus, in their
“Eddington” models the eccentricity of a typical cluster or-
bit increases with galactocentric distance (the time average of
rgc/rp increases with radius), such that globulars spread over
a larger range of current rgc can have more similar rp and asso-
ciated MTO. However, the initial velocity-anisotropy gradient
required to fit the Milky Way GCMF data specifically is only
marginally consistent with the observed kinematics of the GC
system (e.g., Dinescu, Girard, & van Altena 1999).13 Subse-
quently, Vesperini et al. (2003) constructed broadly similar
models for the GCMF of the Virgo elliptical M87 and con-
cluded that there, too, a variable radial velocity anisotropy is
required to match the observed MTO versus rgc; but the model
anisotropy profile in this case is clearly inconsistent with the
true velocity distribution of the GC system, which is observed
to be isotropic out to large rgc (Romanowsky & Kochanek
2001; Côté et al. 2001).
These results certainly suggest that some element is lacking
in rgc-oriented GCMF models developed as outlined above.
But they do not mean that the fault lies with the main hypoth-
esis, that the difference between the mass functions of young
clusters and old GCs is due to the effects of slow, relaxation-
driven disruption in the latter case. Any conclusions about
velocity anisotropy depend on the totality of steps taken to
connect the densities and positions of clusters; and it is possi-
ble that reasonable changes to one or more of these ancillary
assumptions could make the models compatible with the ob-
served kinematics of GCs in both the Milky Way and M87,
without abandoning a basic physical picture of evaporation-
dominated GCMF evolution that is otherwise quite success-
One issue is that previous models have always specified
evaporation rates a priori as functions of cluster density (or or-
bital pericenter), usually normalizing µev so that tdis/trh ≃ 20–
40 as in standard treatments of two-body relaxation. How-
ever, following our discussion in §3.1 and §3.2, it would seem
worthwhile to investigate these models with µev increased at
fixed ρh or rp to allow tdis/trh ≈ 10 (if β ≃ 2 for the low-mass
power-law part of the initial GCMF).
FZ01 and Vesperini et al. (2003) both consider velocity dis-
tributions parametrized by a galactocentric anisotropy radius,
RA, inside of which a cluster system is essentially isotropic
and beyond which it is increasingly dominated by radial or-
13 The fact that clusters on radial orbits are preferentially disrupted lessens
any inconsistency between the radial anisotropy required in the initial veloc-
ity distribution and observational constraints on the present velocity distribu-
tion.
16 McLAUGHLIN & FALL
bits. In these terms, the difficulty with the published models
is that, to reproduce the observed insensitivity of MTO to rgc
given standard normalizations of µev, they require values of
RA that are smaller than allowed by observations (especially
for M87). Increasing RA to more realistic values while keep-
ing the normalization of µev fixed leads to a stronger gra-
dient in MTO: the orbits of GCs at small rgc . RA remain
closely isotropic and the typical rp and MTO are essentially
unchanged, while at large galactocentric distances the clus-
ter orbits are on average less radial than before, with larger
rp, lower densities, and lower evolved MTO for a given rgc.
This effect is illustrated, for example, in Figure 9 of FZ01.
However, it can be compensated at least in part by increas-
ing µev by a common factor for all GCs, with the new, larger
RA fixed, if the initial mass function is assumed to have been
a Schechter function rather than a pure power law extend-
ing to arbitrarily high masses. A faster evaporation rate
will then lead to a (roughly) proportionate increase in the
evolved GCMF peak mass for GCs with relatively low den-
sities, i.e., those at large rgc and rp; but the increase in MTO
will be smaller, and eventually even negligible, for higher-
density clusters at progressively smaller rgc—again because
MTO grows less than linearly with µev ∝ ρ
h when there is
an upper limit MTO < Mc due to an exponential cut-off in the
initial dN/d log M0. Thus, the qualitative effect of increasing
the normalization of µev in models with radially varying GC
velocity anisotropy is to weaken the amount of radial-orbit
bias required to fit an observed MTO versus rgc.
Another point, emphasized by FZ01, has to do with the
standard starting assumption that GCs orbit in galaxies that
are perfectly static and spherical. In reality, galaxies grow
hierarchically. In this case, even if the values of µev are not
changed, much of the burden for the weakening or erasing of
any initial gradients in MTO versus rgc may be transferred from
velocity anisotropy to the time-dependent evolution of the
galaxies themselves. Violent relaxation, major mergers, and
smaller accretion events all work to move clusters between
different parts of galaxies and between different progenitors,
scrambling and combining any number of pericenter–density–
MTO relations. Any position dependences in the GC ρh dis-
tribution and in MTO itself for the final galaxy are therefore
bound to be weaker, more scattered, and more difficult to re-
late accurately to a cluster velocity distribution than in the
case of a monolithic, non-evolving potential. Allowing for
a non-spherical galaxy potential would have qualitatively the
same effect, because in this case every cluster explores a range
of pericenters and different maximum tidal fields on each of
its orbits.
In this situation, it may be important to ask how evapora-
tion rates can still be approximately constant in time—so that
cluster masses still decrease approximately linearly with t as
our models assume—if the tidal field around any given GC
changes significantly over time. Thus, consider first a sys-
tem of GCs in a single, static galaxy potential. The mass-
evolution curve for each cluster is approximately a straight
line, M(t) ≃ M0 −µevt, with µev depending on some measure
of internal density, which may be ρ
h , ρ
t , or Σ
t . The av-
erage mass-evolution curve for the entire system of clusters
is also approximately linear, 〈M(t)〉 ≃ 〈M0〉 − 〈µev〉t. If now
a merger or other event rearranges the clusters in the galaxy,
then after the event the mass-loss rates of some clusters will
be higher than before and the rates of other clusters will be
lower than before. However, if the mean density of the galaxy
as a whole is roughly the same after the event as before, then
so too will be the average of the GC densities, because of
tidal limitation. The average 〈µev〉 ∝ 〈ρ
h 〉 (say) will differ
even less between the pre- and post-merger systems. Thus, al-
though using instantaneous densities to estimate the past µev
of individual clusters may err on the high side for some clus-
ters and on the low side for others, these errors will average
away to a small or even zero net bias. The approximation
µev ≃ constant in time in our GCMF models will then still be
valid in the mean, and the average 〈M(t)〉 dependence of suf-
ficiently large numbers of clusters will remain roughly linear.
This type of scenario might be expected to pertain at least
to galaxies that evolve on the fundamental plane, since this
entails a connection between the total (baryonic plus dark)
masses and circular speeds of galaxies, of the form Mgal ∝V
or Mgal ∝ V
c . By the virial theorem, the average densities
scale as ρgal ∝ V
2, and thus ρgal ∝ M
gal or ρgal ∝ M
gal .
Insofar as 〈ρh〉 ∝ 〈ρt〉 ∝ ρgal for the GCs, the system-wide av-
erage 〈µev〉 ∝ 〈ρ
h 〉 should therefore not change drastically
even after a major merger between two fundamental-plane
galaxies; at most, the ratio of final to initial 〈µev〉 will be
roughly of order the −1/4 power of the ratio of final to ini-
tial Mgal. Note that this line of reasoning is closely related to
that applied by FZ01 to explain the small observed galaxy-
to-galaxy differences in the average turnover masses of entire
GC systems (although non-zero differences do exist, and can
be accomodated in these sorts of arguments; see Jordán et al.
2006, 2007).
A full exploration of questions such as these, about the wide
range of ingredients in current GC-plus-galaxy models, will
most likely require large N-body simulations set in a realistic,
cold dark matter cosmology. Until these can be carried out,
it is our view that the kinematics of globular cluster systems
cannot be used as decisive side constraints on theories for the
GCMF.
4. CONCLUSIONS
We have shown that the mass function dN/d log M of glob-
ular clusters in the Milky Way depends significantly on clus-
ter half-mass density, ρh, with the peak or turnover mass MTO
increasing and the width of the distribution decreasing as ρh
increases. This behavior is expected if the GCMF initially
rose towards masses below the present turnover scale—as the
mass functions of young cluster systems like that in the An-
tennae galaxies do—and has evolved to its current shape via
the slow depletion of low-mass clusters over Gyr timescales,
primarily through relaxation-driven evaporation. The fact that
MTO increases with cluster density favors evaporation over
external gravitational shocks as the primary mechanism of
low-mass cluster disruption, since the mass-loss rates asso-
ciated with shocks depend inversely on cluster density and
directly on cluster mass. Our results therefore add to previ-
ous arguments supporting an interpretation of the GCMF in
terms of evaporation-dominated evolution, based on the fact
that dN/d log M scales as M1−β with β ≃ 0 in the low-mass
limit (Fall & Zhang 2001).
The observed GCMF as a function of ρh is fitted well
by simple models in which the initial distribution was
a Schechter function, dN/d log M0 ∝ M
0 exp
−M0/Mc
with β = 2 and Mc ≃ 10
6 M⊙ assumed, and in which clusters
have been losing mass for a Hubble time at roughly steady
rates that can be estimated from their current half-mass den-
SHAPING THE GCMF BY EVAPORATION 17
sities as µev ∝ ρ
h . We have shown that, although this pre-
scription is approximate, it captures the main physical depen-
dence of relaxation-driven evaporation. In particular, it leads
to model GCMFs that are entirely consistent with those re-
sulting from alternative characterizations of evaporation rates
in terms of cluster tidal densities ρt or mean surface densities
Σt (§3.2). The normalization of µev at a given ρh (or ρt , or Σt)
required to fit the GCMF implies total cluster lifetimes that
are within range of the lifetimes typically obtained in theoret-
ical studies of two-body relaxation, although our values may
be slightly shorter than the theoretical ones if the low-mass,
power-law part of the initial cluster mass function was as steep
as we have assumed.
Taking clusters in various bins of central concentration c
and Galactocentric radius rgc and using their (individual) ob-
served densities as direct input to our models yields dynam-
ically evolved GCMFs as functions of c and rgc that agree
well with all data. This again indicates that the most fun-
damental physical dependence in the GCMF is that on clus-
ter density. Moreover, our models for dN/d log M versus
rgc obtained in this way are consistent in particular with the
well-known insensitivity of the GCMF peak mass to Galac-
tocentric position. This is seen to follow from a significant
variation of MTO with ρh (or ρt , or Σt)—due in our analysis
to evaporation-dominated cluster disruption—combined with
substantial scatter in the GC densities at any Galactocentric
position.
We have not invoked an anisotropic GC velocity distribu-
tion to explain the observed weak variation of MTO with rgc;
indeed, we have made no predictions or assumptions what-
soever about velocity anisotropy. We have emphasized that,
when velocity anisotropy enters other long-term dynamical-
evolution models for the GCMF, it is only in conjunction with
several additional, interrelated assumptions made as part of
larger efforts to derive theoretical density–rgc relations for
GCs—which we have not attempted to do here. The appar-
ent need in some current models for a strong bias towards
high-eccentricity cluster orbits to explain the near-constancy
of MTO versus rgc might well be avoided by changing one or
more ancillary assumptions in the models, without having to
discard the underlying idea that the peak and low-mass shape
of the GCMF are the result of relaxation-driven cluster dis-
ruption.
It clearly will be of interest to test and refine the main
ideas in this paper through modeling of the GCMFs in other
galaxies. For the time being at least, doing so will re-
quire the estimation of approximate mass-loss rates using
cluster half-mass densities rather than tidal quantities, sim-
ply because GC half-light radii can be measured accurately
in many systems beyond the Local Group, whereas tidal
radii are much more model-dependent and difficult to ob-
serve. Chandar, Fall, & McLaughlin (2007) have recently
shown that the peak mass of the GCMF in the Sombrero
galaxy (M104) increases with ρh in a way that is reasonably
well described by sums of evolved Schechter (1976) functions
as in the models presented in this paper. It should be rela-
tively straightforward to pursue similar studies in other nearby
galaxies.
We thank Michele Trenti, Douglas Heggie, Bill Harris, Ru-
pali Chandar, and Bruce Elmegreen for helpful discussions
and comments. SMF acknowledges support from the Am-
brose Monell Foundation and from NASA grant AR-09539.1-
A, awarded by the Space Telescope Science Institute, which is
operated by AURA, Inc., under NASA contract NAS5-26555.
REFERENCES
Aguilar, L., Hut, P., & Ostriker, J. P. 1988, ApJ, 335, 720
Barmby, P., Huchra, J. P., & Brodie, J. P. 2001, AJ, 121, 1482
Barmby, P., McLaughlin, D. E., Harris, W. E., Harris, G. L. H., & Forbes, D.
A. 2007, AJ, 133, 2764
Baumgardt, H. 1998, A&A, 330, 480
Baumgardt, H. 2001, MNRAS, 325, 1323
Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227 (BM03)
Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton
University Press)
Burkert, A., & Smith, G. H. 2000, ApJ, 542, L95
Caputo, F., & Castellani, V. 1984, MNRAS, 207, 185
Chandar, R., Fall, S. M., & McLaughlin, D. E. 2007, ApJ, 668, L119
Chandrasekhar, S. 1942, Principles of Stellar Dynamics (Chicago:
University of Chicago Press)
Chernoff, D. F., & Weinberg, M. D. 1990, ApJ, 351, 121
Côté, P., et al. 2001, ApJ, 559, 828
Dinescu, D. I., Girard, T. M., & van Altena, W. F. 1999, AJ, 117, 1792
Djorgovski, S., & Meylan, G. 1994, AJ, 108, 1292
Elmegreen, B. G., & Efremov, Y. N. 1997, ApJ, 480, 235
Fall, S. M., & Rees, M. J. 1977, MNRAS, 181, 37P
Fall, S. M., & Zhang, Q. 2001, ApJ, 561, 751 (FZ01)
Fukushige, T., & Heggie, D. C. 2000, MNRAS, 318, 753
Giersz, M. 2001, MNRAS, 324, 218
Giersz, M., & Heggie, D. C. 1996, MNRAS, 279, 1037
Gnedin, O. Y. 1997, ApJ, 487, 663
Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223
Gnedin, O. Y., Lee, H. M., & Ostriker, J. P. 1999, ApJ, 522, 935
Harris, W. E. 1996, AJ, 112, 1487
Harris, W.E. 2001, in Star Clusters (28th Saas-Fee Advanced Course) ed. L.
Labhardt & B. Binggeli (Berlin: Springer), 223
Harris, W. E., & Pudritz, R. E. 1994, ApJ, 429, 177
Harris, W. E., Harris, G. L. H., & McLaughlin, D. E. 1998, AJ, 115, 1801
Hénon, M. 1961, Ann. d’Astrophys., 24, 369
Innanen, K. A., Harris, W. E., & Webbink, R. F. 1983, AJ, 88, 338
Johnstone, D. 1993, AJ, 105, 155
Jordán, A., et al. 2005, ApJ, 634, 1002
Jordán, A., et al. 2006, ApJ, 651, L25
Jordán, A., et al. 2007, ApJS, 171, 101
Joshi, K. J., Nave, C. P., & Rasio, F. A. 2001, ApJ, 550, 691
Kavelaars, J. J., & Hanes, D. A. 1997, MNRAS, 285, L31
Lee, H. M., & Goodman, J. 1995, ApJ, 443, 109
King, I. 1958, AJ, 63, 109
King, I. 1959, AJ, 64, 351
King, I. R. 1966, AJ, 71, 64
Lee, H. M., & Ostriker, J. P. 1987, ApJ, 322, 123
McLaughlin, D. E. 2000, ApJ, 539, 618
McLaughlin, D. E., & van der Marel, R. P. 2005, ApJS, 161, 304
Murali, C., & Weinberg, M. D. 1997, MNRAS, 291, 717
Okazaki, T., & Tosa, M. 1995, MNRAS, 274, 48
Ostriker, J. P., & Gnedin, O. Y. 1997, ApJ, 487, 667
Parmentier, G., & Gilmore, G. 2007, MNRAS, 377, 352
Prieto, J. L., & Gnedin, O. Y. 2006, preprint (astro-ph/0608069)
Romanowsky, A. J., & Kochanek, C. S. 2001, ApJ, 553, 722
Schechter, P. 1976, ApJ, 203, 297
Smith, G. H., & Burkert, A. 2002, ApJ, 578, L51
Spitler, L. R., Larsen, S. S., Strader, J., Brodie, J. P., Forbes, D. A., &
Beasley, M. A. 2006, AJ, 132, 1593
Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton:
Princeton Univ. Press)
Takahashi, K., & Portegies Zwart, S. F. 1998, ApJ, 503, L49
Takahashi, K., & Portegies Zwart, S. F. 2000, ApJ, 535, 759
Trenti, M., Heggie, D. C., & Hut, P. 2007, MNRAS, 374, 344
van den Bosch, F. C., Lewis, G. F., Lake, G., & Stadel, J. 1999, ApJ, 515, 50
Vesperini, E. 1997, MNRAS, 287, 915
Vesperini, E. 1998, MNRAS, 299, 1019
Vesperini, E. 2000, MNRAS, 318, 841
Vesperini, E. 2001, MNRAS, 322, 247
Vesperini, E., & Heggie, D. C. 1997, MNRAS, 289, 898
Vesperini, E., & Zepf, S. E. 2003, ApJ, 587, L97
Vesperini, E., Zepf, S. E., Kundu, A., & Ashman, K. M. 2003, ApJ, 593, 760
Waters, C. Z., Zepf, S. E., Lauer, T. R., Baltz, E. A., & Silk, J. 2006, ApJ,
650, 885
Zhang, Q., & Fall, S. M. 1999, ApJ, 527, L81
http://arxiv.org/abs/astro-ph/0608069
