Multiple populations in globular clusters: the distinct kinematic imprints of different formation scenarios

Abstract

1 INTRODUCTION

Once considered as quintessential simple stellar populations (i.e. no dispersion in age or metal content), globular clusters (GCs) now appear to host multiple stellar populations. This has been revealed by ubiquitous star-to-star abundance variations of light elements such as Na, O, and Al (e.g. Gratton, Carretta & Bragaglia 2012), indicating that a significant fraction of cluster stars (∼30–70 per cent; e.g. Carretta et al. 2009a) must have been polluted by material processed in high-temperature stellar interiors via proton capture reactions. It is also apparent from the multiple or extended main sequences, subgiant and red giant branches (RGBs) seen in the colour–magnitude diagrams of a growing number of Galactic GCs (e.g. Piotto et al. 2012), which have been attributed to different levels of helium enrichment (e.g. D'Antona et al. 2002; Norris 2004; Pasquini et al. 2011) and light element abundance variations (e.g. Marino et al. 2008; Sbordone et al. 2011) among the cluster stars.

Several models have been proposed in an attempt to explain the observed chemical abundance anomalies, with different types of stars suggested as the source of helium-enriched material. The yields from these potential enrichment sources must also explain the observed chemical abundance patterns (e.g. the Na–O anticorrelation) without producing a spread in iron abundances (such a spread is absent in almost all but the most massive GCs; e.g. Carretta et al. 2009b). Asymptotic Giant Branch (AGB) stars (e.g. Ventura et al. 2001; D'Ercole et al. 2008), fast rotating massive stars (‘spin-stars’; e.g. Decressin et al. 2007), interacting massive binaries (de Mink et al. 2009) and super-massive stars with M ∼ 10⁴ M_⊙ (Denissenkov & Hartwick 2014) have all been considered as possible polluters.

Obtaining the right chemistry for the polluting material is one aspect of the problem, but finding a way to get this material into stars is another equally important and challenging issue. Many of the models put forward invoke the formation of a second generation of stars from the processed material released by a first generation. For this reason, the terms multiple generations and multiple populations are often interchangeably used in the literature to refer to the presence of polluted and pristine stars, but note that the notion of multiple star formation events within a single cluster is an assumption of these models for which there is no direct observational evidence.

In a scenario that has received significant attention, D'Ercole et al. (2008) explored the possibility that a second generation is formed from the gas ejected by AGB stars. They used 1D hydrodynamical simulations to show that the low-velocity ejecta of AGB stars can be retained by the cluster and form a cooling flow that rapidly collects in the innermost regions. The new generation then forms a more concentrated stellar subsystem, in keeping with observations of a more centrally concentrated polluted (i.e. Na-rich, O-poor) population in a number of Galactic GCs (e.g. Lardo et al. 2011). This scenario was further studied by different authors (e.g. Bekki 2011; Conroy & Spergel 2011), many of which stressed the importance of accreting pristine interstellar gas on to the cluster for it to work. In particular, it was suggested that the total mass of accreted pristine gas must be comparable to that of AGB ejecta in order to reproduce the observed Na–O anticorrelation (e.g. D'Ercole et al. 2010; Ventura et al. 2013). How such pristine interstellar material can remain free of pollution by supernova ejecta from the first generation (as required by the lack of an iron abundance spread) however remains unclear. In addition to the above model based on AGB stars, the main scenario invoking spin-stars as the source of pollution also works under the assumption that the slow ejecta from these stars form the basis of a new stellar generation (e.g. Decressin et al. 2007). Charbonnel et al. (2014) even postulated that, within this framework, all low-mass stars present in GCs today could in fact be second-generation stars.

The models that require multiple star formation events come with many potential shortcomings (for an overview see Renzini 2008; Bastian et al. 2013a; Cabrera-Ziri et al. 2015). For example, they generally suffer from an ‘internal mass budget’ problem, as the amount of processed material returned by either AGB stars or spin-stars is insufficient to account for the large observed fractions of polluted stars in GCs. To circumvent this issue, these models either invoke a non-standard initial mass function (IMF) for the first or second generation (e.g. Decressin et al. 2007) and a 100 per cent star formation efficiency for the second-generation stars (e.g. D'Antona et al. 2013), or imply that GCs were initially ∼10–100 times more massive than observed at the present time and that a large fraction of first-generation (i.e. pristine) stars have since been lost (e.g. D'Ercole et al. 2008; Bekki 2011; Schaerer & Charbonnel 2011). The latter possibility appears to be in conflict with the ‘external mass budget’ constraints imposed by observations of GCs and field stars in nearby dwarf galaxies (Fornax, WLM, IKN), which imply that GCs in these systems could not have been more than about five times more massive initially (Larsen, Strader & Brodie 2012; Larsen et al. 2014). These models are also at odds with observations of young massive clusters forming at the present epoch, for which there is currently no conclusive evidence for age spreads, ongoing star formation, or reservoirs of gas and dust from which a new generation could form (e.g. Kudryavtseva et al. 2012; Bastian & Silva-Villa 2013; Bastian et al. 2013b; Bastian & Strader 2014; Bastian, Hollyhead & Cabrera-Ziri 2014; Cabrera-Ziri et al. 2014).

To address these problems, Bastian et al. (2013a) proposed a new scenario in which processed material shed from interacting massive binaries during the first ∼10 Myr of the cluster's life is accreted on to the circumstellar discs of low-mass stars (still on the pre-main sequence), ultimately polluting those stars. It must be said that the suggestion that pollution could originate from the accretion of ejecta from higher mass stars was first made by D'Antona, Gratton & Chieffi (1983), who at the time considered that the observed abundance patterns were reflecting surface contamination. More recent evidence however indicates that abundance anomalies are present throughout the star, as evolved RGB stars (mixed through convection) show the same Na–O anticorrelations as main-sequence stars (Cohen, Briley & Stetson 2002). The variation proposed by Bastian et al. (2013a) addresses this by polluting stars when they are on the pre-main sequence and convective, and it also maximizes the efficiency of accretion by invoking gas capture by discs.

In the so-called ‘early disc accretion’ scenario, the polluted stars originate from the same generation as the pristine stars, with no need for multiple star formation events. Only the stars that pass through the cluster core to sweep up gas are polluted, as this is where high-mass stars and their ejecta are expected to be preferentially located, either because they form there (e.g. Bonnell, Bate & Zinnecker 1998) or because they segregate soon after their formation (e.g. Portegies Zwart, McMillan & Gieles 2010). Depending slightly on the velocity structure and the gas distribution, this leaves about half of the low-mass stars, which did not cross the core, with normal abundances. The model also predicts an enriched population that is more centrally concentrated than the pristine population. Given that a large fraction of high-mass stars are in binary systems that will interact at some point in their evolution (Sana et al. 2012, 2013), these interacting systems are promising sources of polluting material. It has been argued that their ejecta can display the required chemical properties and also satisfy the internal mass budget constraints (de Mink et al. 2009) without requiring GCs to have been significantly more massive initially.

Among the caveats of the early disc accretion model, an important one is the uncertainty on the expected lifetime of circumstellar discs around low-mass pre-main sequence stars in the dense environment of a GC. These discs need to survive for 5–10 Myr to sweep up enough material for this scenario to work, and Bondi–Hoyle accretion is expected to be inefficient due to the large velocity dispersions of massive clusters. D'Antona et al. (2014) also pointed out that the structure of the seed stars may not remain fully convective for the whole duration of the accretion phase. If this is true, the full mixing of the helium-enriched material required down to the stellar centre would not be achieved. Thermohaline mixing may offer a way around this, but the same authors suggested that this could destroy any lithium surviving in the envelope, which may be in contradiction with observations (e.g. Shen et al. 2010).

Both the better-studied multiple generations model (with AGB stars as the source of pollution; e.g. D'Ercole et al. 2008) and the newer early disc accretion model (Bastian et al. 2013a) have their difficulties, but neither has been conclusively ruled out to date. In this paper, we therefore compare two main scenarios for how enriched material makes its way into stars: one in which a second generation is formed from gas that collects at the centre of the cluster (hereafter referred to as the multiple generations scenario), and another one in which pollution is due to accretion on to low-mass stars from the same generation when they pass through the core of the cluster (hereafter referred to as the accretion scenario). Instead of focusing on the chemistry as is often the case in studies of multiple populations, we approach this problem from another angle and explore the fossil kinematic imprints that may follow from these two formation scenarios after a Hubble time of dynamical evolution. Our results are not tied to a specific source of polluting material. They apply to any ‘multiple generations’ scenario that proceeds as suggested by D'Ercole et al. (2008), whether AGB stars are the polluters or not. Similarly, they would be relevant for any ‘accretion’ scenario that proceeds as outlined by Bastian et al. (2013a), regardless of the source of pollution or whether accretion proceeds via discs or not.

In Section 2, we start by reviewing previous work on the dynamics of multiple populations, which was mainly set in the context of models invoking multiple generations. Section 3 describes the N-body simulations that we performed and the different initial conditions that we chose to capture the distinct kinematic properties implied by the two self-enrichment models considered. We then present in Section 4 the results of these simulations, with an emphasis on the identification of a unique kinematic signature that allows us to distinguish between the two models. In Section 5, we discuss how this kinematic signature could be observed and we identify a large subsample of Galactic GCs in which it is expected to have survived to the present epoch. We present our conclusions in Section 6.

2 PREVIOUS STUDIES OF THE DYNAMICS OF MULTIPLE POPULATIONS

As mentioned previously, the internal mass budget problem of models invoking multiple generations can in principle be avoided by assuming that GCs were initially much more massive. A significant fraction of first-generation stars (∼90 per cent or more) must then be lost (i.e. unbound) by dynamical evolution, while the more concentrated second generation remains relatively unaffected. Decressin, Baumgardt & Kroupa (2008) showed with N-body simulations that two-body relaxation alone cannot cause such a strong preferential loss of first-generation stars, and suggested that a more efficient mechanism such as primordial gas expulsion is needed to unbind the first-generation stars located in the outer regions of the cluster on a short time-scale.

Decressin et al. (2010) found that, under favourable circumstances, the expulsion of gas left over from star formation can lead to a cluster containing 60 per cent of second-generation stars. The radial distribution of the two populations is still expected to be different after gas expulsion, so two-body relaxation could further increase the fraction of second-generation stars. Note however that the importance of early gas expulsion in young massive clusters has recently been questioned by hydrodynamical simulations (Kruijssen et al. 2012) and kinematic observations (e.g. Cottaar et al. 2012; Hénault-Brunet et al. 2012b) suggesting that newly formed clusters are relatively gas-poor at their centre and in virial equilibrium from a young age (a few Myr).

D'Ercole et al. (2008) studied the effect of mass-loss from Type II supernova ejecta (without including primordial gas left over from star formation) and suggested that this alone could trigger rapid early expansion of the cluster and lead to a strong preferential loss of first-generation stars. In any case, it seems that the initial conditions and mass lost in this early phase need to be fine-tuned to explain that the mass ratio of the observed populations is always of order unity. One may otherwise expect GCs experiencing weak tidal fields to contain a much larger fraction of pristine stars. This does not seem to be the case for the remote and massive cluster NGC 2419 (R_G = 94.7 kpc, M_V = −9.5). Multiband photometric observations of this cluster indeed suggest that a significant fraction of its stars belong to a helium-enriched population (Beccari et al. 2013).

In the following sections, we will nevertheless ignore these caveats about early mass-loss and follow the long-term dynamical evolution of multiple populations starting from a cluster in virial equilibrium with a similar number of polluted and pristine stars. This approach was also adopted by Vesperini et al. (2013), who studied the long-term evolution of the relative spatial distribution of first- and second-generation stars (with spherically symmetric distributions) in the context of the scenario presented by D'Ercole et al. (2008). Starting with a second generation that is more centrally concentrated than the first generation, they explored how the time-scale for spatial mixing of the two populations depends on the initial concentration of the second generation. They found that complete spatial mixing is expected only for dynamically evolved clusters which have lost at least 60–70 per cent of their mass due to two-body relaxation (i.e. ignoring the rapid early loss of first-generation stars), irrespective of the initial concentration of the second generation. In Section 5, we comment on the fraction of Galactic GCs expected to have evolved to that stage.

In addition to having different spatial distributions, stars belonging to different populations could also display distinct kinematic properties. Bekki (2010, 2011) simulated the formation of second-generation stars from AGB ejecta within a cluster and showed that if the first generation formed with some small net angular momentum, the second generation can show considerable rotation and even adopt a flattened distribution initially. This is a simple consequence of conservation of angular momentum, along with dissipative processes driving the polluted material to higher densities towards the centre of the cluster. Such properties are reminiscent of the multiple stellar populations (with different ages and chemical properties) observed in nuclear star clusters, the youngest of which often form centrally concentrated stellar discs (e.g. Seth et al. 2006). Bekki (2010, 2011) also showed that the second-generation stars are expected to initially have a lower velocity dispersion than the first-generation stars.

Mastrobuono-Battisti & Perets (2013) were the first to study the long-term dynamical evolution of a GC with nested structures while also taking into account rotation and a flattened spatial distribution for the more concentrated population. Their N-body simulations, tailored to ω Cen, trace the evolution of a relatively low-mass (10⁵ M_⊙) cold disc component embedded in a massive GC (2.5 × 10⁶ M_⊙). They showed that these initial conditions can leave behind signatures even after 12 Gyr of relaxation-driven dynamical evolution.¹ For example, while the disc stars do become more isotropically distributed with time, they do not attain a completely relaxed spherical shape at the end of the simulations. Also, as these disc stars exchange angular momentum with the other stars, the whole cluster can become slightly flattened. Another prediction of this study is a lower velocity dispersion and a more radially anisotropic velocity distribution for the polluted (disc) stars compared to the pristine stars. We will address the uniqueness of these kinematic signatures in Section 4 when comparing the long-term evolution of the accretion and multiple generations scenarios.

Before moving to our own N-body simulations, we note that the work presented here is also motivated by recent observational results. For example, from Hubble Space Telescope (HST) proper motion measurements in 47 Tuc, Richer et al. (2013) found that the stars on the bluer side of the main sequence (presumably more helium-rich/polluted) exhibit the largest proper-motion anisotropy, with preferentially radial orbits, while the presumably pristine stars on the redder side of the main sequence display no velocity anisotropy (in the plane of the sky). The larger radial anisotropy of the polluted stars was interpreted by these authors as a possible effect of two-body relaxation, with the more centrally concentrated polluted population slowly diffusing outward. Kinematic differences were also observed between the stellar populations of 47 Tuc from radial velocity measurements. Kucinskas, Dobrovolskas & Bonifacio (2014) indeed reported that the line-of-sight velocity dispersion of their subsample of more enriched stars is lower than that of their subsample of presumably pristine stars by about 1 km s⁻¹.

Based on radial velocities and Na abundance measurements in 20 Galactic GCs, the study of Bellazzini et al. (2012) is the only one to date that has explored the possible connection between kinematics and multiple chemical populations in a large sample of clusters. Their data generally did not reveal an obvious correlation between either velocity dispersion or systemic cluster rotation and Na abundance (a convenient tracer of self-enrichment). Three of their clusters (NGC 6388, NGC 6441, NGC 2808) show tentative evidence for a drop in the velocity dispersion with increasing Na abundance, but this will need to be verified with larger and more evenly distributed samples. While it was reported that the Na-rich and Na-poor subsamples display similar rotational patterns in NGC 6441 and NGC 6388, there is some indication that Na-poor stars display a larger rotational amplitude than Na-rich stars in NGC 2808 (≃ 2.25 kms⁻¹ and 1.0 kms⁻¹, respectively²). Two other clusters (NGC 6171 and NGC 7078) show hints of Na-poor stars having a larger rotational amplitude by a few kms⁻¹, albeit from small samples.

3 N-BODY SIMULATIONS

3.1 Initial conditions for the pristine and polluted populations

We studied the dynamical evolution of polluted and pristine stars in the context of the accretion (Section 3.1.1) and multiple generations (Section 3.1.2) scenarios by means of N-body simulations. More details about the aspects common to all of our simulations are summarized in the next subsection. Here, we start by describing the initial conditions of the different populations and how these depend on the pollution scenario considered.

3.1.1 Accretion scenario

To set up our systems, positions and velocities of the stars were sampled randomly from an isotropic isochrone model, and the masses of the particles were sampled from a Kroupa mass function (Kroupa 2001, see also Section 3.2), with no correlation between stellar mass and position.

We assumed that the enriched material – expelled by massive stars in the model proposed by Bastian et al. (2013a) – is located within the core radius of the cluster and we split the stars into two categories: stars that do not enter the core of the cluster are considered as pristine stars, while stars on an orbit with a pericentre smaller than the core radius are assumed to accrete and are flagged as members of the polluted population. This yields a cluster containing about 45 per cent of polluted stars at the start of the simulation, and these are by construction more centrally concentrated than the pristine stars.

Note that both in Bastian et al. (2013a) and in this work, the fraction of polluted stars is determined by the initial orbit distribution of low-mass stars. Early segregation (either primordial or dynamical) of the most massive stars is simply invoked as a way to get a source of polluted material in the central regions of the cluster. We do not set up our clusters with primordial mass segregation, but dense clusters that are not primordially mass segregated will quickly segregate by two-body interactions anyway.³

Not including primordial mass segregation could potentially affect the early dynamical evolution of the cluster, but we argue below that it does not change our results about the long-term evolution. Early mass-loss from centrally concentrated massive stars can lead to a stronger expansion than the same mass-loss distributed throughout the body of the cluster. For tidally limited and strongly mass segregated clusters, Vesperini, McMillan & Portegies Zwart (2009) showed with N-body simulations that rapid dissolution may occur as a consequence of this early expansion and the associated larger flow of mass over the tidal boundary. However, they also showed that segregated clusters initially underfilling their Roche lobe survive the early expansion because much of this expansion occurs within the tidal radius and thus does not cause a significant loss of stars. In this case, the subsequent evolution, lifetime, and mass-loss rate of the cluster do not differ significantly from those of an initially non-segregated cluster. The main notable difference is that long-lived initially segregated clusters will tend to have a less concentrated structure, delaying core collapse. Note that the early mass-loss driven expansion of a mass segregated cluster is not homologous, as the massive (segregated) population loses more mass than the lower mass stars in the outer parts. The expansion is therefore more important in the core and has much less severe effects in the outer parts.

Because all our simulated clusters are initially tidally underfilling (as most GCs probably start) and because we are interested in long-term kinematic signatures that are more easily observed (and that more easily survive) in the outer parts of clusters, we conclude from the arguments above that our results would not be significantly affected if we had included primordial mass segregation.

In the scenario proposed by Bastian et al. (2013a), low-mass stars may accrete a significant fraction of their final mass, which could modify the mass function of the polluted population. However, we did not take this into account when building our N-body models for the accretion scenario, and instead assumed a standard mass function for both the polluted and the pristine population initially. This is partly for the sake of simplicity, but also because currently there are no clear predictions or robust observational constraints on the mass function of the polluted population following this gas accretion phase.

We only considered isochrone models with an isotropic velocity distribution. Bastian et al. (2013a) showed that when considering a radially anisotropic velocity distribution with the Osipkov–Merritt prescription (Osipkov 1979; Merritt 1985) and adopting an anisotropy radius equal to the half-mass radius, the fraction of stars that spend very little or no time in the core is in good agreement with the isotropic case. A radially anisotropic velocity distribution in the outer parts may represent more realistic initial conditions if clusters form via mergers of smaller clumps, but we will show in the next section that our simulated clusters rapidly develop radial anisotropy in the outer parts anyway.

To consider clusters with different amounts of angular momentum, we added rotation to the isochrone model by flipping the z-component of the angular momentum vector of a fraction of randomly selected stars having negative z-angular momentum within the cluster (which preserves virial equilibrium). The rotation curve initially has the same shape as the velocity dispersion profile, but as two-body relaxation proceeds it looks similar to what is obtained from self-consistent models including rotation (Varri & Bertin 2012), with a peak located near the half-mass radius of the cluster. We set up three systems (ACC0, ACC1, ACC2) with increasing values of λ matching those used for the initial conditions of the multiple generations scenario (see Table 1 and Section 3.1.2). The stability of these systems was verified by checking that their total angular momentum, their half-mass radius, and their global virial ratio stayed constant for several crossing times after the start of the simulation (for test runs without stellar mass-loss).

3.1.2 Multiple generations scenario

To emulate the initial conditions implied by the model of D'Ercole et al. (2008) and the hydrodynamical simulations of Bekki (2010, 2011), we set up self-gravitating models with a flattened and rotating stellar distribution embedded in a more extended spherical cluster. We preferred a purely dynamical approach over hydrodynamical simulations, partly to avoid the complications inherent to modelling the star formation process, but also because this approach is fast and gives us more control over the initial conditions.

We used the code magalie (Boily, Kroupa & Peñarrubia-Garrido 2001) included in the nemo⁴ Stellar Dynamics Toolbox (Teuben 1995, version 3.3.2). magalie can create multicomponent (e.g. disc+halo) stellar systems in near equilibrium using a method adapted from the buildgal procedure of Hernquist (1993). When embedding a disc within a spherical component, the method starts with the exact spatial density profiles and then approximates the phase-space distribution function of all components (and hence the velocity field) using moments of the collisionless Boltzmann equation. For all our initial conditions generated with magalie (see below), we verified the stability of the constructed systems by checking that the angular momentum in both populations, their half-mass radii and the global virial ratio stayed constant for several crossing times and revolution periods of the disc after the start of the simulation (for test runs without stellar mass-loss).

We considered two distinct populations in our toy models, namely a polluted and a pristine one. For models with net angular momentum, we represented the pristine stars (first generation) by a spherically symmetric halo following a Hernquist profile⁵ (Hernquist 1990) and the embedded polluted population (second generation) by an exponential disc with a scaleheight equal to 20 per cent of its scale-length. Like the isochrone model, the Hernquist model has a power-law density profile with a r⁻⁴ fall off in the outer parts, which describes young massive clusters very well (Elson et al. 1987) and therefore represents a good choice for realistic initial conditions.

We first explored how much angular momentum is expected to be lost from the first generation if ∼90 per cent of the initial cluster mass is rapidly removed from the outer parts. To do so, we set up systems in which the mass of the halo is 10 times larger than the mass of the disc (i.e. before mass-loss). We assumed that half of the mass of the disc is made of ejecta from the first generation, and the other half comes from accretion on to the cluster of diluting (pristine) gas with zero net angular momentum. From conservation of angular momentum, we thus require that the total angular momentum of the halo is 20 times larger than that of the disc. The amount of rotation in the halo is controlled by flipping the z-component of the angular momentum vector of a fraction of randomly selected stars having negative z-angular momentum. A system satisfying the above constraints and having λ = 0.091 is obtained by flipping the z-momentum of 65 per cent of the halo stars having a negative z-momentum and by setting the ratio between the scalelength of the Hernquist halo (a_halo) and the scalelength of the exponential disc (R_disc) to a_halo/R_disc = 9.37.

Fig. 1 shows the cumulative number and angular momentum distribution as a function of radius (in 3D) for the pristine and polluted populations of such a system. The vertical dashed line in each panel indicates the radius beyond which 90 per cent of the cluster stars reside. If we assume that the stars beyond this radius are removed instantaneously, we can see that this would not only preferentially remove the pristine stars (as required in this scenario, leaving the number ratio of the two populations close to unity), but it would also remove a large fraction of the angular momentum of the pristine population. In this particular setup, the resulting cluster would have a larger total angular momentum in the second generation than in the first generation by a factor of ∼3.

Figure 1.

Open in new tabDownload slide

Cumulative number and angular momentum distributions of polluted (blue curves) and pristine stars (red curves), as a function of radius, in an equilibrium ‘disc+halo’ configuration with λ = 0.091, M_halo = 10 M_disc, and J_{z, halo} = 20 J_{z, disc}. This captures the configuration of the multiple generations scenario before early violent mass-loss. The vertical dashed lines indicate the radius beyond which 90 per cent of the cluster stars reside. Top panel: cumulative distribution of angular momentum (z-component) for the polluted (disc) and pristine (halo) stars as a function of radius, normalized to the total angular momentum in each component. Middle panel: cumulative number distribution of the polluted and pristine stars as a function of radius, normalized to the total number of stars in each component. Bottom panel: ratio of the number of polluted and pristine stars within radius r.

We also looked at similar configurations but with different values of λ and a larger initial mass ratio between the first and second generations (a mass ratio of 10 is a lower limit if the internal mass budget problem is to be avoided in scenarios with multiple generations). For a fixed value of λ, configurations with a larger mass ratio have a disc that is less compact with respect to the halo (i.e. smaller a_halo/R_disc). In this case, the preferential loss of first-generation stars is less important when the outermost stars are removed. This can produce a system in which the total angular momentum of the first generation is still larger than that of the second generation, but it may also yield a fraction of second-generation stars that is too low (see Fig. 2). As λ is increased, a_halo/R_disc also gets smaller and it becomes increasingly difficult to remove a large fraction of pristine stars without losing a large fraction of polluted stars. If we take the simple dynamical considerations above at face value and if GCs can form with a large amount of angular momentum, a large initial mass ratio between the polluted and pristine populations (M_halo ≳ 20 M_disc) may be difficult to reconcile with the observed fraction of polluted stars in GCs.

Figure 2.

Open in new tabDownload slide

Top panel: percentage of polluted stars in the cluster after instantaneous removal of all the stars beyond a certain radius (in the multiple generations scenario), as a function of the percentage of the total mass removed. Bottom panel: ratio of total angular momentum in the remaining polluted and pristine stars, as a function of the percentage of the total mass removed from the cluster. In both panels, the different curves illustrate the effect of different initial conditions for the mass ratio between the two populations and the total amount of angular momentum (before mass removal).

In what follows, we focus on systems for which removing a large fraction of stars from the outer parts can leave the cluster with at least ∼50 per cent of polluted stars. In these cases, the total angular momentum in the remaining polluted stars is always larger than the total angular momentum in the remaining pristine stars (Fig. 2). For simplicity, in all our models mimicking the multiple generations scenario, we therefore assumed that the net angular momentum of the cluster is dominated by polluted stars and that the pristine population has no net angular momentum initially (J_{z, pristine}/J_{z, polluted} = 0). We also attributed an equal number of stars to each of the two populations. This captures the initial conditions after the early mass-loss phase implied by this scenario. In the remainder of this paper, when we refer to mass-loss, we are only concerned with the mass lost during the long-term evolution driven by two-body relaxation.

The initial conditions of our different simulations are summarized in Table 1. For a subset of these (models MGEN1, MGEN2, MGEN1-nosev), we set up an exponential disc embedded in a spherical Hernquist halo. The spatial extent of the disc with respect to the halo is then dictated by the choice of λ. One of these simulations (MGEN1-nosev), with λ = 0.091, is the one for which we did not include stellar evolution and which we followed until complete dissolution (see Section 3.2).

For comparison, we also considered another subset of initial conditions where the cluster has no net angular momentum and the second generation is represented by a spherically symmetric Hernquist model embedded in a more extended Hernquist halo. Unlike the rotating models above, the spatial extent of the second generation in these cases is not influenced by rotational support and does not depend on the properties of the first generation. Ultimately, it depends on the star formation process for this second generation. We thus decided to consider a few different concentrations for the second generation. In two cases (models MGEN0a and MGEN0b), the ratio of the radial scales of the two populations (a_pristine/a_polluted) was taken to be the same as the a_halo/R_disc values of our rotating models. We also considered a model in which the second generation is more concentrated (MGEN0c; see Table 1).

3.1.3 Comparison of the initial conditions of the two scenarios

As the polluted stars in the accretion model are assumed to be the ones crossing the core of the cluster, they are expected to be preferentially on radial orbits initially. In the presence of net angular momentum, this preference for radial orbits among the polluted stars also means that the mean rotational amplitude of these stars (at a given radius) will be lower than the mean rotational amplitude of the pristine stars. This is opposite to what is expected from the multiple generations scenario, where the angular momentum in the polluted population and the mean rotational velocity of the polluted stars is expected to be larger (Fig. 3).

Figure 3.

Open in new tabDownload slide

From top to bottom: number density, velocity dispersion in the z direction, velocity dispersion in the x direction, velocity anisotropy, and mean azimuthal velocity as a function of radius (in 3D) for polluted (blue) and pristine (red) stars. The initial conditions for model MGEN1 are shown in the left-hand panels, and those for model ACC1 are shown in the right-hand panels. The half-mass radius is indicated by a grey dashed line.

Note that we have neglected the effect of gas accretion on the orbits of stars in our simulations. In the early disc accretion scenario, the enriched gas is released by massive stars near the centre of the cluster before being quickly swept up (within a few Myr). As these stars located near the centre have low orbital angular momentum compared to the rest of the cluster stars, the enriched material that they expel would also have low specific angular momentum. The specific angular momentum of a star would thus typically be reduced following accretion, and its orbit would shrink (through conservation of energy and angular momentum). If anything, taking into account the effect of gas accretion would make the difference in the mean rotational velocity (at a given radius) between the polluted and pristine populations even larger, enhancing the difference already present between the initial conditions of the multiple generations and accretion scenarios.

In Fig. 3, we illustrate the similarities and differences in the initial conditions of the two scenarios by showing the number density, velocity dispersion (z and x components), velocity anisotropy, and mean azimuthal velocity (around the z-axis) profiles at time t = 0 for the polluted and pristine populations of models MGEN1 and ACC1. Similar plots are shown in Appendix A (available online as supplementary information) to compare the initial conditions of models MGEN0c and ACC0 (Fig. A1) and models MGEN2 and ACC2 (Fig. A2).

The 3D number density profile is already very similar in both of the scenarios considered at the beginning of the simulation. The global density profile of all our models will evolve towards a King-like profile (King 1966) at late times, and any memory of the specific model used to generate the initial conditions will gradually be erased as the density profile is reshaped by two-body relaxation and the tidal field.

The initial conditions of both scenarios are also generally characterized by a dynamically cooler and more centrally concentrated polluted population. Because the velocity dispersion of the cluster is higher towards the centre of the cluster, it may appear counter-intuitive that the more centrally concentrated population has a lower velocity dispersion. This is however a natural consequence of the more rapid fall-off of the density distribution of this population (see e.g. van den Bosch et al. 1999). Another way to see this is that because the polluted population is preferentially located in the central regions, it contains fewer stars on wide orbits that cross the inner regions of the cluster at high velocity.

For models MGEN1 and MGEN2, note however that the z-component of the velocity dispersion is by definition initially smaller for the disc population than for the spherical (pristine) population. The only case where the velocity dispersion of the polluted population can be larger than that of the pristine population is when looking at the component of the velocity dispersion in the plane of the disc (σ_x or equivalently σ_y) for these two models, from about the half-mass radius and beyond.

The anisotropy parameter is defined in the usual way as |$\beta = 1 - \frac{\sigma _{\rm t}^2}{2 \sigma _{\rm r}^2}$|⁠, where σ_t and σ_r are the tangential and radial velocity dispersions. The velocity distribution is thus radially anisotropic for β > 0. For all models with net rotation, the anisotropy profiles of the accretion and multiple generations scenarios are qualitatively similar. In both scenarios, the polluted population has a radially anisotropic velocity distribution, while the pristine population is either isotropic or even tangentially anisotropic in the inner regions for the early disc accretion scenario. Our accretion model without rotation (ACC0) displays the same behaviour, while both populations are (by construction) fully isotropic in the initial conditions of the multiple generations models without rotation (MGEN0a, MGEN0b, MGEN0c).

The mean azimuthal velocity profile (also referred to as the ‘rotation curve’ throughout this paper) is where the difference between the initial conditions of the accretion and multiple generations scenarios is the most striking. Apart from models without rotation, where the mean azimuthal velocity profile of both populations is obviously flat and zero, all of our initial conditions imply significant differential rotation between the polluted and pristine stars. At a given radial distance from the centre of the cluster, the mean azimuthal velocity is always larger for the pristine stars in the accretion scenario, while it is always larger for the polluted stars in the multiple generations scenario. The latter follows from our choice of J_{z, pristine}/J_{z, polluted} = 0 for the initial conditions of the multiple generations scenario. It would however remain true even with J_{z, pristine}/J_{z, polluted} = 1 because the polluted population would still be more centrally concentrated and would thus need to have a larger rotational velocity amplitude for its total angular momentum to be the same as the more extended pristine population.

3.2 Description of the runs

To follow the long-term dynamical evolution of the clusters, we used the code nbody6, a fourth-order Hermite integrator with Ahmad & Cohen (1973) neighbour scheme (Makino & Aarseth 1992; Aarseth 1999, 2003) and accelerated force calculation on NVIDIA Graphical Processing Units (Nitadori & Aarseth 2012).

As we are interested in modelling mixing processes for which the evolution is driven by two-body relaxation (see Vesperini et al. 2013), we did not include primordial binaries. This is a reasonable simplification because the post-core collapse expansion of the cluster and the accompanying evolution of the two-body relaxation time-scale is expected to behave in the same way irrespective of the details of the source of energy driving the expansion (e.g. primordial binaries or mass-loss from stellar evolution; Giersz & Heggie 2011; Gieles, Heggie & Zhao 2011).

With one exception (see below), all of our models include stellar and binary evolution (and associated mass-loss⁶) following the prescriptions of Hurley, Pols & Tout (2000) and Hurley, Tout & Pols (2002). All stars are given a metallicity of Z = 0.0035, but we ignore the effect of a spread in helium abundance among cluster stars on their evolution. For all our models with stellar evolution, we adopted a Kroupa (2001) IMF between 0.1 and 100 M_⊙, with a power-law slope of −1.3 below 0.5 M_⊙ and −2.3 above.

For the multiple generations scenario (Section 3.1.2), including stars up to 100 M_⊙ may seem rather extreme given that we follow the evolution of the cluster after the second generation has formed and the cluster has reached virial equilibrium, in which case high-mass stars (>8 M_⊙) from the first generation should have already exploded as supernovae. By ignoring that first-generation massive stars have already exploded, we actually overestimate the early mass-loss from supernovae by only ∼10 per cent. This may lead to an overestimate of the early expansion of the cluster by ∼10 per cent (e.g. Hills 1980), which is likely not very important when addressing whether signatures of the initial conditions can survive for a time-scale of the order of the relaxation time.

For the model without stellar evolution, stellar masses were sampled from a double power-law IMF between 0.1 and 1.2 M_⊙ (to mimic the present-day mass spectrum of GCs), with a break at 0.5 M_⊙ and power-law slopes of −1.3 and −2.0.

We evolved these clusters up to an age of 10.75 Gyr, which is close to the latest age determination of 47 Tuc from the properties of its white-dwarf cooling sequence (Hansen et al. 2013). At the end of the simulations, the modelled clusters are all reasonably consistent with the present-day characteristics of 47 Tuc (see Giersz & Heggie 2011): half-mass radii between 5 and 7 pc, central line-of-sight velocity dispersions between 10 and 13 kms⁻¹, Jacobi radii between 64 and 69 pc. Table 1 lists the fraction of the initial mass (M_f/M₀) and number of stars (N_f/N₀) left at the end of the simulation for all these models with stellar evolution.

For the model without stellar evolution, a cluster with 10⁵ stars, a virial radius of 2 pc and a mean mass 〈m〉 = 0.3458 M_⊙ was placed on a circular orbit with R_G = 10 kpc. Because the physical time-scale of this model is not set by the stellar evolution time-scale, we are free to scale both the mass and radius of the system to represent clusters with different relaxation times and Galactocentric radii.⁷ We ran this model until complete dissolution of the cluster, which allows us to explore the late stages of the dynamical evolution of multiple populations without being tied to a specific scaling.

The initial half-mass relaxation time is 0.9 Gyr for all our simulations with stellar evolution and 0.3 Gyr for our model without stellar evolution, using the equation from Spitzer (1987), which assumes an equal-mass cluster and spherical symmetry. One should however be cautious when directly comparing the age of the systems to this initial half-mass relaxation time, first because the half-mass relaxation time is not well defined for systems with a disc embedded in a spherical halo. We should expect two-body relaxation to be more efficient in systems in which stars are moving coherently (e.g. the disc component of our multiple generations systems), but this is not captured by the equation used. Our clusters also have a stellar mass spectrum and the mass function changes with time as a result of stellar and dynamical evolution. This will have an effect on the two-body relaxation time at any given point of the evolution.

4 LONG-TERM KINEMATIC IMPRINTS

4.1 Phase-space mixing

We now turn to the long-term dynamical evolution of clusters from the initial conditions described in the previous section. After 10.75 Gyr, our 47 Tuc-like models with stellar evolution have lost between 50 per cent and 58 per cent of their initial mass (∼10–25 per cent of their stars) through a combination of stellar mass-loss and tidal evaporation. At this stage, preferential loss of pristine stars has increased the ratio of polluted to pristine stars to N_polluted/N_pristine ≈ 0.85 for the accretion models and 1.03 < N_polluted/N_pristine < 1.42 for the multiple generations models.

Fig. 4 shows the time evolution of the Lagrangian radii of polluted and pristine stars (along with those of the whole cluster) for all the models with stellar evolution. Overlap of all Lagrangian radii for the two populations would indicate complete spatial mixing, but for none of the models in Fig. 4 are the polluted and pristine stars fully mixed at the end of the simulation. Only in model MGEN0a are both populations almost mixed, but this is the model in which the polluted population started the least concentrated with respect to the pristine population. Interestingly, the effect of the gravogyro instability in accelerating core collapse (e.g. Ernst et al. 2007) is seen when comparing models ACC0, ACC1, and ACC2, for which the initial conditions only differ in their total amount of angular momentum.

Figure 4.

Open in new tabDownload slide

Time evolution of the 1 per cent, 5 per cent, 10 per cent, 25 per cent, 50 per cent, 75 per cent, and 90 per cent Lagrangian radii for the polluted (blue lines), pristine (red lines) and all stars (black lines) for all of our models with stellar evolution (see Table 1).

Vesperini et al. (2013) showed that both internal two-body relaxation and mass-loss drive the evolution of the mixing of multiple populations. The relaxation time-scale is longer in the outer regions of the cluster, making two-body relaxation (and mixing) less efficient at these larger distances from the centre. This is amplified by the expansion of the cluster which slows down the rate of dynamical ageing with time. Mass-loss acts to slow the growth of the relaxation time, which eventually starts to decrease as the remaining mass and the radius of the cluster decrease. The Jacobi radius also decreases as the cluster mass goes down, and the less mixed outer layers are gradually stripped away, which accelerates the evolution towards complete spatial mixing.

Regardless of initial differences in the concentration of the polluted population, Vesperini et al. (2013) found that clusters always approach a state of complete mixing after losing 60–70 per cent of their mass. This is also what we find in model MGEN1-nosev, which was followed until complete dissolution. Fig. 5 shows the evolution of the Lagrangian radii of the polluted and pristine stars for this model as a function of the fraction of the initial mass left (M/M₀). All Lagrangian radii for the two populations overlap when M/M₀ ≈ 0.35.

Figure 5.

Open in new tabDownload slide

Evolution of the 1 per cent, 5 per cent, 10 per cent, 25 per cent, 50 per cent, 75 per cent, and 90 per cent Lagrangian radii for the polluted (blue lines), pristine (red lines) and all stars (black lines) as a function of the fraction of the initial mass left in the cluster for model MGEN1-nosev (see Table 1). The radius scale of this scalable model is expressed in N-body units.

In this work, we are mainly interested in the kinematic properties of multiple populations as a way to distinguish between the proposed scenarios. It is therefore important to verify if these kinematic imprints are erased on the same time-scale as full spatial mixing occurs. In their study of the long-term dynamical evolution of GCs with two stellar populations (one of them being more concentrated initially), Decressin et al. (2008) found that the loss of information on the stellar orbital angular momentum occurs on the same time-scale as spatial homogenization. To check this, we looked at the evolution of the two populations of our models in energy and angular momentum space.

Fig. 6 shows the distribution of energy and z-angular momentum for the two populations of model ACC1 at the beginning and at the end of the simulation. The energy of a star is defined as |$E = \frac{1}{2} m v^2 + m \phi _{\rm c}$|⁠, where ϕ_c – the potential energy due to all other cluster stars – excludes the contribution of the tidal field. One can see the obvious differences in the initial conditions of the two populations for the accretion scenario. Polluted stars typically having a lower energy (i.e. more negative – they are more tightly bound to the cluster due to their smaller distance from the centre and their low velocity). The net angular momentum of this model and the larger angular momentum of the pristine population are also apparent from the distribution of J_z. After 10.75 Gyr, differences are still present both in the energy and angular momentum distributions of the two populations.

Figure 6.

Open in new tabDownload slide

Distribution of energy and z-angular momentum for the polluted (blue) and pristine (red) stars of model ACC1 at t = 0 (top panel) and at the end of the simulation (bottom panel). The solid lines represent isodensity contours in the |E| versus J_z plane for each population.

Fig. 7 shows the same quantities but for model MGEN1. Again the polluted stars typically have a lower energy. The angular momentum distribution of the pristine stars is initially symmetric around J_z = 0, while the polluted stars have an angular momentum distribution skewed towards positive values of J_z for this multiple generations model with net angular momentum. After 10.75 Gyr, the two populations are again not fully mixed in energy and angular momentum space. Some exchange of angular momentum has occurred between the polluted and pristine stars, as the final J_z distribution of pristine stars is slightly skewed towards positive values.⁸

Figure 7.

Open in new tabDownload slide

Distribution of energy and z-angular momentum for the polluted (blue) and pristine (red) stars of model MGEN1 at t = 0 (top panel) and at the end of the simulation (bottom panel). The solid lines represent isodensity contours in the |E| versus J_z plane for each population.

As pointed out by Mastrobuono-Battisti & Perets (2013), partial relaxation and the exchange of angular momentum between the polluted and pristine populations can have consequences on the spatial morphology of the two populations (and the cluster as a whole) when the polluted population starts as a disc/flattened structure. We illustrate this in Fig. 8 by showing the cumulative fraction of stars as a function of cos i (where i is defined as the position angle of the star with respect to the positive z-axis) for the polluted, pristine, and all stars of model MGEN1 at the beginning and at the end of the simulation. The polluted population has not relaxed to a completely spherical shape (the isotropic distribution is represented by the dashed grey line in Fig. 8) by the end of the simulation. Slight flattening of the pristine population (which started isotropic) and of the cluster as a whole also follows from the exchange of angular momentum between the two populations, although the pristine population is obviously still less flattened than the polluted one. Studying the ellipticity profiles of multiple populations in GCs may reveal important clues about their formation. To establish whether this would allow us to distinguish between different scenarios, one should however also consider accretion models in which the cluster is initially flattened, but this is beyond the scope of this work.

Figure 8.

Open in new tabDownload slide

Cumulative fraction of stars as a function of cos i (where i is the position angle of the star with respect to the positive z-axis) for the polluted (blue), pristine (red), and all stars (black) of model MGEN1 at t = 0 (top panel) and at the end of the simulation (bottom panel). Stars in the plane of the disc (the x–y plane) have cos i = 0. The dashed grey line represents an isotropic spatial distribution.

Like the two examples shown in Figs 6 and 7, none of our 47 Tuc-like models are fully mixed in energy and angular momentum space after 10.75 Gyr (see figures in Appendix B, available online as supplementary information). To estimate when full mixing occurs, we turn to model MGEN1-nosev again. Fig. 9 shows the distribution of angular momentum and energy of polluted and pristine stars in this model at three different moments of the evolution of the cluster: when M/M₀ = 1.0, 0.5, and 0.35. After the cluster has lost 50 per cent of its initial mass, the distribution of angular momentum of the polluted stars is still skewed towards larger positive values of J_z compared to that of the pristine stars, similar to what was found for model MGEN1. The energy distribution of the polluted stars is also skewed towards lower (more negative) values compared to the pristine stars. When M/M₀ = 0.35, the shapes of the |E| and J_z distributions are essentially the same for polluted and pristine stars (up to a scaling factor due to the different number of stars in the two populations). This suggests that kinematic differences between polluted and pristine stars are erased on the same time-scale as full spatial mixing occurs.

Figure 9.

Open in new tabDownload slide

Distribution of energy and z-angular momentum for the polluted (blue) and pristine (red) stars of model MGEN1-nosev at t = 0 (top panel), when M/M₀ = 0.50 (middle panel) and when M/M₀ = 0.35 (bottom panel). The solid lines represent isodensity contours in the |E| versus J_z plane for each population.

In the next subsections, we look in more details into the specific kinematic signatures that are expected from each scenario when a cluster is not fully mixed.

4.2 Velocity dispersion

We now briefly discuss the evolution of the velocity dispersion profile of the two populations for our models with stellar evolution. Fig. 10 shows the time evolution of the velocity dispersion (the z-component) profile of the polluted and pristine stars for models MGEN1 and ACC1.

Figure 10.

Open in new tabDownload slide

Velocity dispersion (z-component) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGEN1 (left-hand panels) and model ACC1 at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines.

The variation of the two-body relaxation time as a function of radius and time is clearly seen in Fig. 10. After only about 1 Gyr, the central z-velocity dispersion is the same for the polluted and pristine stars in both models. Nearly 10 Gyr later, the region within which the velocity dispersion profile of the two populations is identical extends further out, but differences of the order of 1 or 2 km s⁻¹ persist in the outer parts of the cluster where the relaxation time is longer. The region where these differences persist appears to extend further in for model MGEN1, but when also considering simulations with other values of λ and different concentrations of the polluted population initially (see Figs C1 and C2 in the Appendix), this is not always apparent. It does not appear to be an unambiguous signature that would allow us to distinguish between the scenarios.

Note that for a dynamically young system, differences in the velocity dispersion of the polluted and pristine stars may persist even in the inner regions of the cluster. This could explain why the core velocity dispersions of the first and second generation stars are still different after 12 Gyr in the simulations of Mastrobuono-Battisti & Perets (2013), as their disc+halo model is tailored to ω Cen.

As discussed by Bekki (2010, 2011) and Mastrobuono-Battisti & Perets (2013), a lower velocity dispersion for the polluted population is one kinematic signature that can be left in a scenario with multiple generations, especially when looking at the component perpendicular to the initial plane of the disc. This is however clearly not unique to such a model, as we find the same signature in the accretion scenario, regardless of the initial value of λ. Therefore, the lower velocity dispersion reported for the enriched stars in 47 Tuc (Richer et al. 2013; Kucinskas et al. 2014) cannot be used to favour or discard one of the two scenarios considered here.

As mentioned in Section 3.1.3, the initial x (or y) velocity dispersion (i.e. the component in the rotation plane of the disc) of the polluted population can be larger than that of the pristine population in a multiple generations model with net angular momentum. Fig. 11 shows the time evolution of the velocity dispersion (the x-component) profile of the polluted and pristine stars for models MGEN1 and ACC1 (see Figs C3 and C4 in the appendix for cases with different values of λ initially).

Figure 11.

Open in new tabDownload slide

Velocity dispersion (x-component) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGEN1 (left-hand panels) and model ACC1 at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines.

At t = 10.75 Gyr, there is almost no trace of the initially larger x velocity dispersion of the polluted population in model MGEN1, as both the polluted and pristine populations have essentially the same x dispersion profile. In model MGEN2 (Fig. C4), a larger x-velocity dispersion for the polluted population is however still apparent at the end of the simulation in the outer parts of the cluster, but the difference is very small (<1 kms⁻¹).

For the accretion models and the multiple generation models without a disc (i.e. no net angular momentum), which are all spherically symmetric, the x and y components of the velocity dispersion are of the same magnitude as the z-component, so the arguments above about the z-dispersion still hold for these systems. The velocity dispersion (any component) is always smaller for the polluted (and more concentrated) population in these cases.

4.3 Velocity anisotropy

Motivated by the finding of a radially anisotropic polluted population in 47 Tuc (Richer et al. 2013), we now look at the evolution of the velocity anisotropy profile of the polluted and pristine stars in the accretion and multiple generations scenarios. Fig. 12 shows the time evolution of this profile for each of the two populations in models MGEN1 and ACC1.

Figure 12.

Open in new tabDownload slide

Velocity anisotropy (β) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGEN1 (left-hand panels) and model ACC1 at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines.

In both models, the polluted population stays more radially anisotropic than the pristine population in the outer parts of the cluster. The anisotropy profile of the polluted stars is characterized by an inner isotropic core followed by radial anisotropy that increases with radius. In the late stages, the radial anisotropy peaks at an intermediate distance from the centre and then decreases outwards, as expected when the effect of tides becomes important and stars on radial orbits are preferentially lost (e.g. Fukushige & Heggie 2000; Baumgardt & Makino 2003) and gain angular momentum due to the tides (Oh & Lin 1992). After 10.75 Gyr, the pristine population is also isotropic in the inner parts (over a larger radial extent than the polluted population) and shows some very mild radial anisotropy at larger radii, followed by a decrease of β and isotropy or mild tangential anisotropy in the outermost regions.

These conclusions do not change in models with different values of λ. Even when λ = 0 and both populations of a model with multiple generations start fully isotropic across the whole cluster, the more rapidly expanding (and initially more centrally concentrated) polluted population rapidly develops stronger radial anisotropy in the outer parts compared to the pristine population. The kinematic signature at the end of the simulation in this case is remarkably similar to all of our other models (Fig. 13). A more radially anisotropic polluted population is therefore not a unique signature of either scenarios, but another natural consequence of starting with a more centrally concentrated polluted population.

Figure 13.

Open in new tabDownload slide

Velocity anisotropy (β) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGEN0c (left-hand panels) and model ACC0 at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines.

4.4 Rotation

We now look at the long-term evolution of the bulk rotation of pristine and polluted stars, a promising diagnostic given the very distinct initial conditions of the accretion and multiple generations scenarios for this kinematic property (Fig. 3). To illustrate the differences between the two scenarios as a function of time, Fig. 14 shows the evolution of the mean azimuthal velocity profile for the two populations in models MGEN1 and ACC1.

Figure 14.

Open in new tabDownload slide

Mean azimuthal velocity as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGEN1 (left-hand panels) and model ACC1 at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines. The inset plot in the bottom right panel shows a zoomed view of a portion of the rotation curve for the final snapshot of model ACC1.

As the clusters evolve due to two-body relaxation and lose mass, angular momentum is transported outwards and the total angular momentum decreases. The peak of the rotation curve (located near the half-mass radius) also moves outward and its amplitude decreases. In all models with net rotation, as mixing progresses, the mean azimuthal velocity of the pristine stars eventually matches that of the polluted stars in the inner regions. The rotational amplitude in these inner regions is very low after 10.75 Gyr, but a small velocity gradient (of the order of 1 km s⁻¹ pc⁻¹) is still present across the core of the cluster, even in post-core collapse clusters (for example model MGEN1, see Fig. 4). This suggests that observational evidence for a velocity gradient/rotation in the inner parts of a cluster cannot be used to argue that the cluster is in the pre-core collapse phase of its evolution (cf. Fabricius et al. 2014).

In all our models with net rotation, the differential rotation of the pristine and polluted stars is maintained in the outer parts of the cluster after 10.75 Gyr of dynamical evolution. For the multiple generations scenario, the polluted stars rotate faster than the pristine stars, while the opposite is true for the accretion scenario. While the rotation curves of the two populations become indistinguishable in the inner regions of the clusters, differences of the order of 1–2 km s⁻¹ or more persist in the outer parts.⁹ A larger initial value of λ yields larger differences in the rotation curves of pristine and polluted stars (see Fig. D2 in the appendix).

In models MGEN1 and MGEN2 (Figs 14 and D2), the exchange of angular momentum from the polluted to the pristine population is seen through the subtle growth of the rotational amplitude of the pristine population (which has zero angular momentum initially) in the inner parts of the clusters. The rotation curve of the pristine stars also rises very slightly in the outer parts of these clusters as they evolve, showing a gain of weak prograde rotation at the level of ∼1 km s⁻¹. This effect is also seen in our models with no initial rotation (Fig. D1), but the effect remains small and is always confined to regions beyond the tidal radius of the cluster. Tidal effects are expected to cause a preferential loss of stars on prograde orbits (Hénon 1970; Keenan & Innanen 1975; Oh & Lin 1992; Fukushige & Heggie 2000), so we are most likely witnessing the stripping of stars on prograde orbits when seeing an increasing prograde rotation outside of the tidal radius. In any case, for our models with initial rotation, such a tidally induced rotational signature never compromises the clear differential rotation signal between the polluted and pristine stars.

5 PROSPECT OF DETECTING THE IDENTIFIED SIGNATURES

Based on the results of the previous section, comparing the velocity dispersion of the polluted and pristine populations may offer a way to distinguish between the scenarios kinematically, although the expected signal is weak. A larger velocity dispersion for the polluted population in the plane of rotation can be explained by a multiple generations model with net angular momentum, but not by the accretion model nor any model without net rotation initially. Observationally, this component of the velocity dispersion is best probed with radial velocities if the rotation axis is perpendicular to the line of sight, but best probed with proper-motion data if the rotation axis is along the line of sight.

If a larger velocity dispersion is measured for the polluted population in a given data set, the multiple generations model would be favoured and the same data should reveal a larger rotational amplitude for this polluted population. On the other hand, if a lower velocity dispersion is measured for the polluted population, it could mean that the accretion model prevails but could also be because the inclination of the rotation axis is not favourable enough to probe the velocity dispersion in the plane of rotation. In that sense, the velocity dispersion is perhaps not the most powerful probe of the formation scenario. On the other hand, it is clear that the differential rotation of polluted and pristine stars represents a unique kinematic signature that can help to constrain scenarios for the formation of multiple populations in GCs.

The expected strength of the signal in a cluster like 47 Tuc with a moderate amount of angular momentum initially is a few km s⁻¹ or less. This could be observable with precise ground-based radial velocity measurements from large high-resolution spectroscopic samples of RGB stars (bearing in mind inclination effects). As discussed in Section 2, tantalizing results in that respect have been obtained by Bellazzini et al. (2012), but they will need to be confirmed with larger data sets or at least with a more in-depth statistical analysis. Interestingly, in the three clusters for which these authors report different rotational amplitudes for the Na-rich and Na-poor samples, the pristine (Na-poor) stars always display a larger rotational amplitude, as expected in the accretion scenario.

The Gaia satellite will also allow us to obtain proper motions with a precision better than 2 km s⁻¹ for RGB stars in the outskirts of GCs, which could be sufficient to see a difference in the rotation of polluted and pristine stars in some clusters. Gaia will not be helpful for proper motions in the dense cores of GCs, but for the application discussed here, this is not a problem because the signal is expected to be stronger in the outer parts of clusters, where the two-body relaxation time-scale is longer. Proper motions from HST (e.g. Richer et al. 2013) may also be useful to uncover the expected signal. In this case, the analysis would be made easier by the differential nature of the signal, which eliminates the need for determining an absolute reference frame and avoids the associated complications, especially when the field of view of the observations is small (e.g. Anderson & King 2003).

To estimate the fraction of the Galactic population of GCs in which we can hope to detect differential rotation of the pristine and polluted stars, we must identify the clusters for which the different populations are not expected to be fully mixed. As we argued previously, these are the clusters which have lost less than ∼60–70 per cent of their initial mass (during the long-term evolution driven by two-body relaxation). Fig. 15 shows that the amplitude of the rotation curve is the same for the polluted and pristine stars of model MGEN1-nosev when the cluster has lost 70 per cent of its mass or more, supporting the idea that the differential rotation signature would be erased beyond that stage.

Figure 15.

Open in new tabDownload slide

Amplitude of the rotation curve (i.e. peak height of the mean azimuthal velocity profile) for the polluted (blue lines), pristine (red lines) and all stars (black lines) as function of the fraction of the initial mass left in the cluster for model MGEN1-nosev.

Figure 16.

Open in new tabDownload slide

Distribution of stellar mass and Galactocentric radius of Galactic GCs from the 2010 version of the Harris catalogue. Clusters above the dashed line are presumably in the expansion-dominated phase of their evolution and not fully mixed, while clusters below the line are expected to be in the evaporation-dominated phase. 47 Tuc is represented by a green square, and the fully spatially mixed cluster NGC 6362 by a red cross.

6 CONCLUSIONS

We focused on identifying a unique long-term kinematic imprint that would allow us to distinguish between different scenarios for how enriched material makes its way into stars in the context of multiple stellar populations in GCs. To do so, we presented a suite of N-body simulations of the long-term dynamical evolution of GCs, with initial conditions chosen to capture the distinct kinematic properties of two main pollution scenarios.

The first scenario (the multiple generations scenario) is one in which gas collects in a cooling flow into the core of the cluster, where a new generation of stars is eventually formed. In the other scenario (the accretion scenario), a fraction of the stars accrete material when passing through the core of the cluster, with no need for multiple star formation events. While in the current paradigm these scenarios are naturally associated with the models put forward by D'Ercole et al. (2008) and Bastian et al. (2013a) and the specific source of polluting material favoured by these authors (AGB stars and massive interacting binaries, respectively), our study does not concern the chemistry and our results would apply to scenarios that proceed in the same way dynamically, regardless of the source of pollution.

We showed that the two scenarios imply different initial conditions for the kinematics of the various populations. The most striking difference between the initial conditions of the accretion and multiple generations scenarios arises in the presence of cluster rotation. As a result of dissipative processes driving the polluted material to higher densities towards the centre of the cluster and conservation of angular momentum, an initially larger rotational amplitude for the polluted stars compared to the pristine stars is expected from the multiple generations scenario. In the accretion scenario, the polluted stars (the ones crossing the core) are preferentially on radial orbits and their initial rotational amplitude is expected to be smaller than that of the pristine stars.

We showed that initial differences in the kinematics of pristine and polluted stars can survive for a Hubble time of relaxation-driven dynamical evolution in a cluster with properties similar to those of 47 Tuc, especially in the outer parts of the cluster where the two-body relaxation time-scale is longer. In particular, some differential rotation of the pristine and polluted stars is preserved, which means that this unique signature can potentially be used to distinguish between the pollution scenarios considered. The expected strength of the differential rotation signal is typically a few km s⁻¹ in the outer parts (the exact value depends on the specific scenario, the initial amount of angular momentum of the cluster, the radius at which this is measured, and on the degree of mixing – i.e. the dynamical age of the cluster). This is challenging to measure but within the reach of current and future surveys of the internal kinematics of GCs (either through radial velocities or proper motions).

Initial differences in the velocity dispersion and velocity anisotropy profiles of polluted and pristine stars are also expected to survive for a Hubble time in a 47 Tuc-like cluster. These kinematic imprints are however remarkably similar in the accretion and multiple generations scenarios. A lower velocity dispersion and a more radially anisotropic velocity distribution for the polluted stars (common signatures of both scenarios) can therefore not be used to distinguish the two scenarios. That said, the multiple generations scenario could be favoured if a larger velocity dispersion is measured for the polluted population. This would stem from the initially larger velocity dispersion of this population in the plane of rotation.

We argued that full kinematic mixing of the two populations happens on the same time-scale as full spatial mixing. It is achieved when ∼60–70 per cent of the initial mass of the cluster has been lost. We identified a large subsample of Galactic GCs having large enough masses and Galactocentric radii to still be in the expansion-dominated phase of their evolution (Gieles et al. 2011). Clusters in this phase are still expanding within their tidal boundary and have not suffered substantial tidal evaporation. They are thus good candidates for observing the signatures of incomplete spatial and kinematic mixing, in particular the differential rotation of pristine and polluted stars.

We should mention that this work does not encompass all suggested frameworks for the formation of multiple populations. For example, Maxwell et al. (2014) suggested that GCs formed in the inner regions of dwarf galaxies can accrete gas on each passage through the centre of these galaxies. The accreted matter, eventually forming a second generation within the cluster, would in this case be a mix of pristine material from the intergalactic medium and surrounding dwarf, plus moderate velocity AGB ejecta retained by the potential well of the galaxy. These GCs would be pushed out on to larger orbits as star formation progresses, and may be stripped during mergers. While it avoids some of the problems of the other models, this scenario cannot explain the presence of multiple stellar populations in the Galactic GCs thought to form in situ (e.g. Forbes & Bridges 2010; Leaman, VandenBerg & Mendel 2013). It would however be interesting to see if this scenario can leave a characteristic kinematic signature. Similarly, we have not investigated the implications of another scenario in which second-generation stars are proposed to form in the decretion discs of fast rotating massive stars (Krause et al. 2013).

Finally, there is definitely scope to improve the toy models that we have set up to study the kinematics of multiple populations in different scenarios. For example, the effect of gas accretion on the orbits of stars could be included in a self-consistent way. It would also be interesting to consider the effect of the different initial helium content of the polluted and pristine stars on their evolution, and indirectly on the dynamical evolution of the cluster.

We thank the anonymous referee for a constructive report. We also thank Nate Bastian for useful comments, and Florent Renaud for pointing us to the code magalie. VHB is grateful to the ‘Fonds the recherche du Québec - Nature et technologies’ (FRQNT) for financial support through a postdoctoral fellowship grant. MG thanks the Royal Society for financial support in the form of a University Research Fellowship (URF) and the European Research Council (ERC-StG-335936, CLUSTERS).

REFERENCES

APPENDIX A: SUPPLEMENTARY FIGURES (INITIAL CONDITIONS)

This section presents two additional figures to complement Section 3.1 and Fig. 3 and illustrate initial conditions for models without net angular momentum initially (MGEN0c and ACC0; Fig. A1), and for models with a larger initial amount of angular momentum (MGEN2 and ACC2; Fig. A2).

Figure A1.

Open in new tabDownload slide

Same as Fig. 3, but comparing the initial conditions of models MGEN0c and ACC0.

Figure A2.

Open in new tabDownload slide

Same as Fig. 3, but comparing the initial conditions of models MGEN2 and ACC2.

APPENDIX B: SUPPLEMENTARY FIGURES (PHASE-SPACE MIXING)

The additional figures presented in this section show, for the models not illustrated in Section 4.1, the initial and final energy/z-angular momentum distributions of the pristine and polluted populations.

Figure B1.

Open in new tabDownload slide

Same as Figs 6 and 7 but for model ACC0.

Figure B2.

Open in new tabDownload slide

Same as Figs 6 and 7 but for model ACC2.

Figure B3.

Open in new tabDownload slide

Same as Figs 6 and 7 but for model MGEN0a.

Figure B4.

Open in new tabDownload slide

Same as Figs 6 and 7 but for model MGEN0b.

Figure B5.

Open in new tabDownload slide

Same as Figs 6 and 7 but for model MGEN0c.

Figure B6.

Open in new tabDownload slide

Same as Figs 6 and 7 but for model MGEN2.

APPENDIX C: SUPPLEMENTARY FIGURES (VELOCITY DISPERSION)

This section presents four additional figures to complement Section 4.2, Figs 10 and 11. These illustrate the evolution of the x and z velocity dispersion profiles for models without net angular momentum initially (MGEN0c and ACC0; Figs C1 and C3), and for models with a larger initial amount of angular momentum (MGEN2 and ACC2; Figs C2 and C4).

Figure C1.

Open in new tabDownload slide

Same as Fig. 10 but for models MGEN0c and ACC0.

Figure C2.

Open in new tabDownload slide

Same as Fig. 10 but for models MGEN2 and ACC2.

Figure C3.

Open in new tabDownload slide

Same as Fig. 11 but for models MGEN0c and ACC0.

Figure C4.

Open in new tabDownload slide

Same as Fig. 11 but for models MGEN2 and ACC2.

APPENDIX D: SUPPLEMENTARY FIGURES (ROTATION)

This section presents two additional figures to complement Section 4.4 and Fig. 14 and illustrate the evolution of the rotation curve for models without net angular momentum initially (MGEN0c and ACC0; Fig. D1), and for models with a larger initial amount of angular momentum (MGEN2 and ACC2; Fig. D2).

Figure D1.

Open in new tabDownload slide

Same as Fig. 14 but for models MGEN0c and ACC0.

Figure D2.

Open in new tabDownload slide

Same as Fig. 14 but for models MGEN2 and ACC2.