Shaping the Globular Cluster Mass Function by Stellar-Dynamical Evaporation

TL;DR

This paper demonstrates that the globular cluster mass function's shape and turnover mass depend on cluster density, with evaporation driven by internal relaxation explaining observed variations without assuming universality.

Contribution

It introduces a simple model linking cluster evaporation rates to density, accurately reproducing observed GCMF features and their dependence on cluster properties.

Findings

GCMF turnover mass increases with cluster density

Relaxation-driven evaporation explains GCMF shape and density dependence

Cluster lifetimes inferred are within standard ranges

Abstract

We show that the globular cluster mass function (GCMF) in the Milky Way depends on cluster half-mass density (rho_h) in the sense that the turnover mass M_TO increases with rho_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 rho_h, central concentration c, and Galactocentric distance r_gc -- and a simple model in which the relaxation-driven mass-loss rates of clusters are approximated by -dM/dt = mu_ev ~ rho_h^{1/2}. In particular, we recover the well-known insensitivity of M_TO to r_gc. This feature does not derive from a literal ``universality'' of the GCMF turnover mass, but rather from a significant variation of M_TO with rho_h -- the expected outcome of relaxation-driven cluster disruption -- plus significant scatter in rho_h as a function of r_gc. 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 mu_ev ~ rho_t^{1/2} or mu_ev ~ Sigma_t^{3/4} -- alternative prescriptions that are physically motivated but involve cluster properties (rho_t and Sigma_t) that are not as well defined or as readily observable as rho_h. In all cases, the normalization of mu_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.

Figures

Tables

Table 1. Table 1 Milky Way GC Properties in Bins of Density and Galactocentric Radius

Bin	$𝒩$	${\hat{\rho }}_h$ aaThe notation $\hat{x}$ represents the median of quantity $x$.	${\hat{r}}_{\mathrm{gc}}$ aaThe notation $\hat{x}$ represents the median of quantity $x$.	$M_{\min }$	$M_{\max }$	$\hat{M}$ aaThe notation $\hat{x}$ represents the median of quantity $x$.	$M_{\mathrm{TO}}$ bb$M_{\mathrm{TO}}$ 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 $\beta =2$, $M_c={10}^6{M}_{\odot }$, and individual $\Delta$ given by the observed ${\rho }_h$ of each cluster through equation (4).
		[$M_{\odot }{\mathrm{pc}}^{-3}$]	[kpc]	[$M_{\odot }$]	[$M_{\odot }$]	[$M_{\odot }$]	[$M_{\odot }$]
${\rho }_h$  bins
$0.034\le {\rho }_h\le 76.5M_{\odot }{\mathrm{pc}}^{-3}$	48	8.48	12.9	$5.63\times {10}^2$	$8.84\times {10}^5$	$4.12\times {10}^4$	$3.98\times {10}^4$
$78.8\le {\rho }_h\le 526M_{\odot }{\mathrm{pc}}^{-3}$	49	232	5.6	$8.37\times {10}^3$	$1.67\times {10}^6$	$1.22\times {10}^5$	$1.58\times {10}^5$
$579\le {\rho }_h\le 5.65\times {10}^4{M}_{\odot }{\mathrm{pc}}^{-3}$	49	973	3.2	$1.93\times {10}^4$	$1.30\times {10}^6$	$2.82\times {10}^5$	$2.88\times {10}^5$
$r_{\mathrm{gc}}$  bins
$0.6\le r_{\mathrm{gc}}\le 3.2$ kpc	47	597	1.9	$4.47\times {10}^3$	$1.02\times {10}^6$	$1.15\times {10}^5$	$2.14\times {10}^5$
$3.3\le r_{\mathrm{gc}}\le 9.4$ kpc	50	261	5.2	$2.02\times {10}^3$	$1.67\times {10}^6$	$1.27\times {10}^5$	$1.66\times {10}^5$
$9.6\le r_{\mathrm{gc}}\le 123$ kpc	49	18.4	18.3	$5.63\times {10}^2$	$1.30\times {10}^6$	$7.42\times {10}^4$	$8.71\times {10}^4$

Equations

$$\frac{dN}{d\log M}=\frac{M}{M_0}\times \frac{dN}{d\log M_0}=\frac{M}{(M+\Delta )}\frac{dN}{d\log \left(M+\Delta \right)}.$$

$$dN/d\log M_0\propto M_0^{1-\beta }\exp \left(-M_0/M_c\right).$$

$$\frac{dN}{d\log M}=\sum _{i=1}^{\mathcal{N}}\frac{A_{i}M}{{\left[M+{\Delta }_i\right]}^{\beta }}\exp \left[-\frac{M+{\Delta }_i}{M_c}\right].$$

$$\Delta =1.45\times 10^4{M}_{\odot }{\left({\rho }_h/M_{\odot }{\mathrm{pc}}^{-3}\right)}^{1/2}$$

$${\mu }_{\mathrm{ev}}\simeq 1100M_{\odot }{\mathrm{Gyr}}^{-1}{\left({\rho }_h/M_{\odot }{\mathrm{pc}}^{-3}\right)}^{1/2}.$$

$$\frac{t_{\mathrm{dis}}}{t_{\mathrm{rh}}}=\frac{M}{{\mu }_{\mathrm{ev}}{t}_{\mathrm{rh}}}\simeq 10\left[\frac{\ln \left(0.57M/M_{\odot }\right)}{\ln \left(0.57\times 10^5\right)}\right].$$

$$t_{\mathrm{dis}}/t_{\mathrm{rh}}\propto {\widehat{{\rm Y}}}_V^{-1/2}(\beta -1{)}^{-1}.$$

$$\Delta =2.1\times 10^5{M}_{\odot }{\left({\rho }_t/M_{\odot }{\mathrm{pc}}^{-3}\right)}^{1/2}.$$

$${\mu }_{\mathrm{ev}}\equiv -dM/dt\propto M/t_{\mathrm{dis}}\propto {\Sigma }_t^{3(1-x)}{\rho }_t^{2(x-3/4)}.$$

$$\Delta =1.45\times 10^4{M}_{\odot }{\left({\Sigma }_t/M_{\odot }{\mathrm{pc}}^{-2}\right)}^{3/4}.$$

$${\mu }_{\mathrm{ev}}=\Delta /(13\mathrm{Gyr})\simeq 1100M_{\odot }{\mathrm{Gyr}}^{-1}({\Sigma }_t/M_{\odot }{\mathrm{pc}}^{-2}{)}^{3/4}$$

$${\mu }_{\mathrm{ev}}(BM03)$$

$$\times {\left[\frac{\ln \left(0.036M_0/M_{\odot }\right)}{\ln \left(0.036\times 10^5\right)}\right]}^{3/4}{\left(\frac{{\Sigma }_{t,0}}{M_{\odot }{\mathrm{pc}}^{-2}}\right)}^{3/4}.$$

Full text

††slugcomment: The Astrophysical Journal, 679:1272–1287, 2008 June 1

1. Introduction

The mass functions of star cluster systems provide an important point of reference for attempts to understand the connection between old globular clusters (GCs) and the young massive clusters that form in local starbursts and galaxy mergers. When expressed as the number per unit logarithmic mass, $dN/d\,\log\,M$, the GC mass function (GCMF) is characterized by a peak, or turnover, at a mass $M_{\rm TO}\approx 1$–$2\times 10^{5}\,M_{\odot}$ that is empirically very similar in most galaxies. By contrast, the mass functions of young clusters show no such feature but instead rise monotonically towards low masses over the full observed range ($10^{6}\,M_{\odot}\gtrsim M\gtrsim 10^{4}\,M_{\odot}$ in the best-studied cases), in a way that is well described by a power law, $dN/d\,\log\,M\propto M^{1-\beta}$ with $\beta\simeq 2$ (e.g., Zhang & Fall, 1999).

At the same time, for high $M>M_{\rm TO}$, 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 complete 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 depletion over many Gyr of globulars from initial mass distributions that were similar to those of young clusters below $M_{\rm TO}$ as well as above.

