The evolution of the stellar mass function in star clusters^{†}^{†}thanks: The models presented in this paper will be publicly available in electronic form at the CDS via anonymous ftp to http://cdsweb.ustrasbg.fr/ (130.79.128.5) and also at http://www.astro.uu.nl/~kruijs.
Key Words.:
stellar dynamics – stars: kinematics – (Galaxy:) globular clusters: general – (Galaxy:) open clusters and associations: general – galaxies: star clusters – galaxies: stellar contentAbstract
Context:The dynamical ejection of stars from star clusters affects the shape of the stellar mass function (MF) in these clusters, because the escape probability of a star depends on its mass. This is found in body simulations and has been approximated in analytical cluster models by fitting the evolution of the MF. Both approaches are naturally restricted to the set of boundary conditions for which the simulations were performed.
Aims:The objective of this paper is to provide and to apply a simple physical model for the evolution of the MF in star clusters for a large range of the parameter space. It should also offer a new perspective on the results from body simulations.
Methods:A simple, physically selfcontained model for the evolution of the stellar MF in star clusters is derived from the basic principles of twobody encounters and energy considerations. It is independent of the adopted mass loss rate or initial mass function (IMF), and contains stellar evolution, stellar remnant retention, dynamical dissolution in a tidal field, and mass segregation.
Results:The MF evolution in star clusters depends on the disruption time, remnant retention fraction, initialfinal stellar mass relation, and IMF. Lowmass stars are preferentially ejected after Myr. Before that time, masses around 15—20% of the maximum stellar mass are lost due to their rapid twobody relaxation with the massive stars that still exist at young ages. The degree of lowmass star depletion grows for increasing disruption times, but can be quenched when a large fraction of massive remnants is retained. The highly depleted MFs of certain Galactic globular clusters are explained by the enhanced lowmass star depletion that occurs for low remnant retention fractions. Unless the retention fraction is exceptionally large, dynamical evolution always decreases the masstolight ratio. The retention of black holes reduces the fraction of the cluster mass in remnants because white dwarfs and neutron stars have masses that are efficiently ejected by black holes.
Conclusions:The modeled evolution of the MF is consistent with body simulations when adopting identical boundary conditions. However, it is found that the results from body simulations only hold for their specific boundary conditions and should not be generalised to all clusters. It is concluded that the model provides an efficient method do understand the evolution of the stellar MF in star clusters under widely varying conditions.
1 Introduction
The evaporation of star clusters is known to change the shape of the underlying stellar mass function^{1}^{1}1Hereafter, ‘mass function’ is referred to as ‘MF’. (Hénon 1969; Chernoff & Weinberg 1990; Vesperini & Heggie 1997; Takahashi & Portegies Zwart 2000; Portegies Zwart et al. 2001; Baumgardt & Makino 2003). This phenomenon has been used to explain the observed MFs in globular clusters (Richer et al. 1991; De Marchi et al. 2007; De Marchi & Pulone 2007), which are flatter than typical initial mass functions (IMFs, e.g. Salpeter 1955; Kroupa 2001). In addition, the effect of a changing MF on cluster photometry has been investigated (Lamers et al. 2006; Kruijssen & Lamers 2008; Anders et al. 2009). This has been shown to explain the low masstolight ratios of globular clusters (Kruijssen 2008; Kruijssen & Mieske 2009) and to have a pronounced effect on the inferred globular cluster mass function (Kruijssen & Portegies Zwart 2009).
The existing parameterised cluster models that incorporate a description of lowmass star depletion are restricted by the physically selfcontained models on which they are based. Some studies (Lamers et al. 2006; Kruijssen & Lamers 2008) assume an increasing lower stellar mass limit to account for the evolving MF, others (Anders et al. 2009) fit a changing MF slope to body simulations. In both cases, the models are accurate for a certain range of boundary conditions, but they do not include a physical model and are therefore lacking flexibility. While body simulations do include the appropriate physics, they are very timeconsuming. As a result, only a limited number of clusters can be simulated and the applicability of the simulations is thus restricted to the specific set of boundary conditions for which they have been run.
It would be desirable to obtain a simple physical model for the evolution of the MF, which would have a short runtime and could be used independently of body simulations. Forty years ago, a pioneering first approach to such a model was made by Hénon (1969), who considered the stellar massdependent escape rate of stars from star clusters. However, the applicability of his model was limited due to a number of assumptions that influenced the results. First of all, Hénon (1969) assumed that the clusters exist in isolation and neglected the tidal field. As a consequence, the ejection of a star could only occur by a single, close encounter and the repeated effect of twobody relaxation was not included. Secondly, the distribution of stars was independent of stellar mass, i.e. mass segregation was not included. Both mass segregation and the influence of a tidal field are observed in real clusters, and can be expected to affect the evolution of the MF.
The aim of this paper is to derive a physical description of the evolution of the stellar MF in star clusters, alleviating the assumptions that were made by Hénon (1969). This should explain the results found in body simulations and observations, while providing the required flexibility to explore the properties of star clusters with simple, physically selfcontained models. The outline of this paper is as follows. In Sect. 2, total mass evolution of star clusters is discussed. A recipe for the evolution of the MF is derived in Sect. 3, covering stellar evolution, the retain of stellar remnants, dynamical dissolution and mass segregation. The model is compared to body simulations in Sect. 4. In Sect. 5, the model is applied to assess the evolution of the MF for different disruption times and remnant retention fractions. The consequences for other cluster properties are also considered. This paper is concluded with a discussion of the results and their implications.
2 The mass evolution of star clusters
The mass of star clusters decreases due to stellar evolution and dynamical dissolution. This is expressed mathematically as
(1) 
with the cluster mass, and the subscripts ‘ev’ and ‘dis’ denoting stellar evolution and dynamical dissolution. The contribution of stellar evolution to the mass loss is derived from the decrease of the maximum stellar mass with time and depends on the adopted stellar evolution model.
The dynamical evaporation of star clusters is increasingly well understood. Over the past years it has become clear that clusters lose mass on a disruption timescale that is proportional to a combination of the halfmass relaxation time and the crossing time as (e.g. Baumgardt 2001; Baumgardt & Makino 2003; Gieles & Baumgardt 2008). It is found that —0.80, depending on the concentration () or King parameter () of the cluster (Baumgardt & Makino 2003). This proportionality leads to a disruption timescale that scales with the present day mass as (Lamers et al. 2005):
(2) 
with the cluster mass, the dissolution timescale parameter which sets the rapidity of dissolution and depends on the cluster environment, and a constant related to . Lamers et al. (2009) find for and for . This timescale implies a mass loss rate due to dissolution that can be described with the simple relation
(3) 
which can be integrated for the mass evolution of the cluster due to dynamical dissolution.
The above formulation of the cluster mass evolution was extended to include stellar remnants, photometric cluster evolution, and a simple description of the MF in the SPACE cluster models (Kruijssen & Lamers 2008). Stellar remnants were accounted for by assuming initialfinal mass relations (similar to Sect. 3.1 of the present paper), while the photometric evolution was computed by integrating stellar isochrones from the Padova group (Bertelli et al. 1994; Girardi et al. 2000). The description of lowmass star depletion followed the simple model from Lamers et al. (2006) in which the minimum stellar mass of the MF increases with time.
The present study provides a new description of the evolution of the MF which is based on fundamental principles, and does not depend on the above prescription for the total mass evolution. In addition, the latest Padova models (Marigo et al. 2008) are incorporated to calculate the photometric cluster evolution. These improvements update the SPACE cluster models.
3 The evolution of the stellar mass function
To describe the evolution of the MF, the effects of stellar evolution, stellar remnant production, and dynamical dissolution need to be included. While the focus of this paper lies with the effects of dissolution, a proper treatment of stellar evolution is essential. This is described first, before presenting a model for cluster dissolution.^{2}^{2}2The model presented in this paper is independent of the mass loss rate and of the form of the IMF , but for explanatory purposes a Kroupa (2001) IMF is adopted later on.
3.1 Stellar evolution
The influence of stellar evolution on the MF is twofold. First of all, the maximum stellar mass decreases, because at any time during cluster evolution the most massive stars reach the end of their lives. Secondly, the stellar remnants that are created upon the death of these massive stars constitute a part of the MF that can only be lost from the cluster by dynamical mechanisms.
The maximum stellar mass in the cluster as a function of its age is taken from the Padova 2008 isochrones (Marigo et al. 2008) for metallicities in the range —0.03. The stellar remnant masses are computed from their progenitor stellar mass using initialfinal mass relations. Following Kruijssen & Lamers (2008), for white dwarfs () the relation from Kalirai et al. (2008) is adopted:
(4) 
which holds for all ages that are covered by the Padova isochrones. For neutron stars () the relation from Nomoto et al. (1988) is used:
(5) 
while for black holes () a simple relation is assumed that is in acceptable agreement with theoretically predicted masses of stellar mass black holes (Fryer & Kalogera 2001):
(6) 
With these relations, the remnant MF is computed from conservation of numbers as
(7) 
with denoting the appropriate remnant type, representing its MF, denoting the cluster massdependent fraction of these remnants that is retained after applying kick velocities, and representing the progenitor MF.
For a given velocity dispersion of remnants, the retention fraction of each remnant type depends on the local escape velocity , which is related to the potential as . Stellar remnants are predominantly produced in the cluster centre in the case of mass segregation, which is reached most rapidly for massive stars (see Sect. 3.2). For a Plummer (1911, also see Eq. 9) potential this implies that upon remnant production , with the gravitational constant and the Plummer radius. Adopting a Maxwellian distribution of velocities that is truncated at , it is straightforward to show that
(8) 
where is a normalisation constant to account for the truncation of the velocity distribution and , with denoting the total velocity dispersion of the produced remnant type, which arises from the central velocity dispersion in the cluster (e.g. Heggie & Hut 2003) and the velocity dispersion of the exerted kick . The normalisation constant then follows as .
Typical values of the kick velocity dispersion are given in literature. White dwarf kicks have recently been proposed to be of order km s (Davis et al. 2008; Fregeau et al. 2009). For neutron stars km s is adopted, which is a somewhat conservative estimate with respect to theory, but it agrees reasonably well with observed neutron star numbers in globular clusters and represents a compromise between single star and binary channels (for estimates of the retention fraction and discussions of the ‘neutron star retention problem’ see Lyne & Lorimer 1994; Drukier 1996; Arzoumanian et al. 2002; Pfahl et al. 2002). Gravitational wave recoils are thought to exert black hole kicks of order km s (Moody & Sigurdsson 2009). This value depends on metallicity, but for simplicity I assume a single, typical value here.
The retention fractions following from Eq. 8 are shown as a function of cluster mass per unit Plummer radius in Fig. 1. This quantity best reflects the retention fraction because in Eq. 8. Open clusters (with initial masses such that typically , Larsen 2004) do not retain any neutron stars or black holes, while globular clusters (—, Harris 1996) retain 0.1—4% of the neutron stars and 0.3—7% of the black holes. These values are in excellent agreement with other studies (e.g. Pfahl et al. 2002; Moody & Sigurdsson 2009), but are still lower than the large observed number of neutron stars in a number of globular clusters (the aforementioned ‘retention problem’).
3.2 Dissolution and the evolution of the mass function
Dissolution alters the shape of the stellar MF in star clusters due to the effects of twobody relaxation and energy equipartition. In a pioneering paper, Hénon (1969) derived the escape rate of stars of different masses from an isolated cluster. The cluster was represented by a Plummer (1911) gravitational potential:
(9) 
where denotes the Plummer radius setting the concentration of the cluster and represents the central potential, with the gravitational constant and the cluster mass. It was argued by Hénon (1960) that the only way for stars to escape such an isolated cluster is by a single, close encounter. The corresponding stellar massdependent escape rate was found to be (Hénon 1969):
(10) 
with the MF, the stellar mass, the total energy of the cluster, and a function related to the ejection probability for a star of mass in a close encounter with a star of mass and a corresponding mass ratio . The expression in Eq. 10 is independent of the adopted IMF. The function will be referred to as the ‘Hénon function’ and is shown in Fig. 2. The original expression consists of several integrals that have to be solved numerically. In Hénon (1969), a table is given for the Hénon function, but it can also be fitted by:
(11) 
This approaches the power law for , as was derived explicitly by Hénon (1969).
The total mass loss rate corresponding to Eq. 10 conflicts with body simulations (as was already noted by Wielen 1971) because only ejections by single, close encounters are included. This restriction implies that the disruption timescale is proportional to the crossing time (), while body simulations show that it scales with a combination of the halfmass relaxation time and the crossing time () due to twobody relaxation, i.e. the repeated effect of soft encounters (e.g. Baumgardt & Makino 2003). Nonetheless, the escape rate from Hénon (1969) does accurately describe what happens if two stars interact and can therefore be used as a starting point for a more complete description of the evolution of the MF. For that purpose, it is convenient to scale Eq. 10 to the mass loss rate found in body simulations and only use the relative or ‘differential’ stellar mass dependence from Hénon (1969). This is allowed if the ratio only depends on global cluster properties. It is straightforward to show (e.g. Spitzer 1987; Heggie & Hut 2003) that indeed this is the case as with the Coulomb logarithm. As such, one can write
(12) 
with the mass loss rate from Eq. 3 (Lamers et al. 2005) and the stellar massdependent escape rate per unit mass loss rate. The quantity is completely independent of the prescription for the total mass evolution. In order to derive , I start from Eq. 10 and express as
(13)  
where represents a correction factor to account for additional physics (see below). The numerator reflects the escape rate, while the denominator is proportional to the mass loss rate.
For mathematical simplicity^{3}^{3}3And because this is the only way to obtain an analytical solution as in Eq. 10. Hénon (1969) made the following assumptions in the derivation of Eq. 10.

