{'fraction_non_alphanumeric': 0.03694295253302366, 'fraction_numerical': 0.015556200706440218, 'mean_word_length': 4.383563427976036, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical friction in multi-component evolving globular clusters

We use the Chandrasekhar formalism and direct N-body simulations to study the effect of dynamical friction on a test object only slightly more massive than the field stars, orbiting a spherically symmetric background of particles with a mass spectrum. The main goal is to verify whether the dynamical friction time (t_DF) develops a non-monotonic radial-dependence that could explain the bimodality of the Blue Straggler radial distributions observed in globular clusters. In these systems, in fact, relaxation effects lead to a mass and velocity radial segregation of the different mass components, so that mass-spectrum effects on t_DF are expected to be dependent on radius. We find that, in spite of the presence of different masses, t_DF is always a monotonic function of radius, at all evolutionary times and independently of the initial concentration of the simulated cluster. This because the radial dependence of t_DF is largely dominated by the total mass density profile of the background stars (which is monotonically decreasing with radius). Hence, a progressive temporal erosion of the BSS population at larger and larger distances from the cluster center remains the simplest and the most likely explanation of the shape of the observed BSS radial distributions, as suggested in previous works. We also confirm the theoretical expectation that approximating a multi-mass globular cluster as made of (averaged) equal-mass stars can lead to significant overestimates of t_DF within the half-mass radius.


Introduction
Blue Straggler Stars (hereafter BSSs) are core hydrogen-burning objects that, in the color-magnitude diagram (CMD) of globular clusters (GCs) and other stellar systems, populate a region along the extrapolation of the main sequence (MS) towards colors bluer than the turn-off point (see e.g., Sandage 1953;Ferraro et al. 2009;Mathieu & Geller 2009;Monelli et al. 2012). This observational evidence and direct measurements imply that the mass of BSSs is larger (in GCs, m BSS ≃ 1.21.4M ⊙ ; Shara et al. 1997;Ferraro et al. 2006a;Lanzoni et al. 2007a;Fiorentino et al. 2014) than that of typical stars (m ave ≃ 0.4M ⊙ in GCs) 1 . In the last decades, the ratio between the number of BSSs and that of a reference population (as red giant or horizontal branch stars) has been derived as a function of radius for more than 20 Galactic GCs and three different shapes of the BSS radial distribution have been observed (see Ferraro et al. 2012): a flat distribution (with constant value throughout the cluster), a bimodal trend (with a central peak, followed by a minimum and an outer rising branch) and a unimodal behavior (with a central peak and a monotonic decrease outward). Some of the features of the observed BSS radial distributions (as the radial position of the minimum in the case of bimodal trends) have been found to be well reproduced in terms of dynamical friction (DF) effects, even if several simplifying assumptions were used to describe the process (see Mapelli et al. 2004Mapelli et al. , 2006Lanzoni et al. 2007a,b).
Building on these considerations, Ferraro et al. (2012) proposed a comprehensive interpretative scenario identifying in the DF process the primary cause of the whole variety of shapes of the observed BSS radial distributions (e.g., Ferraro et al. 1993Ferraro et al. , 1997Ferraro et al. , 2004Ferraro et al. , 2006bLanzoni et al. 2007a,b;Dalessandro et al. 2008;Beccari et al. 2011). These authors linked the three classes of BSS radial profiles to the cluster dynamical age, thus defining the so-called dynamical clock for GCs: a flat, bimodal, and unimodal BSS radial distribution corresponds, respectively, to a dynamically-young, -intermediate, and -old cluster, with the position of the minimum of the distribution acting as a clock hand and marking the stage of internal dynamical evolution reached by the system. The idea behind this scenario is that the position of the minimum corresponds to the distance from the center where the DF time (t DF ) equals the system age. Since, in general, t DF monotonically increases with radius, such a minimum is not yet dug in dynamically-young clusters (

with a flat BSS distribution), it is found at intermediate radii in systems with intermediate dynamical age (displaying a bimodal BSS distribution), and it already reached the farthest outskirts in dynamically-old clusters (characterized by a unimodal BSS distribution).
While the proposed scenario is simple and seems to properly account for the observational data, one cannot exclude other more complicated, yet realistic, possibilities, given the intrinsic complexity of DF. In particular the present study aims at answering the following question: could the DF time-scale develop a non-monotonic radial behavior at some time during the cluster evolution? In fact, if t DF develops a minimum at a given radius r min , the BSSs orbiting at that distance from the center would suffer from an enhanced drag force with respect to the other BSSs orbiting at different radii, and a minimum in the BSS radial distribution would therefore appear at r min (i.e., at a place and time different from what expected by the scenario depicted above, therefore affecting the possibility to use the distribution of BSSs as a simple clock). Note that, while t DF is known to be a monotonically increasing function of radius in most of the astrophysical problematics investigated so far, this could not be the case for BSSs in GCs, because of three combined conditions that could play a relevant role, namely: (i) a test particle only slightly more massive than the field objects, (ii) a background field made up of stars with different masses, and (iii) an evolving and radially dependent field mass function, due to the tendency to energy equipartition. In this paper we aim at investigating the possibility that t DF is not a monotonic function of radius, by using direct N-body simulations and some analytical results pertinent to DF in the presence of a mass spectrum.
In general terms, DF can be described as the slowing-down of a body (the test particle) moving in a sea of background particles, due to the cumulative effect of long-range interactions (e.g., Binney & Tremaine 1987;Spitzer 1987;Bertin 2001). The processes underlying DF have been studied notably by Chandrasekhar (1943), in the case of a test particle moving in an infinite and homogeneous background field. In this view, DF is interpreted as the result of the non-zero vectorial sum of the changes in the parallel component of the relative velocity caused by all the encounters of the test particle with the field particles (each of them treated as an isolated two-body encounter), as a consequence of energy conservation on each relative orbit. Remarkably, in the case of an isotropic velocity distribution function, and under the assumption of impulsive approximation, it can be shown that (at the leading order) only field particles moving slower than the test particle contribute to the deceleration of the test particle. Many authors have extended the Chandrasekhar's framework to other more realistic physical cases, with different assumptions or different and more sophisticated methods (Chandrasekhar & von Neumann 1942, 1943White 1948;Thorne 1968;Lee 1969;Binney 1977;Tremaine & Weinberg 1984;Ostriker 1999;Ciotti & Binney 2004;Nipoti et al. 2008). From the astrophysical point of view, it is clear that DF plays an important role on different scales, from galaxy clusters and their cD galaxies (e.g., Ostriker & Tremaine 1975;White 1976;Binney 1977;Dressler 1979;Kashlinsky 1987;Nipoti et al. 2004;Kim et al. 2005;El-Zant 2008), to galaxies and their GC systems (e.g., Tremaine, Ostriker & Spitzer 1975;Bontekoe & van Albada 1987;Bertin et al. 2003Bertin et al. , 2004Capuzzo-Dolcetta & Vicari 2005;Arena et al. 2006;Mastrobuono-Battisti & Capuzzo-Dolcetta 2012;Arca-Sedda & Capuzzo-Dolcetta 2014), to binary black holes at the center of early type galaxies (e.g., Fukushige et al. 1992;Vecchio et al. 1994;Milosavljević & Merritt 2001).
It is clear that in the cases mentioned above the following assumptions are fully justified: (i) the mass of the test particle is much larger than the mass of the field ones (m t ≫ m). This is a realistic situation, for example, when studying the sinking of GCs in galaxies (m GC ≃ 10 5 M ⊙ ≫ m * ≃ 1M ⊙ , where m GC and m * are, respectively, the typical mass of GCs and that of a star in a galaxy).
(ii) the field particles all have the same mass. This assumption becomes realistic if the previous condition is verified: when the test particle is much more massive than the field ones, the background can be safely approximated by stars with mass equal to their average value.
Clearly the case of BSSs in GCs is a significant exception to points (i) and (ii), being the mass of the test particle only slightly (2-3 times) larger than that of the field stars, so that taking into account a mass spectrum for the background can make significant differences. Moreover, real GCs are composed of stars in a relatively large range of masses (nominally from 0.1 up to 0.8 M ⊙ ), then also assumption (ii) is not strictly valid. Finally, it should be noticed that in general (e.g., Bertin 2001), the ratio between the DF and the two-body relaxation times is given by: t DF /t 2b ∝ 2m/(m t + m). Hence, while in the case of massive objects DF effects manifest on time-scales shorter than the two-body relaxation time of the system (t 2b ), for comparable masses these effects occur on quite long times, with the tendency to mix with two-body relaxation time effects. All these considerations make the problem of modeling the DF action in a GC more complex.
For the reasons above, it is not surprising that the case in which the field particles have a mass spectrum has not been extensively investigated in the literature. However, the considerations presented in Ciotti (2010) for a homogeneous and infinite density background with a mass spectrum, coupled with the well known dynamical evolution of a multi-mass GC (e.g., Spitzer 1987), prompted us to investigate in more detail the problem. In fact, the dynamical evolution of the parent GC leads to a radially dependent stratification of masses for the background stars, so that in practice each radius is characterized by a different mass function, leading to a radially dependent DF strength. Moreover, it should be noted that, if a sort of equipartition is established in the background particles, then some additional non trivial effect due to their velocity distributions may take place (Ciotti 2010, and Sect. 2). For example, in case of a bottom-heavy mass spectrum, on one side the low-mass population can have a higher density than that of high-mass stars (thus providing a proportionally larger contribution to the total DF), but, from the other side, the velocity of the light particles is also higher, so that their contribution to DF is reduced, in a compensating effect. Since DF depends on the "local" conditions of the system, in terms of both density and velocity distribution, it is not trivial to predict the final effect on the test particle due to the interplay between the various masscomponents in the cluster. Furthermore, the dynamical evolution of a GC produces significant changes on these distributions in time (e.g., during and after the core-collapse stage), which, in turn, depend on the component stellar mass, too. Hence the problem of modeling the DF on (the slightly heavier than average) BSSs in a dynamically evolving (multi-mass) GCs is a quite complex task. Therefore, what really happens to the BBS population in a GC can be analyzed only by considering the combined effects, at each radius, of both a radially and a time dependent mass spectrum obtained (for example) by N-body simulations, in a self-consistent way.
Within this context, we investigate here the multifaceted nature of DF of test particles slightly heavier than average field particles in a dynamically evolving system with a mass spectrum, combining the approach described in Ciotti (2010) with a set of numerical N -body simulations. The main goal is to verify whether in the presence of a multi-mass background, the DF time-scale develops a non-monotonic radial behavior in some stages of the system evolution, that could provide an additional explanation to the observed variety of BSS radial distributions.
In Section 2 we introduce the analytic approach to the problem of DF. In Section 3, we outline the mono-mass and multi-mass N -body simulations used to describe the background field component. The results are discussed in Section 4 and summarized in Section 5.

Mono-mass case
In order to introduce the case of a mass spectrum, we begin by recalling the relevant aspects of DF in the standard case of a background made by identical scattering masses m; this case is a useful benchmark which allows to better identify the mass-spectrum effects on DF. The deceleration of a test particle of mass m t , moving with speed v t and modulus v t in a homogeneous background of particles of equal mass m, constant number density n, and isotropic velocity distribution f (v), under the effect of DF can be written as where G is the gravitational constant, Ξ(v t ) is the fraction of particles slower than v t and ln Λ is the velocity weighted Coulomb logarithm (Chandrasekhar 1943(Chandrasekhar , 1960Binney & Tremaine 1987;Bertin 2001). By definition, Ξ(0) = 0 and Ξ(∞) = 1: we recall that the sharp truncation of the function Ξ for velocities larger than v t is an approximation (not affecting our discussion) due to the assumption of velocity isotropy of the background and to the use of the lowest order term in the impulsive approximation adopted to compute the two-body interactions. It is well known that several problems affect the direct application of eq. (1) to spherical systems, as a local description of DF (Bontekoe & van Albada 1987; Bertin et al. 2003Bertin et al. , 2004Arena et al. 2006). Nonetheless several studies have been based on the applications of eq. (1) to spherical systems.
Here we follow the same approach, for the case of interest (i.e., that of the evolution of the BSS population of a GC), and take into account the radial dependence of the number density by replacing n with n(r). In addition, in all our discussion we assume the test particle to be on circular orbit in the GC potential well, so that also v t depends on r: where M (r) is the total mass of the system enclosed within a sphere of radius r. From these assumptions, it follows that also the function Ξ depends only on the radial distance from the center through v t (r) and the local velocity distribution of the field particles. As usual, from eq. (1) we can define the characteristic DF times-scale as so that in our case where ρ(r) = n(r)m is the mass density of the background. This expression explicitly shows that the radial dependence of the DF time-scale is shaped by three functions: the velocity of the test particle, v t (r), the mass density of the field stars, ρ(r), and the relative number of stars moving slower than v t at any radius, Ξ(r). For simplicity, we neglect the possible dependence of Λ on radius: at the present level of approximation this seems a reasonable assumption, due to the logarithmic nature of the associated term. Note that the function Ξ(r) is sometimes evaluated analytically by assuming a local Maxwellian velocity distribution for the background particles, determined by the value of the local velocity dispersion (e.g., Binney & Tremaine 1987). Here however we avoid this additional assumption, as we compute the function Ξ(r) directly by counting the number of particles moving slower than v t in the N-body simulation outputs (see Section 3).

Multi-mass case
As described in the Introduction, this study is focused on the case of a system made by the superposition of different mass components, such as a GC with stars distributed according to a prescribed mass spectrum. Ciotti (2010) showed that in such circumstances the DF experienced by a test particle can be significantly different from that experienced in a single-mass background. Therefore, the natural question arises of what happens in a GC, where the mass spectrum is associated to the initial mass function (IMF), and the dynamical evolution of the GC leads to a redistribution (through the tendency to equipartition, e.g. Spitzer 1987) of the density and velocity profiles of the different mass components. It is easy to realize that all these trends, weighted by the ratio between the mass of a BSS and the mass of the field stars in each subcomponent of the cluster, could lead to a quite complicate radial trend of the total DF.
In the presence of field particles with a mass spectrum, the total DF deceleration can be split in the individual contributions due to each single population. The i th population causes a deceleration of the test particle which is again expressed by eq. (1); where now ρ i (r) = n i (r)m i is the local density of the background component with stellar mass m i . Due to the additive nature of scattering effects in the Chandrasekhar treatment of DF, the total deceleration is obtained from the sum of all the contributions. In particular, from eqs. (3) and (5) it follows that: where In eq. (6) we assumed that the mass spectrum of the background stars is represented by the sum of a finite number of components (N pop ). However eq. (6) can be easily written as an integral in the case of a continuous mass spectrum (Ciotti 2010). As a general comment, note that eq. (6) indicates that the local value of the DF time is roughly determined by the smallest among the various t i DF (r). So far, the description in this Section just reflects the "standard" approach to our problem: namely, one could construct an equilibrium model for a GC (for example by solving numerically the Poisson equation for one-and multi-component King (1966) or Wilson (1975) models), take the radial profile of n i (r) and the local velocity distribution and compute eq. (7). Here, we follow a more realistic approach, namely we make use of (collisional) N-body simulations of a set of mono-and multi-mass GC models, capable to provide us with the self-consistent radial behaviors of all the various quantities needed to evaluate eq. (7). This approach, at variance with the solution of the Poisson problem, allows us to take into account also the time evolution of the t DF radial profile.

N -body simulations
Here we describe the set-up of the N-body simulations performed in order to determine the time evolution of the GC models hosting the BSS population. Since GCs are collisional systems, very accurate and specifically designed numerical methods are required to properly describe their time evolution. In practice, we use N-body simulations to obtain a self-consistent description of the phase-space density distribution of the different components of the GC, and then we apply the equations in Section 2 to estimate the radial trend of t DF . We also perform some simulation where the background stars are all characterized by the same mass, so that the effects of a mass spectrum can be better appreciated through comparison. For our simulations we used the direct N -body code NBODY6 (Aarseth 2003). In all cases, the simulations are not meant to describe the evolution of a GC from its formation to the present days. They just provide a simplified "picture" of a current GC, which could be in a pre-or in a post-core collapse phase.

Mono-mass simulations
The mono-mass system is composed of N = 10 4 particles with mass m. The initial conditions (particle positions and velocities) have been generated from a King (1966) model with central dimensionless potential W 0 = 4. We followed the dynamical evolution of the cluster up to a final time t f = 1000 in N-body time units (i.e. in units in which the total mass of the cluster is M = 1, G = 1 and the total energy is E = −1/4; see Aarseth 2003, Sect. 7.4). A characteristic time scale for two-body relaxation is given by t rh (0), the half-mass relaxation time of the initial conditions (e.g., Spitzer 1987): where r h (0) is the half mass radius (i.e., the radius containing half the total mass of the system) at the initial time t = 0. In N-body units t rh (0) ≃ 210 so that in practice we follow the GC evolution up to t f ≃ 4.8t rh (0). To improve the statistics, 20 different sets of initial conditions have been generated by changing only the random seed from which the positions and velocities of the particles are extracted starting from the distribution function. All the realizations have then been combined at each time-step, shifting the center of mass of each system to a common origin. This procedure generated "supersnapshots" containing N super = 2 × 10 5 particles. The simulations have been run on a dedicated workstation, and each simulation in the mono-mass case took approximately 3 hours (thanks to the use of a GPU card).

Multi-mass simulations
As described in Section 2.2, in order to determine the time evolution of the radial trend of t DF for the population of BSS in a multi-mass GC, we need the distribution function of each stellar component of the cluster, i.e., the associated density profile and its velocity distribution. Of course, some educated guess could be used to describe the radial profile of the velocity distribution (e.g., to solve the associated Jeans equations and use the resulting velocity dispersion in the local Maxwellian approximation), but here we prefer to use direct N-body simulations that allow to compute in a self-consistent way the evolution of the structural and dynamical properties of the different components.
The modelization of a mass-spectrum with a necessarily limited number of particles imposes some constrain in the choice of the number of mass bins. Of course, the larger is the number of bins, the finer is the spectrum. However, to avoid too little numbers of particles in each bin as a consequence of an excessive partition among the field particles, as well as to understand more clearly the contribution of the various mass ranges (at the various radii) to the resulting DF, we represent the mass spectrum as the superposition of three different populations with a total number of stars N 1 , N 2 , and N 3 . The masses of the individual stars in each population are m 1 , m 2 , and m 3 respectively, with m 1 being the smallest value, m 2 = 2m 1 and m 3 = 3m 1 . Therefore, N = N 1 + N 2 + N 3 and m 1 = M/(N 1 + 2N 2 + 3N 3 ), where M is the total mass of the cluster. The three populations are aimed at grossly representing the main populations of field stars in a present-day GC, namely low-MS stars (with mass m 1 ≃ 0.3M ⊙ ), intermediate-MS objects (m 2 ≃ 0.6M ⊙ ), and TO and giant stars (m 3 ≃ 0.9M ⊙ ). As in the mono-mass simulations, the total initial number of stars is N = 10 4 , with the three components counting N 1 = 8500, N 2 = 1200 and N 3 = 300 particles, respectively. This choice is a reasonable compromise between the description of a realistic case, and the need of a large statistical sampling for the less numerous population (m 3 ) in order to avoid noise-dominated results (see Sect. 5 for a discussion of the impact of such approximation on the results obtained).
The initial positions and velocities of all particles in each of the three groups are randomly extracted from the same King (1966) distribution function, i.e., no initial mass segregation is assumed. This choice is made both to avoid an additional degree of freedom (an unconstrained amount of mass segregation) and because it is justified by the flat BSS distribution observed in ω Cen, Palomar 14 and NGC 2419 (see Ferraro et al. 2012). The system is also fully isolated, with no primordial populations of binaries or multiple systems (see Sect. 5 for the discussion of these approximations).
In order to improve the statistics, as in the mono-mass case, we generated 20 sets of initial conditions by varying only the random seed and we combined all the 19 runs together at every extracted snapshot, thus generating supersnapshots made of N super ≃ 2 × 10 5 stars (this number varies slightly during the evolution due to the loss of unbound particles). All the simulations have been stopped at the time t f ≃ 10t rh (0), with t rh (0) = 190 in N -body units. In these multi-mass cases, the half-mass relaxation time is again computed by using eq. (8). We extracted a snapshot every 5 × 10 −3 t rh (0), thus we have guaranteed a good accuracy in tracking the cluster evolution. Three different values of the King dimensionless potential (W 0 = 4, 6, 8) have been considered for the multi-mass simulations, so in total we ran 60 simulations in the multi-mass case.

Mono-mass case
In the mono-mass case, the DF time-scale of a test particle of mass m t = 4m (m being the mass of the background stars) has been evaluated by using eq. (4). In order to extract the radial profiles of v t (r), n(r) and Ξ(r) entering eq. (4) from the simulations, we considered equally populated radial bins (i.e., 100 concentric spherical shells, each enclosing 2000 particles of a given supersnapshot). This choice has the useful property of maintaining constant the error bars of the number density over the whole radial range, thus reducing the effects of random fluctuations. This choice also implies a finer sampling of the innermost regions as the time passes, because the density progressively increases in the core of the GC. The number density profile of the field component, n(r), is then given by the number of particles in each radial bin divided by the volume of the shell. The circular velocity of the test particle, v t (r), is trivially computed by using eq. (2). Finally, the local estimate of Ξ(r) in each radial bin is obtained by normalizing the number of background particles in the given radial bin that are slower than v t (r), to the total number of background particle in the same bin.
The resulting radial trend of t DF at four representative times is shown in Figure 1, where the radial distance is expressed in units of the half-mass radius at that time r h (t), and the time is normalized to the instantaneous half-mass relaxation time, t rh (t). Figure 1 clearly shows that t DF maintains a monotonic radial trend at increasing time, in agreement with simple expectations. We see also that, as the dynamical evolution of the system proceeds, t DF decreases (i.e., DF becomes more efficient) in the central and external regions, while the effect is opposite at intermediate radii. This behavior is mainly due to the time evolution of the mass density (see Figure 2), which progressively increases with time in the innermost and outermost regions, (the former being due to the core-contraction, the latter due to the related cluster halo expansion), while it tends to decrease on the radial interval −0.3 log(r/r h ) 0.3 (consistently with an evolution towards a higher cluster concentration, resulting in a decrease of r h ). Indeed the density profile is the primary driver of the shape of t DF (r) at all times, while the other terms entering eq. (4) provide a negligible contribution. This is due to a more significant evolution of ρ(r) = m n(r) with respect to those of Ξ(r) and v 3 t (r), as is clearly apparent in Figure 3. Therefore, the mono-mass simulations confirm the scenario adopted by Ferraro et al. (2012) to explain the different radial distributions of BSSs as a function of the cluster age, due to the monotonic increase of the radius at which the cluster age coincides with t DF (r).

Multi-mass case
In the multi-mass case, the DF time-scale has been computed by evaluating each term in eq. (6) associated to each one of the three stellar components considered in the simulations described in Sect. 3.2, and by setting m t = 4m 1 . To avoid large fluctuations due to the low number statistics of the most massive component (which counts only 3% of the total number of particles), the radial binning has been chosen by imposing that at least ≃ 200 stars of mass m 3 are included in each shell and that the total number of bins is always N bin ≥ 19. Thus, at odds with the mono-mass case, the total number of particles in each shell is not constant. The number density profile of each mass component, n i (r), has been computed as in the mono-mass case (see Sect. 4.1), while the circular velocity of the test particle, is again trivially computed from eq. (2). The velocity factors Ξ i (r) in each radial bin have been computed by normalizing the number of stars in the i th mass group that move slower than v t (r), to the total number of particles of the same mass group in that radial bin.
The results obtained for the three considered values of the dimensionless potential (W 0 = 4, 6, 8) and four different times are plotted in Figure 4. In particular, we show t DF (r) at the initial time of the simulation (t = 0) and for three snapshots around the "core-collapse time" t cc , defined as the time at which the Lagrangian radius r 10 (i.e., the radius containing 10% of the total mass) reaches its minimum value. For reference, we notice that t cc ≃ 3.6t rh (0), 2.5t rh (0), 0.6t rh (0) for W 0 = 4, 6, 8, respectively. The radial behavior of the DF time-scale is plotted both for each mass-component separately (color lines), and combining the effects of the three background mass components according to eq. (6) (thick black lines). Since the initial conditions for all mass-groups were built from the same King model properly scaled only for the adopted number of stars, at t = 0 the radial profile of t DF has the same qualitative behavior for any group and its value at any radial distance decreases for decreasing particle mass, because lighter stars are by far the most numerous and thus dominate the local mass density.
At all evolutionary times, as for the mono-mass simulations, also in this case the mass density is the main driver in shaping t DF (r): the local value of t DF (r) is essentially determined by that of the dominant contributor to ρ(r). In fact, due to energy equipartition, the most massive particles progressively migrate toward the center, while the lightest component expands outward. As a consequence the total mass density is mainly contributed by the heaviest mass-group at small radii, while it is dominated by the m 1 particles in the outskirts, and the radial profile of t DF is thus driven by 1/ρ 3 (r) and 1/ρ 1 (r) in the two respective radial regions. In addition, from t = 0 to t = t cc , DF becomes increasingly more efficient in the center (i.e., the value of t DF at fixed small radii becomes smaller) because also the central mass density increases, due to the segregation first of m 3 particles and then also of m 2 stars 2 . The increase of t DF (r) at t = 2.5t cc , is due to a late re-expansion of the system (corresponding to a decrease of the mass density).
Figure 4 also shows that the radial profile of t DF is qualitatively the same for all values of W 0 (i.e., independently of the cluster concentration) and it always monotonically increases with radius, as in the mono-mass case. Our simple analysis then shows that the presence of a mass spectrum does not induce non-monotonic behaviors in the DF time-scale.

Equivalent Classical System
The previous analysis showed that a mass spectrum in a GC does not leads to a non-monotonic behavior of t DF (r), so that the results are in qualitative agreement with the mono-mass case. But what about the absolute value of t DF ? In fact, Ciotti (2010) showed that a significant overestimate of the DF time-scale can arise if approximating a multi-mass system with a single component only. In our case this check can be done by introducing the definition of the Equivalent Classical System (ECS): • the number density in the ECS is equal to the total number density of the multi-mass case: • the mass of the field particles in the ECS (m ECS ) is equal to the average field mass of the multicomponent case: where N is the total number of stars in the multi-mass system; • the velocity dispersion of the ECS is equal to the equipartition velocity of the multi-component case.
The last property is non relevant in our case, since we compute the function Ξ(r) directly from the simulation outputs, as the fraction of all the stars slower than the test particle. In order to compare the DF time-scale obtained in the ECS approximation with the exact determination for the multi-mass case discussed in the previous section, we used eq. (4) by assuming n(r) = n ECS (r), m t = 4 m 1 and m = m ECS ≃ 1.2 m 1 .
To compute M (r) entering the definition of v t , we counted the number of particles with the radius r and multiplied it by m ECS . The comparison between the DF time-scale in the multi-mass and in the ECS cases for W 0 = 8 is shown in Figure 5 (solid lines). As expected, no differences are found at t = 0, when masssegregation has not played any role yet. This is due to the fact that the mass density is nearly the same in the two cases, because the three components have equal radial distributions and m ECS well corresponds to the true mean mass at all radii. When mass-segregation takes place, a systematic underestimate of DF in the ECS case becomes apparent in the central regions, where t DF can be almost a factor of 3 larger than in the true multi-mass system. This discrepancy is due to the systematic underestimate of the ECS local mass density with respect to the multi-mass one. In fact, as a consequence of mass segregation, the number density at small radii is mainly contributed by the heaviest stars, and since m 3 is significantly larger than m ECS , the ECS mass density in these regions is considerably smaller than the multi-mass one. A smaller density implies a lower efficiency of DF, thus bringing to larger values of t DF in the ECS with respect to the multi-mass case. Beyond the half-mass radius, the ratio between the values of t DF in the two cases becomes nearly flat around 1. This is because beyond the half mass radius no appreciable mass-segregation has taken place during the simulated time, and so the local average mass is nearly equal to m ECS in that region, as in the initial conditions.
As a final test we modified the definition of ECS in requiring that the mass density in this approximation is equal to the true mass density of the multi-component system at any radius: As expected the difference between the two representations is considerably reduced and the underestimate of t DF is about half that previously found (see the dashed lines in Figure 5). Therefore, we conclude that special attention should be paid to arguments based on specific requests about the value of t DF , as estimate based on the average properties of the GC can be off by a factor of 2 or 3, factors that are quite important for phenomena happening on time scales of the order of the Hubble time.

Summary and discussion
We studied the radial behavior of the DF time-scale computed by following the Chandrasekhar's formalism in a background of field particles with a "simplified" mass spectrum, evolved by means of direct and collisional N-body simulations. To our knowledge, this is the first numerical study of DF in the presence of a self-consistently evolving field of particles with unequal masses and for a test star with a mass comparable to that of the field (see Ciotti 2010 for an analytical approach to the problem in the case of a spatially homogeneous background). We explored both the mono-mass case (where all particles have unit mass and the test star decelerated by DF is four times more massive), and the multi-mass case (with three field components of masses m 2 = 2 m 1 , m 3 = 3m 1 and total numbers N 1 , N 2 , N 3 , and with the test particle having m t = 4m 1 ). Each simulation was run with N = 10 4 particles initially distributed as a spherical and isotropic King model with dimensionless potential W 0 = 4 in the mono-mass case, and W 0 = 4, 6, 8 in the multi-mass one. To improve the statistics, we combined 20 simulations differing only for the seed used to generate the initial conditions. We find that the radial behavior of t DF is always monotonic (t DF increasing with r), both in the monomass background and in the presence of a mass spectrum, independently of the evolutionary time and the initial concentration of the system. In all cases, the largely dominant factor in determining the shape of t DF (r) is the total mass density profile.
We also find that approximating a multi-mass system as single-mass cluster (made of stars with masses equal to the average particle mass) can lead to a systematic overestimate of t DF within the half-mass radius, up to a factor of 3 in the innermost regions, so that some care must be used when using the average properties of background populations to obtain quantitative estimates of t DF .
The N −body simulations and the overall approach presented in this work are certainly a rough simplification of the much more complex problem of DF in real GCs. First of all, the total number of simulated particles is much smaller than the number of stars in a cluster. However, while this affects the overall cluster evolution time-scale, it is not expected to impact our conclusions about the radial monotonicity of t DF . The same considerations apply to the assumptions adopted for the mass spectrum (only three bins with m 2 = 2m 1 and m 3 = 3m 1 ): using a larger number of mass groups is not expected to induce a non-monotonic behavior on t DF (r), since the density profile ρ(r) would still be a monotonic decreasing function of radius. Adding a population of black holes, neutron stars and massive white dwarfs (i.e., stellar remnants more massive than BSSs, which are certainly present in current GCs) could make some difference, especially in the cluster central regions where they are expected to be concentrated because of mass segregation. However, given the uncertainties on their typical retention fraction, we preferred to neglect such a component in this first approach to the problem, postponing this issue to a future paper. In any case, still because t DF ∝ 1/ρ(r), we do not expect that including massive dark remnants would modify our conclusions substantially. This also holds for the effect of an external tidal field, which is mainly expected to favor the evaporation of the lightest stars from the cluster outskirts. The impact of a population of primordial binaries and multiple systems, instead, is more difficult to predict (especially if also BSSs are modeled as binary systems) and we plan to explore this problem in a forthcoming paper. Fig. 1.-Time evolution of the radial profile of the DF time-scale, t DF , for the mono-mass system with central dimensionless potential W 0 = 4. The four curves correspond to four different times of the simulation, from t = 0 to the final time t f = 4.8t rh (0). As explained in the text, the DF time-scale is normalized to the instantaneous half-mass relaxation time, while the radius is in units of the instantaneous half-mass radius r h . The increase of the probed radial range at large radii with time is due to the decrease of r h , while in the central regions the different extension is due to the fixed number of stars assumed to define the radial bins (see the text).   Figures 1 and 2. In the plots we also show the three radially dependent terms contributing to the final value of t DF (v 3 t : dotted line, 1/ρ: dashed line, 1/Ξ: thin solid line). All quantities are given here in code units, and from eq. (4) it follows that log(t DF ) is just given (modulo an additive constant) by the sum of the three quantities. Fig. 4.-Radial profile of t DF in the multi-mass simulations, plotted for each mass-component separately (red: m 1 , blue: m 2 , green: m 3 ) and from the combined effects of the three background mass components according to eq. (6) (black lines). The different panels are snapshots taken at t = 0, 0.5, 1, and 2.5t cc , where t cc is the fiducial core-collapse time of the system (see Section 4.2). Panels from left to right, models with: W 0 = 4, 6, 8. The radial distance from the center and t DF are normalized, respectively, to the instantaneous half-mass radius and the instantaneous half-mass relaxation time computed for the system as a whole. . -Ratio between the DF time-scale computed in the ECS approximation and for the multi-mass system, as a function of radius (in units of the total half-mass radius), for three evolutionary times (see labels) and for the model with initial conditions W 0 = 8. The thick solid lines show the results obtained under the assumption that the number density of the ECS system equals that of the multi-mass case: n ECS (r) = n(r). The dashed lines refer to the approximation in which the mass density of the ECS system is equal to that of the multi-mass case: ρ ECS (r) = ρ(r). The results are essentially the same W 0 = 4 and W 0 = 6. Notice how the number average estimate can be off by a factor up to 3 in the central regions, while a much better estimate is obtained by using the mass density average.