Our chief purpose in this paper is to establish and interpret an aspect of the Galactic GCMF that appears fundamental but has gone largely unnoticed to date: $dN/d\,\log\,M$ has a strong and systematic dependence on GC half-mass density, $\rho_{h}\equiv 3M/8\pi r_{h}^{3}$ ($r_{h}$ being the cluster half-mass radius), in the sense that the turnover mass $M_{\rm TO}$ increases and the width of the distribution decreases with increasing $\rho_{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 dynamical evolution from a mass function that initially increased towards $M<M_{\rm TO}$, if the long-term mass loss from surviving 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 evolution removes (through supernovae and winds) the same fraction of mass from all clusters of a given age, and so cannot change the shape of $dN/d\,\log\,M$ (unless special initial conditions 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 during disk crossings and bulge passages is generally less than that due to evaporation for $M<M_{\rm TO}$ (FZ01; Prieto & Gnedin 2006).111It is possible that there existed a past population of GCs with low densities 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.

As we discuss further in §2 below, the evaporation of tidally limited clusters proceeds at a rate, $\mu_{\rm ev}\equiv-dM/dt$, that is approximately 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\propto M^{1-\beta}$ with $\beta\rightarrow 0$ (corresponding to a flat distribution of clusters per unit linear mass) at sufficiently low $M<\mu_{\rm ev}t$ in the evolved mass function of coeval GCs that began with any nontrivial initial $dN/d\,\log\,M_{0}$.222Throughout this paper, we use “initial” to mean at a relatively early time in the development of long-lived clusters, after they have dispersed any remnants 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. To accommodate this when $dN/d\,\log\,M_{0}$ originally increased towards low masses as a power law, a time-dependent peak must develop in the GCMF at a mass of order $M_{\rm TO}(t)\sim\mu_{\rm ev}t$ (FZ01). But then, since $\mu_{\rm ev}$ depends fundamentally on cluster density, so too must $M_{\rm TO}$.

A $\beta\simeq 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 $\mu_{\rm 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 $M_{\rm TO}$ generally differs very little among galaxies (e.g., Harris, 2001; Jordán et al., 2006).

If this picture is basically correct, it implies that, even though $M_{\rm TO}$ 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 different densities should have different turnovers. In §2, we show—working for definiteness and relatively easy observability with the half-mass density, $\rho_{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 $\rho_{h}$ with models assuming that (1) the initial distribution increased as a $\beta=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 $\mu_{\rm ev}\propto\rho_{h}^{1/2}$. In §3 we discuss the validity of this prescription for $\mu_{\rm ev}$, which is certainly approximate but captures the main physical dependence of relaxation-driven mass loss. In particular, we show that the alternative mass-loss laws $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ and $\mu_{\rm ev}\propto\Sigma_{t}^{3/4}$—where $\rho_{t}$ and $\Sigma_{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 $\mu_{\rm ev}\propto\rho_{h}^{1/2}$. The normalization of $\mu_{\rm ev}$ required to fit the observed GCMF implies cluster lifetimes that are within a factor of $\approx\!2$ (perhaps slightly on the low side, if the initial power-law exponent at low masses was $\beta=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 individual clusters are used in our models to predict GCMFs in different bins of Galactocentric radius ($r_{\rm gc}$), they fit the much weaker variation of $dN/d\,\log\,M$ and $M_{\rm TO}$ as functions of $r_{\rm gc}$, 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 measured $\rho_{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 density, 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 $M_{\rm TO}$ with $r_{\rm gc}$ can be achieved only in GC systems with strongly radially anisotropic velocity distributions, which are not observed (see especially Vesperini et al. 2003). This apparent inconsistency has been cited to bolster some recent attempts to identify a mechanism by which a “universal” peak at $M_{\rm TO}\sim 10^{5}\,M_{\odot}$ might have been imprinted 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 successes of an evaporation-dominated evolutionary scenario for the origin of $M_{\rm TO}$, 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 models must now be reconsidered in light of the non-universality of $M_{\rm TO}$ 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 $\rho_{h}$ versus $r_{\rm gc}$ in the Milky Way as a given, to show consistency with the observed behavior of $M_{\rm TO}$ as a function of $r_{\rm gc}$. Most other models (FZ01; Vesperini et al. 2003; and references therein) predict dynamically evolved GCMFs directly in terms of $r_{\rm gc}$, and in doing so are forced also to derive theoretical dependences of cluster density on $r_{\rm gc}$. It is only at this stage that GC orbital distributions 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 distributions 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 evolution.

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, $\mu_{\rm ev}\equiv-dM/dt\simeq{\rm constant}$ in time. Thus, the total cluster mass decreases approximately linearly, as $M(t)\simeq M_{0}-\mu_{\rm ev}t$. 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 calculations (e.g., Lee & Ostriker 1987; Chernoff & Weinberg 1990; 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 without external gravitational shocks) and with different internal processes and ingredients (with or without stellar mass spectra, 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\gtrsim M(t)/M_{0}\gtrsim 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 $\mu_{\rm 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 $\Delta\equiv\mu_{\rm ev}t$ is then associated naturally with any system of coeval clusters having a common mass-loss rate: all those with initial $M_{0}\leq\Delta$ are disrupted by time $t$, and replaced with the remnants of objects that began with $M_{0}>\Delta$. As we mentioned in §1, if the initial GCMF increased towards low masses as a power law, then $\Delta$ is closely related to a peak in the evolved distribution, which eventually decreases towards low $M<\Delta$ as $dN/d\,\log\,M\propto M^{1-\beta}$ with $\beta=0$ (FZ01).

In standard theory (e.g., Spitzer 1987; Binney & Tremaine 1987, Section 8.3), the lifetime of a cluster against evaporation is a multiple of its two-body relaxation time, $t_{\rm rlx}$. For a total mass $M$ of stars within a radius $r$, this scales to first order (ignoring a weak mass dependence in the Coulomb logarithm) as $t_{\rm rlx}(r)\propto(Mr^{3})^{1/2}\propto M/\rho^{1/2}$, where $\rho\propto M/r^{3}$. In a concentrated cluster with an internal density gradient, $t_{\rm rlx}(r)$ of course varies throughout the cluster, and the global relaxation timescale is an average of the local values (see the early discussion by King 1958). This can still be written as $t_{\rm rlx}\propto M/\overline{\rho}^{1/2}$, with $M$ the total cluster mass and $\overline{\rho}$ an appropriate reference density. We then have for the instantaneous mass-loss rate, $\mu_{\rm ev}\equiv-dM/dt\propto M/t_{\rm rlx}\propto\overline{\rho}^{1/2}$. Insofar as this is approximately constant in time, a GCMF evolving from an initial $\beta>1$ power law at low masses should therefore develop a peak at a mass that depends on cluster density and age through the parameter $\Delta\propto\overline{\rho}^{1/2}t$.

It remains to identify the best measure of $\overline{\rho}$ in this context. A standard choice in the literature, and the one that we eventually make to derive our main results in this paper, is the half-mass density $\rho_{h}=3M/8\pi r_{h}^{3}$. However, in a steady tidal field, the mean density $\rho_{t}$ inside the tidal radius of a cluster is constant by definition, and thus choosing $\overline{\rho}=\rho_{t}$ instead is the simplest way to ensure that $\mu_{\rm ev}\propto\overline{\rho}^{1/2}$ and $\mu_{\rm ev}\simeq{\rm 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 $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ 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 $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ is routinely used to set the GC mass-loss rates in models for the dynamical evolution of the GCMF, although such studies normally express $\mu_{\rm ev}$ immediately in terms of orbital pericenters, $r_{p}$, most often by assuming $\rho_{t}\propto r_{p}^{-2}$ as for GCs in galaxies whose total mass distributions follow 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 paper. 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 difficult if not impossible to measure in distant galaxies. We deal with this here by focusing on the GCMF as a function of cluster density $\rho_{h}$ inside the less ambiguous, empirically better determined and more robust half-mass radius, asking how simple models with $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ fare against the data.

Taking $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ in place of $\mu_{\rm ev}\propto\rho_{t}^{1/2}$, which we do to construct evaporation-evolved model GCMFs in §2.2, is most appropriate if the ratio $\rho_{t}/\rho_{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 justified. However, real clusters are not homologous ($\rho_{t}/\rho_{h}$ differs among clusters) and they do not evolve self-similarly ($\rho_{h}$ may vary in time even if $\rho_{t}$ does not). The key assumption in our models is that $\mu_{\rm ev}$ is approximately independent of time for any GC, which is well-founded in any case. By using current $\rho_{h}$ values to estimate $\mu_{\rm 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 $(\rho_{t}/\rho_{h})^{1/2}$. Equivalently, we ignore a dependence on cluster concentration in the normalization of $\mu_{\rm ev}\propto\rho_{h}^{1/2}$. As we discuss further in §3, it is reasonable to neglect this complication in a first approximation because $(\rho_{t}/\rho_{h})^{1/2}$ varies much less among Galactic globulars than $\rho_{t}$ and $\rho_{h}$ do separately. We demonstrate this explicitly by repeating our analysis with $\rho_{h}$ replaced by $\rho_{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 $t_{\rm rlx}$ (Baumgardt 2001; Baumgardt & Makino 2003). We point out that this implies that $\mu_{\rm ev}$ may increase as a modest power of the average surface 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),333Feb. 2003 version; see http://physwww.mcmaster.ca/$\sim$harris/mwgc.dat  . along with the distribution of $\rho_{h}$ versus $r_{\rm gc}$ linking the two mass plots. The Harris catalogue actually records the absolute $V$ magnitudes of the GCs. We obtain masses from these by applying the population-synthesis model mass-to-light ratios $\Upsilon_{V}$ computed by McLaughlin & van der Marel (2005) for individual clusters based on their metallicities and an assumed age of 13 Gyr. However, we first multiplied all of the McLaughlin & van der Marel $\Upsilon_{V}$ values by a factor of 0.8 so as to obtain a median $\widehat{\Upsilon}_{V}\simeq 1.5\,M_{\odot}\,L_{\odot}^{-1}$ in the end,444Throughout this paper, we use $\widehat{x}$ to denote the median of any quantity $x$. 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 different metallicities. Our application of a corrective factor to the population-synthesis values, $\Upsilon_{V}^{\rm pop}$, is motivated empirically by the fact that their distribution among Galactic GCs is strongly peaked around a median $\widehat{\Upsilon}_{V}^{\rm pop}\simeq 1.9\ M_{\odot}\,L_{\odot}^{-1}$, while the observed (dynamical) $\Upsilon_{V}^{\rm dyn}$ lie in a fairly narrow range around $\widehat{\Upsilon}_{V}^{\rm dyn}\simeq 1.5\ M_{\odot}\,L_{\odot}^{-1}$ (McLaughlin & van der Marel, 2005). However, it is worth noting that the size of this difference 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, $\Upsilon_{V}^{\rm dyn}$ falls below $\Upsilon_{V}^{\rm pop}$ due to the preferential 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 $\widehat{\Upsilon}_{V}^{\rm dyn}<\widehat{\Upsilon}_{V}^{\rm pop}$ 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 systematically lower ratios of $\Upsilon_{V}^{\rm dyn}/\Upsilon_{V}^{\rm pop}$. 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 correlation between $\Upsilon_{V}^{\rm dyn}/\Upsilon_{V}^{\rm pop}$ and $t_{\rm rh}$ for Galactic globulars). We therefore proceed, as stated, with a single $\Upsilon_{V}^{\rm dyn}/\Upsilon_{V}^{\rm pop}=0.8$ assumed for all GCs.

Harris (1996) gives the projected half-light radius $R_{h}$ 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 $r_{h}=(4/3)R_{h}$ (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 $r_{h}$ equal to the median value for those of the other 141 GCs having masses within a factor two of the one with unknown $r_{h}$. In all cases, the half-mass density is $\rho_{h}\equiv 3M/8\pi r_{h}^{3}$.

The leftmost panel in Figure 1 shows immediately that the cluster mass distribution has a strong dependence on half-mass density: the median $\widehat{M}$ increases with $\rho_{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 $r_{h}$ 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 $\rho_{h}$, is behind the finding (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 $\rho_{h}$ as we have done here, is that any apparent correlation might only be a trivial reflection of the definition $\rho_{h}\propto M/r_{h}^{3}$. 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$–$r_{h}$ correlation does not imply that all GCs have the same $r_{h}$, even within the unavoidable measurement errors. The root-mean-square (rms) scatter of $\log\,r_{h}$ about its average value is $\pm 0.3$ for Galactic GCs, and the 68-percentile spread in $\log\,r_{h}$ 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, $\chi^{2}$ fitting uncertainties) of $\delta(\log\,r_{h})\approx 0.05$, or about 10% relative error; and an rms systematic measurement error (i.e., differences in the $r_{h}$ inferred from fitting different structural models to a single cluster) of perhaps $\delta(\log\,r_{h})\lesssim 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\propto\rho_{h}^{1/2}$, or $Mr_{h}^{3}={\rm constant}$. Insofar as the half-mass relaxation time scales as $t_{\rm rh}\propto(Mr_{h}^{3})^{1/2}$, and to the extent that $\mu_{\rm ev}\propto M/t_{\rm rh}\propto\rho_{h}^{1/2}$ approximates the average 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 itself a strong hint that relaxation-driven cluster disruption has significantly modified the GCMF at low masses (recall that $Mr_{h}^{3}={\rm constant}$ defines one side of the GC “survival triangle” when the $M$–$\rho_{h}$ plot is recast as $r_{h}$ 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 $r_{h}$ with $M$ is due to significant and real differences in cluster radii, since if $r_{h}$ were intrinsically the same for all GCs, then we would see $M\propto\rho_{h}$ instead.

The middle panel of Figure 1 shows the well-known result that the typical GC mass depends weakly if at all on Galactocentric radius, at least until large $r_{\rm gc}\gtrsim 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 position, the large scatter about it is such that convolving the observed $M$ versus $\rho_{h}$ with the observed $\rho_{h}$ versus $r_{\rm gc}$ results in an almost null dependence of $M$ on $r_{\rm gc}$.

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 $\rho_{h}$ and $r_{\rm gc}$ bins are defined in Table 1, which also gives a few summary 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 $\rho_{h}$ and then as a function of $r_{\rm gc}$. 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 normalization, and is proportional to the GCMF for the whole sample of 146 GCs. (In the middle-left panel of Figure 2, which pertains to clusters distributed tightly around the median $\rho_{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 $\rho_{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 statistical significance of this is very high, and qualitatively it is the behavior expected if $M_{\rm TO}$ owes its existence to cluster disruption at a rate that increases with $\rho_{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 Galactocentric position. In fact, the observed distributions in the two $r_{\rm gc}$ bins inside $\simeq\!10$ kpc are statistically indistinguishable in their entirety, and the main difference at larger $r_{\rm gc}\gtrsim 10$ kpc is a slightly higher proportion of low-mass clusters rather than a large change in $M_{\rm TO}$. All of this is consistent with the primary dependence of the GCMF being that on $\rho_{h}$, since Figure 1 shows that the GC density distribution is not sensitive to Galactocentric position for $r_{\rm gc}\lesssim 10$–20 kpc but has a substantial 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 $\rho_{h}$ in Figure 1 but do not attempt this for the distribution of $\rho_{h}$ over $r_{\rm gc}$—the details of which likely depend on a complicated interplay between the tidal field of the Galaxy, the present and past orbital parameters of clusters, and the structural nonhomology of GCs. To compare our models to the current GCMF as a function of $r_{\rm gc}$, we simply calculate them using the observed $\rho_{h}$ of individual clusters in different ranges of Galactocentric radius.

We assume that the initial GCMF was independent of cluster 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 $\mu_{\rm ev}\propto\rho_{h}^{1/2}$. As we discussed earlier, an approximately time-independent $\mu_{\rm 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 $\rho_{h}$, rather than other plausible measures of cluster density, to estimate $\mu_{\rm ev}$.

Consider first a group of coeval GCs with an initial mass function $dN/d\,\log\,M_{0}$ and a single, time-independent mass-loss rate $\mu_{\rm ev}$. The mass of every cluster decreases linearly as $M(t)=M_{0}-\mu_{\rm ev}t$, and at any later time each has lost the same amount $\Delta\equiv M_{0}-M(t)=\mu_{\rm ev}t$. FZ01 show rigorously that in this case, the evolved and initial GCMFs are related by

[TABLE]

This is the basis for the claim that the mass function scales generically as $dN/d\,\log\,M\propto M^{+1}$ (a $\beta=0$ power law) at low enough $M(t)<\Delta$—that is, for the surviving remnants of clusters with $M_{0}\approx\Delta$—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:

[TABLE]

With $\beta\simeq 2$, this distribution describes the power-law mass functions of young massive clusters in systems like the Antennae galaxies (e.g., Zhang & Fall, 1999). An exponential cut-off at $M_{c}\sim 10^{6}\ M_{\odot}$ 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 density that a single GC with known evaporation rate and age has an instantaneous mass $M$. The time-dependent GCMF of a system of ${\cal N}$ GCs with a range of $\mu_{\rm ev}$ (or ages, or both) is then just the sum of all such individual probability densities:

[TABLE]

Here the total mass losses $\Delta_{i}=(\mu_{\rm ev}t)_{i}$ may differ from cluster to cluster ($t_{i}$ being the age of a single GC) but both $\beta$ and $M_{c}$ are assumed to be constants, independent of $\rho_{h}$ in particular.555Note that $M_{c}$ appears to take on different values in the GCMFs of other galaxies, varying systematically with the total luminosity $L_{\rm gal}$ (Jordán et al., 2007). The reasons for this are unclear, as is the origin of this mass scale in the first place. Given each $\Delta_{i}$, the normalizations $A_{i}$ 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 $\Delta$. They refer to this as an evolved Schechter function and describe its properties in detail (including giving a formula for the turnover mass $M_{\rm TO}$ as a function of $\Delta$ and $M_{c}$) for the case $\beta=2$. Here we note only that, at very young cluster ages or for slow mass-loss rates, such that $\Delta\ll M_{c}$ and only the low-mass, power-law part of the initial GCMF is significantly eroded, any one evolved Schechter function has a peak at $M_{\rm TO}\simeq\Delta/(\beta-1)$ (for $\beta>1$). As $\Delta$ increases relative to $M_{c}$, the turnover at first increases proportionately and the width of the distribution decreases (since the high-mass end at $M\gtrsim M_{\rm TO}$ is largely unchanged). For large $\Delta\gg M_{c}$, however, the peak is bounded above by $M_{\rm TO}\rightarrow M_{c}$ and the width approaches a lower limit.666The increase of $M_{\rm TO}$ and the decrease of the full width of $dN/d\,\log\,M$ for increasing $\Delta$ 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\propto M^{+1}\,\exp(-M/M_{c})$ is a self-similar solution to equation (1). Thus, the dependence of $M_{\rm TO}$ on $\Delta$ is weaker than linear when $M_{c}$ is finite in the initial GCMF of equation (2). Any peak in the full equation (3) for a system of GCs with individual $\Delta$ values is an average of ${\cal N}$ different turnovers and must be calculated numerically.

In their modeling of the Milky Way GC system, FZ01 effectively compute mass functions of the type (3)—based on the same initial conditions and dynamical evolution—with a distribution of $\Delta$ values determined by the orbital parameters of clusters in an idealized, spherical and static logarithmic Galaxy potential (used both to fix $\mu_{\rm ev}$ in terms of cluster 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 $\Delta$. They thus estimate the dynamical mass loss from typical clusters in these systems. Here, we construct models for the Milky Way GCMF using $\Delta$ values given directly by the observed half-mass densities of individual GCs.

We adopt $\beta=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 $\beta=2$ and a single $\Delta\equiv\widehat{\Delta}$ for all surviving globulars. They find $M_{c}\simeq 10^{6}\,M_{\odot}$ and $\widehat{\Delta}=2.3\times 10^{5}\,M_{\odot}$. We use this value of $M_{c}$ in equation (3) and we associate $\widehat{\Delta}$ with the mass loss from clusters at the median half-mass density of the entire GC system, which is $\widehat{\rho}_{h}=246\ M_{\odot}\,{\rm pc}^{-3}$ from the data in Figure 1. Since we are assuming that $\Delta=\mu_{\rm ev}t\propto\rho_{h}^{1/2}t$ for coeval GCs, we therefore stipulate

[TABLE]

for globulars with arbitrary $\rho_{h}$. Assuming a typical GC age of $t=13$ Gyr, this corresponds to a mass-loss rate of

[TABLE]

In §3 we discuss the cluster lifetimes implied by this value of $\mu_{\rm ev}$. We emphasize here that the scaling of $\mu_{\rm ev}$ and $\Delta$ with $\rho_{h}^{1/2}$ follows rather generically from our hypothesis of evaporation-dominated cluster evolution, while the numerical coefficients in equations (4) and (5) are specific to the assumption of $\beta=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 $M_{\rm TO}\simeq 1.6\times 10^{5}\,M_{\odot}$ (magnitude $M_{V}\simeq-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\lesssim\rho_{h}\lesssim 530\ M_{\odot}\,{\rm pc}^{-3}$, and in the two inner radius bins, $r_{\rm gc}\leq 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 $\widehat{\rho}_{h}=246\ M_{\odot}\,{\rm pc}^{-3}$ (see Table 1). Even in the outermost $r_{\rm gc}$ bin, a Kolmogorov-Smirnov (KS) test only marginally rejects the dashed-line model (at the $\simeq\!95\%$ level), because this subsample still includes many GCs at or near the global median $\widehat{\rho}_{h}$ (see Figure 1). By contrast, 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 clusters 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 $\widehat{\Delta}=\Delta(\widehat{\rho}_{h})$.

The solid curves in Figure 2, which are different in every panel, are the superpositions of many different evolved Schechter functions, as in equation (3), with distinct $\Delta$ values given by equation (4) using the observed $\rho_{h}$ of each cluster in the corresponding subsample. These models provide excellent matches to the observed $dN/d\,\log\,M$ in every $\rho_{h}$ and $r_{\rm gc}$ bin, with $\chi^{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 $M_{\rm TO}$ at which each of the solid model GCMFs in Figure 2 peaks. We note that these turnovers increase roughly as $M_{\rm TO}\sim\widehat{\rho}_{h}^{~{}0.3-0.4}$ for our specific binnings in $\rho_{h}$ and $r_{\rm gc}$, somewhat shallower than the $\rho_{h}^{1/2}$ scaling of the cluster mass-loss rate that defines the models. This is partly because of the averaging over individual 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 $\Delta\propto\rho_{h}^{1/2}$, but has a strict upper limit of $M_{\rm TO}\leq M_{c}$.

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 discussion immediately after equation (3), we described how it follows from the increase of $M_{\rm TO}$ with $\Delta\propto\rho_{h}^{1/2}$ for a single evolved Schechter function. In addition, the superposition of many such functions with separate, density-dependent turnovers and widths results in wider GCMFs for cluster subsamples spanning larger ranges of $\rho_{h}$. This accounts in particular for the breadth of the mass function at $r_{\rm gc}\geq 9.4$ kpc. The globulars at these radii have $0.034\leq\rho_{h}\leq 4.1\times 10^{3}\ M_{\odot}\,{\rm pc}^{-3}$, corresponding to individual evolved Schechter functions with turnovers at $2.7\times 10^{3}\lesssim M_{\rm TO}\lesssim 4.0\times 10^{5}\ M_{\odot}$. The composite GCMF in the lower-right panel of Figure 2 is therefore extremely broad and shows a very flat peak, such that an overall $M_{\rm TO}$ 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 consistent with the full-Galaxy average, but a larger dispersion (see also Gnedin 1997).

Finally, if the GCMF evolved dynamically from initial conditions 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 internal 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=\mu_{\rm sh}\propto M/\rho_{h}$, differs significantly from that caused by evaporation alone, $-dM/dt=\mu_{\rm ev}\propto\rho_{h}^{1/2}$. The direct dependence of $\mu_{\rm sh}$ on $M$ ensures that shocks become progressively less important compared to evaporation as clusters lose mass (at a given $\rho_{h}$), and consequently shocks are not likely to have had much effect on the observed GCMF for $M<M_{\rm TO}$. Furthermore, the inverse dependence of $\mu_{\rm sh}$ on $\rho_{h}$ is contrary to the direct dependence of $M_{\rm TO}$ on $\rho_{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 result of long-term cluster disruption according to a mass-loss rule like $\mu_{\rm ev}\propto\rho_{h}^{1/2}$, then it should be possible to reproduce the distribution as a function of any other cluster attribute by using the observed $\rho_{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 $r_{\rm gc}$ binning of §2.2. Here we explore one example in which differences in the GCMFs of two groups of globulars can be seen in this way to follow from differences in their $\rho_{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\geq 0.99$. [Here $c\equiv\log\,(r_{t}/r_{0})$, where $r_{t}$ is the fitted tidal radius and $r_{0}$ 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\propto 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\geq 0.99$. They discuss various options to explain these results, including a suggestion that, if the mass functions of both low- and high-concentration clusters 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 $\rho_{h}$, which either derives 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 $\rho_{h}$ distribution is offset 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 $M_{\rm TO}$, 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\geq 0.99$ (open circles). The curves are again given by equation (3) with $\beta=2$, $M_{c}=10^{6}\ M_{\odot}$, and individual $\Delta$ calculated from the observed cluster $\rho_{h}$ through equation (4). These models peak at $M_{\rm TO}\simeq 4.3\times 10^{4}\ M_{\odot}$ for the $c<0.99$ subsample but at $M_{\rm TO}\simeq 1.8\times 10^{5}\ M_{\odot}$ for $c\geq 0.99$, entirely as a result of the different $\rho_{h}$ involved. The larger width of $dN/d\,\log\,M$ and its shallower slope at any $M\lesssim 10^{5}\ M_{\odot}$ 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\gtrsim 2\times 10^{5}\,M_{\odot}$, above the nominal turnover of the full GCMF (as Smith & Burkert 2002 noted). But this is not surprising, 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 $\simeq\!3$ high-mass clusters with $c<0.99$, where none is found.

The apparent variation of the Milky Way GCMF with internal 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 $\rho_{h}$ and $r_{\rm gc}$ in §2.2. To show this, we have made use of the densities $\rho_{h}$ exactly as observed within the two concentration bins indicated in Figure 3—just as we also took $\rho_{h}$ directly from the data for GCs in different ranges of $r_{\rm gc}$ to construct models for comparison with the observed $dN/d\,\log\,M$ in the right-hand panels of Figure 2. Of course, this is not the same as explaining the distribution of $\rho_{h}$ versus $r_{\rm gc}$ 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 equation (5) above implies cluster lifetimes that compare favorably with those expected from relaxation-driven evaporation. Then we discuss why it is reasonable to approximate $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ in the first place. Finally, we address the issue of possible conflict, in some other models for evaporation-dominated GCMF evolution, between the near-constancy of $M_{\rm TO}$ as a function of $r_{\rm gc}$ 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 $\mu_{\rm ev}$ is just $t_{\rm dis}=M/\mu_{\rm ev}$. It is convenient, for purposes of comparison with evaporation times in the literature, to normalize $t_{\rm dis}$ to the relaxation time of a cluster at its half-mass radius. In general, this is $t_{\rm rh}=0.138M^{1/2}r_{h}^{3/2}/\left[G^{1/2}m_{*}\ln\,(\gamma M/m_{*})\right]$, where $m_{*}$ is the mean stellar mass. For clusters of stars with a single mass, $m_{*}\simeq 0.7\,M_{\odot}$ and $\gamma\simeq 0.4$ are appropriate (Spitzer 1987; Binney & Tremaine 1987, equation [8-72]), in which case equation (5) for $\mu_{\rm ev}$ from our GCMF modeling implies

[TABLE]

Clusters with realistic stellar mass spectra will have slightly different values of $m_{*}$ and a smaller $\gamma$ in the calculation of the relaxation time (Giersz & Heggie, 1996), which changes the numerical value of $t_{\rm dis}/t_{\rm rh}$ somewhat but does not alter any scalings.

We obtained the normalization of $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ in §2.2 by fitting to observed GCMFs constructed by applying a specific mass-to-light ratio $\Upsilon_{V}$ to every cluster, with models assuming a specific form for the initial $dN/d\,\log\,M_{0}$. Thus, the result in equation (6) depends both on the median $\widehat{\Upsilon}_{V}$ and on the power-law index $\beta$ at low masses in the original Schechter-function GCMF. The net scaling, for either single- or multiple-mass clusters, is

[TABLE]

To see the dependence of this dimensionless lifetime on $\widehat{\Upsilon}_{V}$, note that we require $\mu_{\rm ev}\propto\Delta\propto\Upsilon_{V}$ to fit the mass losses of clusters with a given distribution of luminosities (the direct observables), whereas $M/t_{\rm rh}$ is proportional to $\rho_{h}^{1/2}\propto\Upsilon_{V}^{1/2}(L/r_{h}^{3})^{1/2}$. Therefore, $t_{\rm dis}/t_{\rm rh}\propto(M/t_{\rm rh})/\mu_{\rm ev}\propto\Upsilon_{V}^{-1/2}$. The mass-to-light ratios adopted in this paper, with a median value $\widehat{\Upsilon}_{V}\simeq 1.5\ M_{\odot}\,L_{\odot}^{-1}$, are tied directly to dynamical determinations (§2.1).

To understand the dependence on $\beta$ in equation (7), recall first that the coefficients in our expressions for $\Delta$ and $\mu_{\rm ev}$ as functions of $\rho_{h}$ (eqs. [4] and [5]) followed from choosing $\beta=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 function with any $\beta>1$ is $M_{\rm TO}\simeq\Delta/(\beta-1)$ in the limit of low $\Delta\propto\rho_{h}^{1/2}$, and $M_{\rm TO}\rightarrow M_{c}$ for very high $\Delta$. In this sense, the strongest observational constraints on the normalizations of $\Delta$ and $\mu_{\rm ev}$ come from the low-density clusters. All other things being equal, their GCMF can be reproduced with $\beta\neq 2$ if $\Delta$ and $\mu_{\rm ev}$ are multiplied by $(\beta-1)$ at fixed $\rho_{h}$. Therefore, $t_{\rm dis}\propto 1/\mu_{\rm ev}\propto 1/(\beta-1)$. Observations of young massive clusters (e.g., Zhang & Fall, 1999) indicate that $\beta$ is near 2; but if it were slightly shallower, then the cluster lifetimes we infer from the old GCMF would increase accordingly. Even a relatively minor change to $\beta=1.5$ would double $t_{\rm dis}/t_{\rm rh}$ from $\approx\!10$ to $\approx\!20$.

In the model of Hénon (1961) for single-mass clusters evolving self-similarly (fixed ratio $\rho_{t}/\rho_{h}$ of mean densities inside the tidal and half-mass radii) in a steady tidal field ($\rho_{t}$ constant in time), a cluster loses 4.5% of its remaining mass every half-mass relaxation time. The time to complete disruption is therefore $t_{\rm dis}/t_{\rm rh}=1/0.045\simeq 22$. For non-homologous clusters in a steady tidal field, $t_{\rm dis}/t_{\rm rh}$ is a function of central concentration and can differ from the Hénon value by factors of about two. From one-dimensional Fokker-Planck calculations, Gnedin & Ostriker (1997) find $t_{\rm dis}/t_{\rm rh}\simeq 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 theoretical calculations, it is certainly within the range of such calculations. 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 typical lifetimes in these particular theoretical treatments and our estimate of $t_{\rm dis}/t_{\rm rh}$ 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 stellar masses in the clusters may also have contributed to an increase 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 $t_{\rm dis}/t_{\rm rh}$ 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 $t_{\rm dis}$ or $\mu_{\rm 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 $\mu_{\rm ev}\propto\rho_{h}^{1/2}$

3.2.1 Half-mass versus Tidal Density

The dimensionless disruption time in equation (6) is independent of any cluster property other than the Coulomb logarithm because we have used GC half-mass densities to estimate $t_{\rm dis}=M/\mu_{\rm ev}\propto M/\rho_{h}^{1/2}$, while $t_{\rm rh}$ also scales as $M/\rho_{h}^{1/2}$. However, as we mentioned above, the Fokker-Planck calculations of Gnedin & Ostriker (1997) in particular show that $t_{\rm dis}/t_{\rm rh}$ is actually a function of central concentration, $c$, for King (1966) model clusters in steady tidal fields. The constant of proportionality in $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ 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 $t_{\rm dis}/t_{\rm rh}$ on $c$ for single-mass King models, as given by equation (30) of Gnedin & Ostriker (1997). The solid curve is proportional to $(\rho_{h}/\rho_{t})^{1/2}=(r_{t}^{3}/2r_{h}^{3})^{1/2}$, which we have calculated as a function of $c$ for these models and multiplied by a constant to compare directly with $t_{\rm dis}/t_{\rm rh}$. Evidently, there is an approximate equality $t_{\rm dis}/t_{\rm rh}\approx 2.15(\rho_{h}/\rho_{t})^{1/2}$, which holds to within $<\!15\%$ over the range of concentrations shown in Figure 4 (note that all but 6 Galactic GCs have $0.7\leq c\leq 2.5$, corresponding to central potentials $3\lesssim W_{0}\lesssim 11$). Thus, if the evaporation time is written as $t_{\rm dis}\propto t_{\rm rh}(\rho_{h}/\rho_{t})^{1/2}\propto M/\rho_{t}^{1/2}$, then the constant of proportionality in the mass-loss rate $\mu_{\rm ev}\propto M/t_{\rm dis}\propto\rho_{t}^{1/2}$ 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 $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ with only a weak dependence on $c$. An essentially concentration-independent scaling of $\mu_{\rm ev}$ with $\rho_{t}^{1/2}$ 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 assumptions specific to the calculations of either King (1966) or Gnedin & Ostriker (1997).

This suggests that it might have been more natural to specify cluster evaporation rates proportional to $\rho_{t}^{1/2}$ rather than $\rho_{h}^{1/2}$ when developing our GCMF models in §2. For any cluster in a steady tidal field, with a constant $\rho_{t}$, such a choice would also have been automatically consistent with an approximately time-independent $\mu_{\rm ev}$ and the corresponding linear $M(t)$ dependence that we have adopted throughout this paper. As we discussed at the beginning of §2, our decision to work with $\rho_{h}$ rather than $\rho_{t}$ was motivated by the fact that the half-mass density is much better defined in principle and more accurately observed in practice. Nevertheless, re-writing $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ as $\mu_{\rm ev}\propto(\rho_{t}/\rho_{h})^{1/2}\times\rho_{h}^{1/2}$ makes it clear that the validity of our models, with a fixed coefficient in $\mu_{\rm ev}\propto\rho_{h}^{1/2}$, depends on the extent to which variations in $(\rho_{t}/\rho_{h})^{1/2}$ can safely be ignored.

Figure 4 shows that the full range of possible values for $(\rho_{h}/\rho_{t})^{1/2}$ in King-model clusters with $c\geq 0.7$ is only a factor of $\simeq\!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 $\Delta$ and $\mu_{\rm 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 $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ are good approximations to GCMF evolution under a mass-loss law $\mu_{\rm ev}\propto\rho_{t}^{1/2}$, we have repeated the analysis of §2 in full but using the GC tidal densities $\rho_{t}$ (derived from the values of $r_{t}$ listed by Harris 1996) in place of $\rho_{h}$ throughout. All of our main results persist.

For example, the two panels of Figure 5, which are analogous to the left- and rightmost panels of Figure 1 above, show that (1) the GC mass distribution has a clear dependence on $\rho_{t}$, with a lower envelope that is well matched by a line of constant evaporation time, $M\propto\rho_{t}^{1/2}$ (the dashed line in the plot); and (2) although the scatter in the distribution of $\rho_{t}$ over Galactocentric radius is smaller than the scatter in $\rho_{h}$ versus $r_{\rm gc}$, it is still significant. Because the $M$–$r_{\rm gc}$ distribution can now be viewed as the convolution of the $M$–$\rho_{t}$ distribution with the $\rho_{t}$–$r_{\rm gc}$ distribution, the scatter in $\rho_{t}$ versus $r_{\rm gc}$ is again critical in explaining the weak or null dependence of the GCMF on Galactocentric radius. (The $M$–$r_{\rm gc}$ distribution is, of course, unchanged from that shown in the middle panel of Figure 1.)777As was also the case with our earlier plots involving $\rho_{h}$ in Figure 1, the scatter and structure in both panels of Figure 5 are real, since the rms scatter of $\log\,r_{t}$ about the best-fit lines to either of $\log\,M$ or $\log\,r_{\rm gc}$ is 0.3–0.35 while the rms errorbars based on formal fitting uncertainties are in the range $\delta(\log\,r_{t})\simeq 0.05$–0.15 for a variety of models (McLaughlin & van der Marel, 2005).

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 models for these distributions are based as before on equation (3) with $\beta=2$, but now the total mass lost from any GC is estimated from its tidal density rather than its half-mass density. Specifically, we take

[TABLE]

The numerical coefficient in equation (8) is such that it gives a $\Delta$ identical to that in equation (4) for a GC with $\rho_{h}/\rho_{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 $\rho_{t}$, just as it does with increasing $\rho_{h}$. Meanwhile, the solid curves are subsample-specific model GCMFs, obtained by using the observed tidal density of each cluster in any $\rho_{t}$ or $r_{\rm gc}$ bin to specify individual $\Delta$ values via equation (8) for each of the evolved Schechter functions in the summation of equation (3). As expected, there is no appreciable difference, in terms of the fits to any of the observed GCMFs, between these models based on evaporation rates $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ and our original models with $\mu_{\rm ev}\propto\rho_{h}^{1/2}$.

3.2.2 Retarded Evaporation

Another potential concern comes from recent arguments (see especially Baumgardt 2001; Baumgardt & Makino 2003) that the total evaporation time of a tidally limited cluster is not simply a multiple of an internal two-body relaxation time, $t_{\rm rlx}\propto(Mr^{3})^{1/2}$, but depends on both $t_{\rm rlx}$ and the crossing time $t_{\rm cr}\propto(M/r^{3})^{-1/2}$ through the combination $t_{\rm dis}\propto t_{\rm rlx}^{x}t_{\rm cr}^{1-x}$, with $x<1$. The mass-loss rate $\mu_{\rm ev}\propto M/t_{\rm dis}$ then scales as $M^{3/2-x}r^{-3/2}$, which for $x\neq 1$ differs from the rates $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ and $\mu_{\rm ev}\propto\rho_{t}^{1/2}$ that we have so far adopted. However, our GCMF models are still meaningful, because postulating $t_{\rm dis}\propto t_{\rm rlx}^{x}t_{\rm cr}^{1-x}$ implies a dependence of $\mu_{\rm ev}$ on a measure of cluster density that is, once again, well approximated by $\rho_{h}^{1/2}$ for Galactic GCs. Before showing this, we briefly discuss the reasons and the evidence for a possible dependence of $t_{\rm dis}$ on both $t_{\rm rlx}$ and $t_{\rm cr}$.

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 $t_{\rm dis}\propto t_{\rm rlx}$ is expected (and confirmed by $N$-body simulations; e.g., Baumgardt 2001). However, 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, although the stars are still scattered to near- and above-escape energies on the timescale $t_{\rm rlx}$, they require some additional time to actually leave the cluster. This escape timescale is related fundamentally to $t_{\rm cr}$ (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 (“retardation”) of the overall evaporation rate (Chandrasekhar, 1942; King, 1959; Takahashi & Portegies Zwart, 1998, 2000; Fukushige & Heggie, 2000; Baumgardt, 2001) and a lengthening of the cluster lifetime $t_{\rm dis}$, by a factor that can be expected to increase with the ratio $t_{\rm cr}/t_{\rm rlx}$. If this factor scales as $(t_{\rm cr}/t_{\rm rlx})^{1-x}$ for some $x<1$, then $t_{\rm dis}\propto t_{\rm rlx}\,(t_{\rm cr}/t_{\rm rlx})^{1-x}=t_{\rm rlx}^{x}t_{\rm cr}^{1-x}$.

While such a retardation of evaporation can be expected to occur at some level in all clusters, there are physical subtleties in the effect that are probably not captured adequately by a simple re-parametrization of lifetimes as $t_{\rm dis}\propto t_{\rm rlx}^{x}t_{\rm cr}^{1-x}$. In particular, it is unlikely that this expression can hold for clusters of all masses with a single value of $x<1$. Since $t_{\rm cr}/t_{\rm rlx}\propto M^{-1}$, very massive clusters have $t_{\rm cr}\ll t_{\rm rlx}$, and stars scattered to greater than escape energies by relaxation cross the tidal boundary effectively instantaneously—implying that the standard $t_{\rm dis}\propto t_{\rm rlx}$, or $x\rightarrow 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 $t_{\rm dis}<t_{\rm rlx}$ would obtain at high enough masses; see Baumgardt (2001) for further discussion. Unfortunately, “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 circular or moderately eccentric orbits in galactic potentials that are static and spherical. This means that any tidal perturbations 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 realistic situations, the galactic potential would be time-dependent and non-spherical and there might be additional tidal perturbations, 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 effect to some degree. Further study is therefore needed to determine the regime of validity of the formula $t_{\rm dis}\propto t_{\rm rlx}^{x}t_{\rm cr}^{1-x}$ and its possible modification outside this regime.

In the meantime, Baumgardt (2001) and Baumgardt & Makino (2003; hereafter BM03) have fitted this formula to the lifetimes of a suite of $N$-body clusters with initial masses $M_{0}\lesssim 7\times 10^{4}\,M_{\odot}$ and several different initial concentrations and orbital eccentricities. BM03 at first write $t_{\rm dis}$ in terms of the relaxation and crossing times of clusters at their half-mass radii, so that $t_{\rm rlx}\propto(Mr_{h}^{3})^{1/2}$, $t_{\rm cr}\propto(M/r_{h}^{3})^{-1/2}$, and $t_{\rm dis}\propto M^{x-1/2}r_{h}^{3/2}$ (see their equation [5]). However, they immediately take a factor of $(r_{t}/r_{h})^{3/2}$ out from the normalization of this scaling—in effect to obtain $t_{\rm dis}\propto M^{x-1/2}r_{t}^{3/2}$ with a different constant of proportionality—and then use a simple definition of the tidal radius (their equation [1], $r_{t}^{3}=GMr_{p}^{2}/2V_{c}^{2}$, which is appropriate for a circular orbit of radius $r_{p}$ in a logarithmic potential with circular speed $V_{c}$; see Innanen, Harris, & Webbink 1983) to obtain the total lifetime of a cluster as a function of its initial mass, perigalactic distance, and $V_{c}$ (their equation [7]). A single exponent $x\simeq 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 $t_{\rm rlx}$ and $t_{\rm cr}$ were fixed at $r_{h}$ rather than $r_{t}$, then $t_{\rm dis}$ would have an additional concentration dependence, related to the ratio $(r_{t}/r_{h})^{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 $t_{\rm dis}\propto M^{x-1/2}r_{t}^{3/2}$ for general $x\leq 1$; but we do not substitute a potential- and orbit-specific formula for $r_{t}$ in terms of $r_{p}$ and galactic properties such as $V_{c}$. Instead, to keep the emphasis entirely on cluster densities, we re-write the scaling of the lifetime in terms of the mean surface density inside the tidal radius, $\Sigma_{t}\equiv M/\pi r_{t}^{2}$, and the corresponding volume density $\rho_{t}=3M/4\pi r_{t}^{3}$. This leads to $t_{\rm dis}\propto M\,\Sigma_{t}^{-3(1-x)}\rho_{t}^{-2(x-3/4)}$, which then implies

[TABLE]

Clearly, the standard $\mu_{\rm ev}\propto\rho_{t}^{1/2}$, which we have already discussed, is recovered for $x=1$; while for $x=0.75$, we have the equally straightforward $\mu_{\rm ev}\propto\Sigma_{t}^{3/4}$.

BM03 find that, even with the retarded evaporation implied by $x\simeq 0.75$, the masses of their simulated clusters still decrease approximately linearly with time after stellar-evolution effects (which are only important for the first few $10^{8}$ 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 evaporation we might expect the current mass function to depend fundamentally on $\Sigma_{t}$ rather than $\rho_{h}$ or $\rho_{t}$. But because $M(t)$ still decreases nearly linearly with $t$, only now with $\mu_{\rm ev}\propto\Sigma_{t}^{3/4}$ for each cluster, the shape of the evolved GCMF and its dependence on $\Sigma_{t}$ should resemble our earlier results for $\rho_{h}$ and $\rho_{t}$.

We have confirmed this expectation by repeating all of our analyses in §2 again, now using $\mu_{\rm ev}\propto\Sigma_{t}^{3/4}$ to estimate cluster mass-loss rates. As before, we calculate $\Sigma_{t}$ from the data in the Harris (1996) catalogue, although we caution once more that the tidal radii, and thus the derived $\Sigma_{t}$, are more uncertain than $r_{h}$ and $\rho_{h}$.

Figure 7, which should be compared to Figures 1 and 5 above, shows that the average Galactic GC mass increases systematically with $\Sigma_{t}$; that the lower envelope of the $M$–$\Sigma_{t}$ distribution is described well by $M\propto\Sigma_{t}^{3/4}$ (the dashed line in the left-hand panel of Figure 7), which is a locus of constant lifetime against evaporation for $\mu_{\rm ev}\propto\Sigma_{t}^{3/4}$; and that the scatter in the distribution of cluster $\Sigma_{t}$ versus Galactocentric radius (right-hand panel of the figure) is substantial, as required to account for the almost non-existent correlation between $M$ and $r_{\rm gc}$.

The left-hand side of Figure 8 shows the mass functions of globulars in three bins of $\Sigma_{t}$, as defined in each panel. The right-hand side of the figure shows $dN/d\,\log\,M$ in the same three intervals of $r_{\rm gc}$ 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 system. Once again, compared to the average $M_{\rm TO}$, the observed turnover mass is significantly lower for clusters in the lowest $\Sigma_{t}$ bin and higher for clusters in the highest $\Sigma_{t}$ bin, while the width of $dN/d\,\log\,M$ decreases noticeably as $\Sigma_{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 $\beta=2$ assumed but with total mass losses estimated individually for each GC in any $\Sigma_{t}$ or $r_{\rm gc}$ bin according to $\Delta\propto\Sigma_{t}^{3/4}$ rather than $\Delta\propto\rho_{h}^{1/2}$ or $\Delta\propto\rho_{t}^{1/2}$. However, it turns out not to be necessary to change the normalization of $\Delta\propto\rho_{h}^{1/2}$ in equation (4) to achieve good fits to the observed GCMF as a function of either $\Sigma_{t}$ or $r_{\rm gc}$. Thus, in Figure 8 we have simply used

[TABLE]

The fits of these models, based on $t_{\rm dis}\propto t_{\rm rlx}^{x}t_{\rm cr}^{1-x}$ with $x\simeq 0.75$, are indistinguishable from the fits of our original models based on the standard $t_{\rm dis}\propto t_{\rm rlx}$, i.e., $x=1$. (We have confirmed that adopting individual $\Delta$ 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 numerical coefficient, but we note that this follows empirically from the fact that the measured $\rho_{h}$ and $\Sigma_{t}$ of Galactic GCs are consistent with the simple near-equality, $\rho_{h}/M_{\odot}\,{\rm pc}^{-3}\approx(\Sigma_{t}/M_{\odot}\,{\rm pc}^{-2})^{1.5}$ in the mean. This is illustrated in Figure 9, which also shows that there is significant scatter about the relation.888Although it may be only a coincidence that the constant of proportionality in $\rho_{h}\propto\Sigma_{t}^{1.5}$ is so near unity, the basic scaling itself holds because combining the observed correlation between cluster mass and central concentration (Djorgovski & Meylan, 1994; McLaughlin, 2000) with the intrinsic dependence of $r_{t}/r_{h}$ on $c$ in King models leads roughly to $(r_{t}/r_{h})\propto M^{1/6}$. However, this scatter does not correlate with cluster mass or Galactocentric radius. From a pragmatic point of view, therefore, $\rho_{h}^{1/2}$ and $\Sigma_{t}^{3/4}$ are near enough to interchangeable for our purposes, and there is no practical difference between 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

[TABLE]

to a formula implicit in BM03. Starting with their equation (7) for the lifetime $t_{\rm dis}$ as a function of initial cluster mass and perigalactic distance and circular speed in a logarithmic halo potential; using their $x=0.75$ and their normalization of $1.91\times 10^{6}$ yr, multiplied as in their equation (9) by $(1+e)$ to allow for eccentric orbits with apo- and perigalactic distances related by $e\equiv(r_{a}-r_{p})/(r_{a}+r_{p})$; inserting their equation (1) for $r_{t}$; taking the mean mass of cluster stars to be $m_{*}=0.55\,M_{\odot}$, as they do; using $\gamma=0.02$ as they do in the Coulomb logarithm, $\ln(\gamma M_{0}/m_{*})$; and defining $\Sigma_{t,0}\equiv M_{0}/\pi r_{t,0}^{2}$ (the subscript 0 denoting initial values), we obtain

[TABLE]

This is appropriate for clusters that just fill their Roche lobes at perigalacticon, which is where $\Sigma_{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 $\mu_{\rm ev}$ is a factor of $\approx\!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 $\mu_{\rm 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\simeq 0.5$–0.6 is typical for tracers with an isotropic velocity distribution in a logarithmic 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 $\mu_{\rm ev}\propto\rho_{h}^{1/2}$. Moreover, our new estimate of $\mu_{\rm ev}$ and that in BM03 are still subject to their own, separate uncertainties and reflect different idealizations and assumptions. For example, our rate still depends on the exact power-law exponent $\beta$ at low masses in the initial GCMF, as discussed after equation (7); while the rate from BM03 still neglects gravitational shocks from disk crossings and passages by a discrete galactic bulge, and may additionally be biased low for $M_{0}>10^{5}\,M_{\odot}$ if $x>0.75$ at such masses. All of this—not to mention again the large uncertainties and possible systematics in the estimates of tidal radii needed to calculate $\Sigma_{t}$—makes the near agreement between equations (11) and (12) more striking than any apparent discrepancy.

In summary, although the relation $\mu_{\rm ev}\propto\rho_{h}^{1/2}\simeq{\rm constant}$ in time is rigorously correct only in rather specific circumstances, our GCMF models based on it in §2 are good proxies, in all respects, for models based on other plausible characterizations 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 to adopt procedures based on $\rho_{h}$ rather than $\rho_{t}$ or $\Sigma_{t}$ because of the difficulty or impossibility of estimating tidal radii.

3.3. $M_{\rm TO}$ versus $r_{\rm gc}$, 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 $\rho_{h}$ (or $\rho_{t}$, or $\Sigma_{t}$) of clusters in relatively narrow ranges of Galactocentric position to show that such models are consistent with the current near-constancy of the GCMF as a function of $r_{\rm gc}$. 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 $r_{\rm gc}$ 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 galaxies of GCs. Taking these to be spherical and static for a Hubble time allows the use of standard tidal-limitation formulae to write GC densities ab initio in terms of the (fixed) pericenters $r_{p}$ of unique orbits in the adopted potentials. Cluster relaxation times and mass-loss rates $\mu_{\rm ev}$ then follow as functions of $r_{p}$ 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 $r_{p}$; how many clusters with a given $r_{p}$ survive; and what the distributions of $r_{p}$ and all dependent cluster properties are at any instantaneous position $r_{\rm gc}$.

In this approach, if the GCMF began with a power-law rise towards low masses and its current peak is due entirely to cluster disruption, then a dependence of $M_{\rm TO}$ on $r_{p}$ is expected in general, because the densities of tidally limited GCs decrease with increasing $r_{p}$. Thus, models along these lines that assume 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 ratio $r_{\rm gc}/r_{p}$ is independent of position) have typically had difficulty in accounting for the observed weak or non-correlation between $M_{\rm TO}$ and present $r_{\rm gc}$ in large galaxies. This is particularly 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 potentially less of a concern if $dN/d\,\log\,M$ started as a Schechter function with an exponential cut-off at masses $M>M_{c}$, as we have assumed, since then the existence of a strict upper bound $M_{\rm TO}\leq M_{c}$ (§2.2) means that the dependence of an evaporation-evolved $M_{\rm TO}$ on $r_{p}$ and $r_{\rm gc}$ must saturate for small enough galactocentric radii (high enough GC densities). Even so, the “scale-free” models of FZ01, in which $M_{c}\simeq 10^{6}\,M_{\odot}$ and all GCs in a Milky Way-like galaxy potential have the same time-averaged $r_{\rm gc}/r_{p}$, still predict a gradient in $M_{\rm TO}$ versus $r_{\rm gc}$ 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 $r_{\rm gc}$ 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 orbit increases with galactocentric distance (the time average of $r_{\rm gc}/r_{p}$ increases with radius), such that globulars spread over a larger range of current $r_{\rm gc}$ can have more similar $r_{p}$ and associated $M_{\rm TO}$. 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).999The fact that clusters on radial orbits are preferentially disrupted lessens any inconsistency between the radial anisotropy required in the initial velocity distribution and observational constraints on the present velocity distribution. Subsequently, Vesperini et al. (2003) constructed broadly similar models for the GCMF of the Virgo elliptical M87 and concluded that there, too, a variable radial velocity anisotropy is required to match the observed $M_{\rm TO}$ versus $r_{\rm gc}$; 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 $r_{\rm gc}$ (Romanowsky & Kochanek, 2001; Côté et al., 2001).

These results certainly suggest that some element is lacking in $r_{\rm gc}$-oriented GCMF models developed as outlined above. But they do not mean that the fault lies with the main hypothesis, 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 possible that reasonable changes to one or more of these ancillary assumptions could make the models compatible with the observed 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 successful.

One issue is that previous models have always specified evaporation rates a priori as functions of cluster density (or orbital pericenter), usually normalizing $\mu_{\rm ev}$ so that $t_{\rm dis}/t_{\rm rh}\simeq 20$–40 as in standard treatments of two-body relaxation. However, following our discussion in §3.1 and §3.2, it would seem worthwhile to investigate these models with $\mu_{\rm ev}$ increased at fixed $\rho_{h}$ or $r_{p}$ to allow $t_{\rm dis}/t_{\rm rh}\approx 10$ (if $\beta\simeq 2$ for the low-mass power-law part of the initial GCMF).

FZ01 and Vesperini et al. (2003) both consider velocity distributions parametrized by a galactocentric anisotropy radius, $R_{A}$, inside of which a cluster system is essentially isotropic and beyond which it is increasingly dominated by radial orbits. In these terms, the difficulty with the published models is that, to reproduce the observed insensitivity of $M_{\rm TO}$ to $r_{\rm gc}$ given standard normalizations of $\mu_{\rm ev}$, they require values of $R_{A}$ that are smaller than allowed by observations (especially for M87). Increasing $R_{A}$ to more realistic values while keeping the normalization of $\mu_{\rm ev}$ fixed leads to a stronger gradient in $M_{\rm TO}$: the orbits of GCs at small $r_{\rm gc}\lesssim R_{A}$ remain closely isotropic and the typical $r_{p}$ and $M_{\rm TO}$ are essentially unchanged, while at large galactocentric distances the cluster orbits are on average less radial than before, with larger $r_{p}$, lower densities, and lower evolved $M_{\rm TO}$ for a given $r_{\rm gc}$. This effect is illustrated, for example, in Figure 9 of FZ01. However, it can be compensated at least in part by increasing $\mu_{\rm ev}$ by a common factor for all GCs, with the new, larger $R_{A}$ fixed, if the initial mass function is assumed to have been a Schechter function rather than a pure power law extending 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 densities, i.e., those at large $r_{\rm gc}$ and $r_{p}$; but the increase in $M_{\rm TO}$ will be smaller, and eventually even negligible, for higher-density clusters at progressively smaller $r_{\rm gc}$—again because $M_{\rm TO}$ grows less than linearly with $\mu_{\rm ev}\propto\rho_{h}^{1/2}$ when there is an upper limit $M_{\rm TO}<M_{c}$ due to an exponential cut-off in the initial $dN/d\,\log\,M_{0}$. Thus, the qualitative effect of increasing the normalization of $\mu_{\rm ev}$ in models with radially varying GC velocity anisotropy is to weaken the amount of radial-orbit bias required to fit an observed $M_{\rm TO}$ versus $r_{\rm gc}$.

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 $\mu_{\rm ev}$ are not changed, much of the burden for the weakening or erasing of any initial gradients in $M_{\rm TO}$ versus $r_{\rm gc}$ 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–$M_{\rm TO}$ relations. Any position dependences in the GC $\rho_{h}$ distribution and in $M_{\rm TO}$ itself for the final galaxy are therefore bound to be weaker, more scattered, and more difficult to relate 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 evaporation 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 system of GCs in a single, static galaxy potential. The mass-evolution curve for each cluster is approximately a straight line, $M(t)\simeq M_{0}-\mu_{\rm ev}t$, with $\mu_{\rm ev}$ depending on some measure of internal density, which may be $\rho_{h}^{1/2}$, $\rho_{t}^{1/2}$, or $\Sigma_{t}^{3/4}$. The average mass-evolution curve for the entire system of clusters is also approximately linear, $\langle M(t)\rangle\simeq\langle M_{0}\rangle-\langle\mu_{\rm ev}\rangle 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 $\langle\mu_{\rm ev}\rangle\propto\langle\rho_{h}^{1/2}\rangle$ (say) will differ even less between the pre- and post-merger systems. Thus, although using instantaneous densities to estimate the past $\mu_{\rm ev}$ of individual clusters may err on the high side for some clusters and on the low side for others, these errors will average away to a small or even zero net bias. The approximation $\mu_{\rm ev}\simeq{\rm constant}$ in time in our GCMF models will then still be valid in the mean, and the average $\langle M(t)\rangle$ dependence of sufficiently 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 $M_{\rm gal}\propto V_{c}^{3}$ or $M_{\rm gal}\propto V_{c}^{4}$. By the virial theorem, the average densities scale as $\rho_{\rm gal}\propto V_{c}^{6}/M^{2}$, and thus $\rho_{\rm gal}\propto M_{\rm gal}^{0}$ or $\rho_{\rm gal}\propto M_{\rm gal}^{-1/2}$. Insofar as $\langle\rho_{h}\rangle\propto\langle\rho_{t}\rangle\propto\rho_{\rm gal}$ for the GCs, the system-wide average $\langle\mu_{\rm ev}\rangle\propto\langle\rho_{h}^{1/2}\rangle$ should therefore not change drastically even after a major merger between two fundamental-plane galaxies; at most, the ratio of final to initial $\langle\mu_{\rm ev}\rangle$ will be roughly of order the $-1/4$ power of the ratio of final to initial $M_{\rm gal}$. 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 globular clusters in the Milky Way depends significantly on cluster half-mass density, $\rho_{h}$, with the peak or turnover mass $M_{\rm TO}$ increasing and the width of the distribution decreasing as $\rho_{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 Antennae 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 $M_{\rm TO}$ increases with cluster density favors evaporation over external gravitational shocks as the primary mechanism of low-mass cluster disruption, since the mass-loss rates associated with shocks depend inversely on cluster density and directly on cluster mass. Our results therefore add to previous arguments supporting an interpretation of the GCMF in terms of evaporation-dominated evolution, based on the fact that $dN/d\,\log\,M$ scales as $M^{1-\beta}$ with $\beta\simeq 0$ in the low-mass limit (Fall & Zhang, 2001).

The observed GCMF as a function of $\rho_{h}$ is fitted well by simple models in which the initial distribution was a Schechter function, $dN/d\,\log\,M_{0}\propto M_{0}^{1-\beta}\,\exp\,\left(-M_{0}/M_{c}\right)$ with $\beta=2$ and $M_{c}\simeq 10^{6}\,M_{\odot}$ 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 densities as $\mu_{\rm ev}\propto\rho_{h}^{1/2}$. We have shown that, although this prescription is approximate, it captures the main physical dependence of relaxation-driven evaporation. In particular, it leads to model GCMFs that are entirely consistent with those resulting from alternative characterizations of evaporation rates in terms of cluster tidal densities $\rho_{t}$ or mean surface densities $\Sigma_{t}$ (§3.2). The normalization of $\mu_{\rm ev}$ at a given $\rho_{h}$ (or $\rho_{t}$, or $\Sigma_{t}$) required to fit the GCMF implies total cluster lifetimes that are within range of the lifetimes typically obtained in theoretical 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 $r_{\rm gc}$ and using their (individual) observed densities as direct input to our models yields dynamically evolved GCMFs as functions of $c$ and $r_{\rm gc}$ that agree well with all data. This again indicates that the most fundamental physical dependence in the GCMF is that on cluster density. Moreover, our models for $dN/d\,\log\,M$ versus $r_{\rm gc}$ obtained in this way are consistent in particular with the well-known insensitivity of the GCMF peak mass to Galactocentric position. This is seen to follow from a significant variation of $M_{\rm TO}$ with $\rho_{h}$ (or $\rho_{t}$, or $\Sigma_{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 distribution to explain the observed weak variation of $M_{\rm TO}$ with $r_{\rm gc}$; indeed, we have made no predictions or assumptions whatsoever 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–$r_{\rm gc}$ relations for GCs—which we have not attempted to do here. The apparent need in some current models for a strong bias towards high-eccentricity cluster orbits to explain the near-constancy of $M_{\rm TO}$ versus $r_{\rm gc}$ 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 disruption.

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 require the estimation of approximate mass-loss rates using cluster half-mass densities rather than tidal quantities, simply 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 observe. Chandar, Fall, & McLaughlin (2007) have recently shown that the peak mass of the GCMF in the Sombrero galaxy (M104) increases with $\rho_{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 relatively straightforward to pursue similar studies in other nearby galaxies.

We thank Michele Trenti, Douglas Heggie, Bill Harris, Rupali Chandar, and Bruce Elmegreen for helpful discussions and comments. SMF acknowledges support from the Ambrose 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.