The cluster exists in isolation and the tidal field is neglected. Therefore, ejection can only occur by a single, close encounter and the repeated effect of soft encounters (twobody relaxation) is not accounted for. This underestimates the escape rate of massive stars.

The distribution of stars is independent of stellar mass, i.e. mass segregation is not included. Depending on the balance between their encounter rate and their proximity to the escape energy, this over or underestimates the escape rate of lowmass stars from Hénon (1969). Considering the results from Baumgardt & Makino (2003), the latter seems to be the case.
The remainder of this section concerns the derivation of the factor in Eq. 3.2 that corrects for the above assumptions.
Let us assume that the distribution of stars over radius and velocity space is initially independent of their mass. This implies that mass segregation is dynamically created and not primordial, which is discussed in Sect. 6. For such an initial distribution, the separation from the escape energy is independent of mass. As the cluster evolves, energy equipartition is reached between the stars and the radius, velocity and proximity to the escape energy become a function of stellar mass. I first consider this effect on the escape rate before including the timescale on which twobody relaxation occurs for different stellar masses. Please note that the formulation of Eq. 3.2 with appearing in the numerator and the denominator implies that only the proportionality of is important. Its exact value is determined by constants that drop out when substituting in Eq. 3.2.
It is intuitive to express the dependence of the escape rate on the energy needed for escape as . The energy that is required for escape is related to the position and velocity of the star.^{4}^{4}4The energy difference that is discussed here concerns the energy that needs to be added to reach the escape energy. As such, it differs from the separation from the escape energy in Fukushige & Heggie (2000) and Baumgardt (2001), who are considering the excess energy of stars and its relation to the escape time, resulting in the aformentioned relation . For the potential in Eq. 9 it is given by
(15) 
with and the radial position and velocity of the star, and its escape velocity. If the cluster is in ‘perfect’ energy equipartition and correspondingly perfect mass segregation, the radius and velocity become a monotonous function of stellar mass (Heggie & Hut 2003, Ch. 16). Mass segregation is strongest in the cluster centre, which for a Plummer (1911) potential can be approximated with a harmonic potential . For a cluster in a tidal field the potential is truncated, and the harmonic approximation serves as a crude but reasonable approximation for the entire cluster (Heggie & Hut 2003, Ch. 16). Energy equipartition yields
(16) 
with the mean speed of all stars squared and the mean stellar mass. For the harmonic potential, this translates to a similar relation for the radial position:
(17) 
where represents the typical radius of the system, in this case the Plummer radius. This relation assumes that there is no particular stellar mass which dominates the mass spectrum. The decrease of radial position with stellar mass implied by Eq. 17 is a direct consequence of the energy loss endured by massive stars^{5}^{5}5And the energy gain experienced by lowmass stars. as the system evolves towards energy equipartition. Substituting Eqs. 16 and 17 into Eq. 15 and dividing out the proportionality gives an expression for :
(18) 
with denoting the ratio of the mean speed squared to the central escape velocity squared. This constant mainly depends on the degree of mass segregation. Consequently, it will depend on the IMF. By comparing the models to the body simulations with a mass spectrum by Baumgardt & Makino (2003) the value is constrained to for a Kroupa IMF, using King (1966) potentials with King parameter —7 (see Sect. 4). For reference, an unevolved Plummer (1911) potential has .
By comparing the models to body simulations (provided by M. Gieles, private communication) with different IMF power law slopes and a ratio between the maximum and minimum mass of 10, the approximate relation is found^{6}^{6}6This prescription for implies that the condition for the stars in the cluster to be physically bound is satisfied for all . for a MF . Fitting the Kroupa IMF with a single power law in the mass range 0.08— (as used by Baumgardt & Makino 2003) yields , resulting in as mentioned earlier.^{7}^{7}7Nonetheless, the relation for should be expected to exhibit some variation for different mass ranges. The comparison with body simulations also showed that a single value of suffices to determine the MF evolution, even though it does not remain constant over the full cluster lifetime.
Because , Eq. 18 indicates that the escape rate of lowmass stars is increased if a cluster is in complete energy equipartition. However, the timescale on which twobody relaxation occurs between different stellar masses has not yet been considered. For a cluster starting with a stellar massindependent distribution of radial positions and velocities, the equipartition timescale scales as
(19) 
for equipartition between stars of masses and (Heggie & Hut 2003). This is a modified version of the relaxation timescale, which shows a very similar proportionality (). It illustrates that twobody relaxation occurs on a shorter timescale for massive stars than for lowmass stars, increasing their escape rate .
The correction factor for the escape rate that appears in the integrals of Eq. 3.2 now follows from Eqs. 18 and 19 as
(20)  
It was mentioned before that the proportionalities of and rather than their exact values suffice for the computation of due to the renormalisation of the total mass loss rate that appears in Eq. 3.2: only the stellar massdependence is important.
The influence of the tidal field is now included in two ways. First of all, the ejection of stars no longer occurs by a single, close encounter but arises due to twobody relaxation on the equipartition timescale, representing the repeated effect of soft encounters. Secondly, the above derivation of the separation from the escape energy assumes a potential which approximates tidally limited clusters. As a result, the escape rate of massive stars is increased with respect to clusters in the model of Hénon (1969), which was derived for an isolated cluster. On the other hand, the effect of mass segregation is included by introducing a stellar massdependence for the energy needed by stars to reach the escape velocity. Lowmass stars are closer to the tidal radius than massive stars, leading to a lower energy that is needed for escape and an increased escape rate. It depends on the shape of the MF which mechanism dominates.
The evolution of the MF of various cluster components is obtained from Eqs. 12, 3.2 and 20 by writing
(21) 
where the MFs of stars, white dwarfs, neutron stars and black holes are represented by , with . The overall cluster evolution is computed by combining the results of this section with the prescription for stellar evolution from Sect. 3.1.
If stellar evolution is included, the resulting mass loss causes an expansion of the cluster, during which stars are lost independently of their masses. This delays the onset of mass segregation and the stellar massdependent mass loss that is described above. The moment of transition to stellar massdependent mass loss can be characterised by a certain fraction of the initial cluster mass that has been lost by dissolution . It is assumed that the fraction of the mass loss for which the ejection rate depends on the stellar mass grows exponentially^{8}^{8}8This form assumes that the increase of the fraction of the mass loss that is stellar massdependent scales with the total dynamical mass loss, which is a compromise between a step function and a linear increase. between 0 and 1 as
(22) 
where the subscript ‘smd’ denotes ‘stellar massdependent’, is the fraction of the initial mass that has been lost by dissolution at which mass segregation is reached, and is a constant to normalise at the reference value . For , per definition , indicating that all mass loss is stellar massdependent. The timescale on which mass segregation is reached and the transition to stellar massdependent mass loss is completed is proportional to the initial halfmass relaxation time (). It has been shown in several studies that for Roche lobefilling clusters the disruption timescale (Vesperini & Heggie 1997; Baumgardt & Makino 2003; Gieles & Baumgardt 2008), implying that . The expression for in Eq. 2 then leads to . Assuming that the cluster mass evolution is close to linear, the firstorder relation is obtained, implying
(23) 
for a King parameter of , with the dissolution timescale at the solar galactocentric radius Myr. For a King parameter of , the exponent of the initial cluster mass becomes 0.23 and Myr (Kruijssen & Mieske 2009). In this relation, represents a constant that is fixed by comparing the model to the results of body simulations from Baumgardt & Makino (2003), giving for and for (see Sect. 4). The variation with King parameter arises because twobody relaxation is faster for more concentrated clusters. If stellar evolution were neglected, at all ages and .
The modeled MF slope change in the mass range — is shown in Fig. 3 for different IMFs covering —. Evidently, is a function of the remaining mass fraction and is insensitive to the slope of the IMF, as long as that the ratio between the maximum and minimum mass is kept fixed and stellar evolution is excluded. This is an interesting observation in view of the MF evolution of globular clusters, in which — and stellar evolution only plays a minor role. Figure 3 shows that the slope of the MF in globular clusters could be a possible indicator for the mass fraction that has been lost due to dissolution, provided that the IMF does not vary and the remnant retention fractions were not substantially dissimilar during the early evolution of different globular clusters (see Sect. 5.2 and Fig. 19).
For the particular example of a Kroupa MF that is truncated at different maximum masses , the relative escape rate per unit mass loss rate (see Eqs. 12 and 3.2) is shown in Fig. 4. This quantity is proportional to and reflects the probability that a star of a certain mass is ejected. Figure 4 illustrates that the mass of the highest relative ejection rate is related to the maximum mass of the MF. The peak occurs at intermediate masses if the MF is truncated at a high mass. This implies that there is a typical mass where the stars are not too far from the escape energy and have an equipartition timescale with the massive stars that is short enough to eject them efficiently. This ‘sweet spot’ depends on the maximum mass of the MF. If the MF is truncated at an intermediate mass, the combination of quick twobody relaxation and proximity to the escape energy favours the ejection rate of stars at the lowest masses.
The maximum stellar mass at which the transition from ‘sweet spot’depletion to lowmass star depletion happens, is determined by the proximity of the lowmass stars to the escape energy. In Fig. 5, the mass of the peak relative ejection rate is shown as a function of the maximum stellar mass. At low truncation masses, the peak occurs at the minimum mass, indicating strong lowmass star depletion. Around , the relative ejection rate at becomes larger than its value at the lowest masses, which causes a jump in Fig. 5. For even higher values of , the peak relative ejection rate typically occurs at 15—20% of the maximum mass, approximately following the relation
(24) 
Even though its quantitative properties only hold for a Kroupa MF, the variation of the relative escape rate with the maximum mass of the MF has several implications for star cluster evolution. The change of in Figs. 4 and 5 can be interpreted as an example of what happens when stellar evolution removes the most massive stars in the cluster, provided that the remnants are all ejected by their kick velocities. If dynamical evolution does not affect the shape of the MF too much before , or Myr, the subsequent evolution of the MF will be dominated by lowmass star depletion. If substantial dissolution occurs earlier on, it is dominated by the ‘sweet spot’ depletion of intermediate masses. Only the retention of massive stellar remnants will make the evolution of the MF deviate from these basic estimates, because remnant retention can provide a fixed maximum (remnant) mass of the MF. This is treated in more detail in Sect. 5.
4 Comparison to body simulations
The model described in Sect. 3 can be easily verified by running it for the exact same boundary conditions as the body simulations^{9}^{9}9These were performed using NBODY4 (Aarseth 1999). by Baumgardt & Makino (2003) and comparing the results. They conducted simulations of Roche lobefilling clusters between 8k and 128k particles, which were evolved in the Galactic tidal field at galactocentric radii in the range 2.833—15 kpc. The boundary conditions for the body runs of Baumgardt & Makino (2003) differ from those described in Sect. 3 by neglecting kick velocities and defining the Kroupa stellar IMF between 0.1 and 15 M, thereby excluding black holes. For this particular comparison, the same IMF, stellar evolution prescription, and initialfinal mass relation for stellar remnants are used in the model that is presented in this paper.
In Fig. 6, the modeled evolution of the (luminous) stellar MF is compared to the body runs with King parameter for a range of cluster masses and total disruption times. As time progresses, the maximum stellar mass decreases due to stellar evolution and the MF is lowered due to the dynamical dissolution of the star cluster. The slope of the MF changes due to the preferential ejection of lowmass stars, which have energies closer to their escape energies, even to the extent that it dominates over their relatively slow twobody relaxation. For both the models and the body simulations, the MF develops a slight bend at when approaching total disruption. The bend arises as an optimum between on the one hand high energies but slow relaxation for the lowest stellar masses, and on the other hand quick relaxation but low energies for the highest stellar masses (see the discussion at the end of Sect. 3).
In all cases, the resemblance of the models and the body simulations is striking. The models reproduce all key aspects of the body runs, such as the amount of lowmass star depletion, the changing slope at for clusters close to dissolution, the survival of the Kroupa bend at , and the dependence of the lowmass depletion on the total lifetime of the cluster (compare the three 32k runs). The only difference occurs at the highmass end of the MF, where the maximum stellar masses do not match at young ages. This is due to a minor dissimilarity of the total mass evolution (also see Lamers et al. 2005; Kruijssen & Lamers 2008). Because the maximum stellar mass only depends on the age of the cluster, this causes a difference in maximum stellar mass when showing the MFs at fixed remaining cluster mass fractions. The contrast is clearest at young ages, since there the maximum stellar mass most rapidly decreases.
In the description of the model in Sect. 3, two constants have been determined from the body simulations by Baumgardt & Makino (2003). These constants are the ratio of the mean speed squared to the central escape velocity squared (, see Eq. 18) and the proportionality constant for the relation marking the transition to stellar massdependent mass loss (, see Eq. 23). As mentioned in Sect. 3, for a Kroupa IMF and King parameter one obtains and . To illustrate the robustness of the models, in Fig. 7 they are compared to a 64k body run with . For such a cluster with a higher concentration, the early mass segregation implies . Again, the model and the simulation are in excellent agreement.
The dependence of the MF evolution on both constants is considered in Fig. 8. For , the dependence of the evolution of the MF on its value is shown in the upper panel of Fig. 8, while for it is shown in the bottom panel of Fig. 8. Both panels show the evolution of the MF for the 64k cluster in Fig. 6 for different values of and .
The ratio of the mean speed squared to the central escape velocity squared affects the ejection probability of the stars with the lowest masses. Because these stars are closest to their escape energies in a masssegregated cluster, they are most strongly influenced by the value of . For higher , the MF gets more depleted in lowmass stars due to their closer proximity to the escape energy, while for lower more lowmass stars are retained as the balance between close proximity to the escape energy and slow relaxation shifts to the latter.
The proportionality constant for the transition to stellar massdependent dissolution in Eq. 23 affects the MF as a whole. For lower , the transition occurs earlier and more lowmass stars are lost, while for higher the onset of the depletion is delayed and the slope of the MF remains closer to its initial value. If one were to assume a constant , which is contrary to the adopted relation with cluster mass in Eq. 23, this would therefore yield a stellar MF in massive clusters that is underpopulated in lowmass stars, and a MF in lowmass clusters that is overabundant in lowmass stars.
5 Star cluster evolution
In this section, the described model is applied to compute the evolution of clusters for a variety of boundary conditions. The stellar content as well as integrated photometry are addressed, using the boundary conditions from Sect. 3 instead of those that were adopted to compare the model to body simulations in Sect. 4. The most important differences are the mass range of the IMF, the inclusion of remnant kick velocities, and the initialfinal mass relation.
The model that will be referred to as the ‘standard model’ uses a metallicity (which is typical of globular clusters), a King parameter^{10}^{10}10For , or , the results vary only marginally. of (corresponding to in Eq. 2), a dissolution timescale parameter Myr, and a Kroupa IMF between and the maximum stellar mass given by the Padova isochrones at , which is typically . For the computation of the retained remnant fraction (see Eq. 8), the Plummer radius is related to the halfmass radius as . The halfmass radius is assumed to remain constant during the cluster lifetime (e.g. Aarseth & Heggie 1998). For the relation between and initial cluster mass the expression from Larsen (2004) is adopted:
(25) 
The models that are used in this section are computed from yr to yr (the maximum age of the Padova isochrones) for initial masses between and , spaced by 0.25 dex intervals.
5.1 The influence of the disruption time
The disruption time of a cluster affects the evolution of the MF and of the integrated photometric properties. To assess the influence of the disruption time on cluster evolution, clusters with low and high remnant retention fractions should be treated separately, because the presence of massive remnants also has a pronounced effect on the results (see Sect. 5.2). As shown in Fig. 1, for a given kick velocity dispersion the remnant retention fraction is set by the cluster mass. This means that the division between low and high remnant retention fractions can be made by making a cut in initial cluster mass.
In Fig. 9, the impact of the disruption time on the evolution of the MF is shown for a cluster with initial mass , representing the evolution for low remnant retention fractions.^{11}^{11}11High remnant retention fractions will be treated in the discussion of the influence of the retention fraction in Sect. 5.2. The range of the dissolution timescale parameter and resulting total disruption times that are considered in Fig. 9 cover two orders of magnitude. As the total lifetime increases, the depletion of the lowmass stellar MF close to total disruption becomes more prominent. Conversely, the MF of shortlived clusters is depleted around . As introduced in the last paragraphs of Sect. 3, this difference is caused by the fixed timescale on which stellar evolution decreases the maximum stellar mass, implying that the masses of the most massive stars are larger in quickly dissolving clusters than in slowly dissolving ones. Because in shortlived clusters the massive stars are still present when the bulk of the dissolution occurs, their rapid twobody relaxation with intermediatemass stars dominates over the relatively close proximity to the escape energy of lowmass stars, yielding a depletion at intermediate masses. In longlived clusters, this cannot occur because the very massive stars have disappeared before the mass loss by dissolution becomes important, thus resulting in the depletion of the very lowmass end of the MF. As a rule of thumb, for Myr (which is the lifetime of a star) the depletion typically occurs around 15—20% of the mass of the most massive star (see Sect. 3). In terms of the total disruption time, the transition from intermediatemass star depletion to lowmass star depletion occurs around Gyr.
A quantifiable way to look at the evolution of the stellar MF in star clusters is to consider the slope of the MF in certain mass intervals (Richer et al. 1991; De Marchi et al. 2007; De Marchi & Pulone 2007; Vesperini et al. 2009). For the commonly used mass intervals () and (), Fig. 10 shows the evolution of the slope for the same clusters as before. Like Fig. 9, this illustrates that for short disruption times the slope steepens as the cluster dissolves, while for long disruption times the slope flattens with time. The presented models and other model runs indicate that increases with time for Gyr and decreases for Gyr. For total disruption times in between these values, the slope first increases and then decreases. The slope in the second mass interval shows the same behaviour. It increases for Gyr and decreases for Gyr.
The masstolight () ratio evolution of star clusters is affected by the evolution of the MF due to the large variations in ratio between stars of different masses. Massive stars have lower ratios than lowmass stars, implying that a cluster with a MF that is depleted in lowmass stars will have a reduced ratio (Baumgardt & Makino 2003; Kruijssen 2008; Kruijssen & Lamers 2008). As such, one would also expect a correlation between the slope of the MF and ratio.
In Fig. 11, the evolution of the ratio of the band to the masstolight ratio due to stellar evolution is shown for the same clusters as in Figs. 9 and 10. This quantity reflects the relative ratio change induced by dynamical evolution with respect to evolutionary fading only. If the ejection rate would be independent of stellar mass, the evolution would follow a horizontal line at . However, when accounting for dynamical evolution, the ratio is always smaller than that for stellar evolution only. Somewhat surprisingly, this is also the case for clusters for which the slope of the MF increases (see Fig. 10). This is explained by looking at the evolution of the entire MF in Fig. 9. Even though the slope at low masses increases for short disruption times due to the ejection of intermediatemass stars, the most massive stars that dominate the cluster light are still retained. Because stars of intermediate masses are lost instead, the ratio decreases.
Because the slope of the stellar MF either increases or decreases at masses , the decrease of the ratio implies a large range of MF slopes that can occur at low ratios. This is shown in Fig. 12, where the relation between and the ratio drop is presented. The slope of the stellar MF in a certain mass range does not necessarily reflect the ratio of the entire cluster. Considering the aforementioned rule of thumb stating that for total disruption times Gyr the depletion of the MF occurs around 15—20% of the mass of the most massive star , it is useful to define the slope in a mass range that is related to . In Fig. 12, the relation between slope and ratio is also shown for the slope in the stellar mass range . In such a relative mass range, the slope follows a much narrower relation with ratio. The range between 30% and 80% of was chosen to maximise this effect.
For the slopes in the fixed stellar mass ranges ( and , see above), the relation with the ratio becomes better defined for longlived clusters. It is shown in Figs. 10—12 that both the slope and the ratio decrease for clusters with long disruption times, indicating that both quantities are more clearly related for globular clusterlike lifetimes.
The colour of star clusters is also influenced by the evolution of the MF, due to the colour differences between stars of different masses. The magnitude difference with respect to the value that a cluster would have if dynamical evolution were neglected is shown in Fig. 13. As the clusters dissolve, their colours become redder due to the ejection of main sequence stars. The magnitude difference in exceeds mag for total disruption times Gyr. In redder passbands (e.g. the colour), the difference grows to several tenths of magnitudes. For longer total disruption times only stars of the lowest masses are ejected (see Fig. 9), which hardly contribute to the cluster light and colour, implying that the colours are only marginally affected.
5.2 The influence of the remnant retention fraction
The formation of stellar remnants introduces massive bodies in the MF that do not end their lives due to stellar evolution like massive stars do. Depending on their kick velocities, stellar remnants can be retained in (massive) clusters. If they are retained, they keep affecting the evolution of the stellar MF until the cluster is disrupted. Especially black holes can have a pronounced effect on cluster evolution.
The remnant retention fraction arises from the cluster mass, radius and the kick velocity dispersion (see Eq. 8). In this section, the massradius relation from Eq. 25 is used. Although the results will differ for other relations, it has been verified that for commonly used alternatives,^{12}^{12}12Such as a constant radius or density. the change is only marginal and does not affect the nature of the conclusions. To separate the effect of remnant retention from that of the disruption time, a fixed initial cluster mass of is assumed while independently varying the velocity dispersion of the remnant kick velocities and the disruption time. The corresponding evolution of the stellar MF is shown in Fig. 14, for the standard model (see the beginning of this section) with black hole kick velocity dispersions km s, equivalent to for a cluster, and for dissolution timescale parameters Myr, which for a cluster implies Gyr. Assuming an age of 12 Gyr, the presentday mass in the case of Myr is about , comparable to globular clusters. The remaining fraction of the initial mass is .
If the velocity dispersion of black hole kicks is low and a relatively large fraction of black holes is retained, then the ejection rate of massive stars is increased with respect to high kick velocity dispersions. This arises due to the quick twobody relaxation between the massive stars and the black holes, which will have masses larger than the most massive stars after a few Myr of stellar evolution. As a result, the ejection rate of lowmass stars is largest in clusters containing only few black holes. This happens for clusters with either long or short disruption times, but the effect is largest for longlived clusters (the solid lines in Fig. 14). In these clusters the maximum stellar mass is more strongly decreased by stellar evolution than in shortlived clusters, implying that the black hole masses are larger compared to the most massive stars in these clusters. For long disruption times, the presence of massive remnants therefore has a more pronounced effect on the ejection rate of massive stars than for short disruption times. If these longlived clusters retain a sufficiently large fraction of the stellar remnants, their stellar MF may even become depleted in massive stars.
The top panel of Fig. 14 also shows that for a cluster with a high remnant retention fraction, the impact of the disruption time on the MF evolution is similar to that of clusters with low retention fractions (see Fig. 9). However, the influence of the disruption time becomes smaller when more remnants are retained. This explains why Baumgardt & Makino (2003) only found a very weak dependence of the evolution of the MF on the disruption time (also see Fig. 6), since they neglected remnant kick velocities and retained all remnants in their simulations.
Analogous to Fig. 10 in Sect. 5.1, the evolution of the MF slope in different mass ranges is shown in Fig. 15 for the clusters with Myr from Fig. 14.^{13}^{13}13For the clusters with relatively long disruption times that are considered in this section, the variable stellar mass range that was introduced in Sect. 5.1 to trace the relation between MF slope and ratio gives an evolution of the slope that is comparable that for the fixed mass ranges. It is omitted from the figures in this section to improve their clarity. The kick velocity dispersion has an effect on the MF that is more uniform than the consequences of variations in the disruption time, leading to very similar slope evolutions in the two different stellar mass ranges. Independent of the mass range, an increase in remnant retention fraction is reflected by an increase of . The model that is displayed for km s, Myr, and (the upper dashed and solid lines in Fig. 15) marks the transition between an increase or decrease of the MF slope by dynamical evolution. For an initial , lowmass stars are preferentially ejected during most of the cluster lifetime, while for mainly the massive stars escape. For shorter disruption times, the transition is located at a smaller black hole retention fraction.
Because the black hole retention fraction affects the overall slope of the stellar MF, the changes in are matched by corresponding changes in the ratio. In Fig. 16, the relative ratio change due to dynamical evolution is shown for same clusters as in Fig. 15. Contrary to the clusters with low remnant retention fractions in Sect. 5.1, the ratio of the clusters in Fig. 16 does not monotonously decrease. Close to total disruption, the massive remnants are the last bodies to be ejected. During that short phase of cluster evolution, the ratio is increased by dynamical evolution and exceeds the value it would have due to stellar evolution alone.
The behaviour of ratio for different black hole kick velocity dispersions has interesting implications for the relation between stellar MF slope and ratio, which is shown in Fig. 17. In combination with Fig. 12 (note the different axes), it shows possible evolutionary tracks of star clusters in this plane, indicating that nearly every location may be reached. However, when limiting ourselves to longlived clusters, Fig. 17 illustrates that these clusters will follow a trend of decreasing slope with decreasing ratio, albeit with excursions to high ratios and slightly higher close to their total disruption. This explains the trend that was found by Kruijssen & Mieske (2009), who considered the relation between the observed MF slopes and ratios of Galactic globular clusters.
The colour change due to dynamical evolution is only very small for clusters with Gyr (see Sect. 5.1). Because clusters in which remnants are retained are massive, their lifetimes are correspondingly long. As a result, the colour evolution is largely unaffected for the clusters in which the remnant retention fraction could play a role ( mag). The colour change is even smaller if more massive remnants are retained, because then the stellar MF more closely resembles its initial form (see the upper panel of Fig. 14). Longlived clusters generally appear mag bluer in due to dynamical evolution during the last —20% of their lifetimes and reach a similar reddening upon their total disruption, which is well within observational errors. The colours of old clusters are thus only marginally affected by dynamical evolution.
The evolution of the total remnant mass fraction is shown in Fig. 18 for different black hole kick velocity dispersions. The seemingly counterintuitive result is that the fraction of the cluster mass that is constituted by remnants is smaller when more black holes are retained. As shown in Fig. 14, the retention of black holes suppresses the depletion of the lowmass end of the MF due to the ‘sweet spot’ ejection (see Sect. 3) of massive () stars by the black holes. After Gyr, white dwarfs and neutron stars have masses that are similar to those of the massive stars, implying that their ejection rate is also increased when more black holes are retained. Because the total mass constituted by white dwarfs and neutron stars is larger than the combined mass of all black holes, the fraction of the total cluster mass that is constituted by remnants decreases if these lowmass remnants are ejected by the more massive black holes.
6 Discussion and applications
The results of this paper show that the stellar MFs in star clusters differ strongly from their initial forms due to dynamical cluster evolution. The specific kinds of these differences depend on the properties of the star clusters and their tidal environment, most importantly on the disruption time, remnant retention fraction, and IMF.^{14}^{14}14Although not specifically shown in this paper (but not surprisingly), the differences also depend on the initialfinal stellar mass relation.
A physical model for the evolution of the stellar MF is presented in which twobody relaxation leads to a stellar mass dependence of the ejection rate. For any particular stellar mass, the ejection rate is determined by the typical proximity of that mass to the escape energy and by the timescale on which the twobody relaxation with the other stars takes place. Combined with a prescription for stellar evolution, stellar remnant production, and remnant retention using kick velocity dispersions, this provides a description for the total evolution of the MF. This description is independent of the adopted total mass evolution. The model shows that the slope of the mass function is a possible indicator for the mass fraction that has been lost due to dissolution, provided that the IMF does not vary and the remnant retention fraction has been fairly similar for young globular clusters.^{15}^{15}15Any variability of the retention fraction would induce substantial scatter, see Sect. 5.2 and Fig. 19.
For the exact same initial conditions, the model shows excellent agreement with body simulations of the evolving MF by Baumgardt & Makino (2003). However, an important advantage of the presented model compared to the (more accurate) body simulations is its short runtime and corresponding flexibility. It can be easily applied to compute the evolution of clusters for a large range of initial conditions. The results can then be used to identify interesting cases for more detailed and less simplified calculations with body or Monte Carlo models.
The most important simplification of the model is neglecting the effect of binary encounters on the stellar mass dependence of the ejection rate. To incorporate binaries, a conclusive census of the binary population in star clusters would be required, which is not yet available. Nonetheless, it is possible to make a qualitative estimate for the effect of binaries. The encounter rate of binaries would typically be higher than that of individual stars, because the cross section of binaries is larger. This would increase the relative escape rate at the stellar mass for which the binary fraction^{16}^{16}16The fraction of stars residing in binary or multiple systems. peaks. This binary fraction is found to increase with primary mass (see e.g. Kouwenhoven et al. 2009). Because massive stars are removed by stellar evolution, this implies that the binary fraction decreases with age, which is in agreement with the low binary fraction observed in globular clusters (, e.g. Richer et al. 2004). The effect of binaries on the evolution of the mass function would thus be most notable if the majority of the dynamical mass loss occurs at ages Myr (the typical lifetime of an star), in which case it would somewhat enhance the relative escape rate of the most massive stars.
The model is applied to investigate the influence of the disruption time and remnant retention on the evolution of the MF and integrated photometric properties of star clusters. For total disruption times Gyr, the modeled relative ejection rate is highest at a certain ‘sweet spot’ mass that is typically 15—20% of the mass of the most massive objects in the cluster. For longer lifetimes, the evolution of the MF is dominated by lowmass star depletion, unless the retention fraction of massive stellar remnants is larger than 0.25. Only in the particular case of such a high retention fraction, the ratio is increased by dynamical evolution when the cluster approaches total disruption. In all other scenarios, the ratio decreases because the most massive (luminous) stars are kept.^{17}^{17}17This process differs from a possible variability of the proportionality between the velocity dispersion and the cluster mass, which concerns a much shorter timescale (e.g. Boily et al. 2009). When defining the slope of the MF in the range 30—80% of the maximum stellar mass, this gives a clear relation between the MF slope and the ratio. For slopes that are defined in fixed mass ranges, there is not necessarily a correlation between slope and ratio if Gyr. In clusters with a longer total disruption time, both quantities are related. Dynamical cluster evolution is found to induce some reddening of the integrated cluster colours, amounting up to 0.1—0.2 mag in for total disruption times Gyr. The fraction of the cluster mass that is constituted by remnants surprisingly decreases if more black holes are retained, because the black holes preferentially eject bodies around the masses of white dwarfs and neutron stars, which contain most of the total remnant mass.
Contrary to what is suggested by other studies (e.g. Baumgardt & Makino 2003; Anders et al. 2009), the evolution of the MF is not homologous. The reason that these studies concluded that its evolution is very similar for all clusters (also see Figs. 6 and 7), is that they assumed that all remnants were retained. It is illustrated in Fig. 14 that the differences between clusters with dissimilar disruption times disappear when the retention fraction increases. For realistic retention fractions, differences do arise. If two clusters with different initial masses have the same total disruption time, their MF evolution will be dissimilar due to their different remnant retention fractions and the impact of the retained remnants on the dynamical cluster evolution. Alternatively, if two clusters have equal initial masses but different total disruption times, for instance due to differences in their galactic location or environment, their MF evolution will be dissimilar due to the dynamical impact of the evolution of the maximum stellar mass.
The larger variation of MF evolution that is found with presented model may also be able to explain observations of globular clusters in which the MF cannot be characterised by a single power law (De Marchi et al. 2000). If the evolution of the MF were homologous, these features would likely be primordial (Baumgardt & Makino 2003), but this is not necessarily the case when using realistic retention fractions. Most other differences between the results presented in Sect. 5 and those from Baumgardt & Makino (2003) are also due to their assumption of full remnant retention. For example, their ratio evolution shows a smaller decrease than in Fig. 11. This is explained in Fig. 16, where it is shown that dynamical evolution reduces the ratio by a smaller amount if the retention fraction is larger.
Studies on the fractal nature of cluster formation show that star clusters are initially substructured (Elmegreen 2000; Bonnell et al. 2003). Even though this substructure is typically erased on a crossing time, it can induce primordial mass segregation in star clusters (McMillan et al. 2007; Allison et al. 2009). The influence of primordial mass segregation on the evolution of the MF has recently been investigated by Baumgardt et al. (2008) and Vesperini et al. (2009). While Baumgardt et al. (2008) do not include stellar evolution and concentrate on twobody relaxation, Vesperini et al. (2009) do include stellar evolution. They show that for some degrees of primordial mass segregation, the mass loss by stellar evolution can induce additional dynamical mass loss that strongly decreases the total disruption time. For clusters that survive for a Hubble time, the MF evolution in the case of primordial mass segregation is very similar to an initially unsegregated cluster. Vesperini et al. (2009) conclude that the evolution of the MF is only affected by primordial mass segregation for clusters in which the total disruption time is sufficiently decreased by the induced mass loss. In that case, the slope of the MF remains much closer to its initial value than it would in clusters without primordial mass segregation. Their conclusion is consistent with the model presented in this paper, because the evolution of the MF is determined by the most massive stars at the time when the largest mass loss occurs (see Figs. 4 and 9). This induced mass loss enters the model in terms of the absolute mass loss rate in Eq. 3, not in the stellar massdependent escape rate per unit mass loss rate of Eq. 3.2.
A change in total mass loss rate is not the only consequence of primordial mass segregation. Baumgardt et al. (2008) have shown that lowmass star depletion is enhanced for clusters without stellar evolution that are primordially masssegregated. This occurs because energy equipartition is reached on a shorter timescale and because of their use of a fixed () maximum stellar mass. As a result, there are no massive bodies to increase the ejection rate of intermediate mass stars (see Fig. 5), implying that only the lowmass stars are preferentially lost. In the present paper, mass segregation is assumed to arise dynamically, but the model could in principle be adapted to cover primordial mass segregation by setting and adjusting to the initial velocity distribution until it is erased by dynamical evolution (see Eq. 23), after which the values from Sect. 3 can be used.^{18}^{18}18As explained in Sect 3, represents the ratio of the mean speed squared to the central escape velocity squared that depends on the degree of mass segregation (and thus on the IMF). On the other hand, is a proportionality constant in the expression for the onset of the stellar massdependent ejection of stars, which depends on the concentration or King parameter. This does not necessarily yield enhanced lowmass star depletion for clusters with a complete IMF (including masses ) because of the presence of massive stars or remnants.