Kinematical evolution of multiple stellar populations in star clusters

TL;DR

This paper uses N-body simulations to explore how internal relaxation and external tidal forces influence the long-term kinematic evolution of multiple stellar populations in globular clusters, highlighting the importance of observational studies.

Contribution

The study provides new insights into the combined effects of internal dynamics and external perturbations on the kinematic properties of multiple populations in star clusters.

Findings

Differential rotation evolves under various spin-orbit conditions.

Kinematic and spatial mixing of populations is quantified.

Projection effects influence the interpretation of internal kinematics.

Abstract

We present the results of a suite of \Nbody simulations aimed at understanding the fundamental aspects of the long-term evolution of the internal kinematics of multiple stellar populations in globular clusters. Our models enable us to study the cooperative effects of internal, relaxation-driven processes and external, tidally-induced perturbations on the structural and kinematic properties of multiple-population globular clusters. To analyse the dynamical behaviour of the multiple stellar populations in a variety of spin-orbit coupling conditions, we have considered three reference cases in which the tidally perturbed star cluster rotates along an axis oriented in different directions with respect to the orbital angular momentum vector. We focus specifically on the characterisation of the evolution of the degree of differential rotation and anisotropy in the velocity space, and we quantify the process of spatial and kinematic mixing of the two populations. In light of recent and forthcoming explorations of the internal kinematics of this class of stellar systems by means of line-of sight and astrometric measurements, we also investigate the implications of projection effects and spatial distribution of the stars adopted as tracers. The kinematic and structural richness emerging from our models further emphasises the need and the importance of observational studies aimed at building a complete kinematical picture of the multiple population phenomenon.

Tables

Table 1. Table 1: Nomenclature of the N 𝑁 N -body simulations presented in this study. The symbol θ ω , i subscript 𝜃 𝜔 i \theta_{\rm\omega,i} denotes the initial value of the angle between the direction of the global angular velocity vector and the z 𝑧 z -axis, measured starting from the z 𝑧 z -axis (which points outwards from the orbital plane). For further details, see Section 2 .

Model ID	${\theta }_{\omega ,i}$ (degrees)
Theta0	0
Theta45	45
Theta90	90

Full text

††pubyear: 2015††pagerange: Kinematical evolution of multiple stellar populations in star clusters–Kinematical evolution of multiple stellar populations in star clusters

1 Introduction

A new picture of globular clusters is emerging from numerous observational studies which are revealing that these stellar systems are characterised by a rich spectrum of chemical and dynamical properties.

Spectroscopic and photometric studies have shown the presence of multiple stellar populations differing in the abundances of a number of light elements (see, e.g., Carretta et al. 2009a; Carretta et al. 2009b; Piotto et al. 2015; Milone et al. 2017). At the same time, large ground-based radial velocity surveys (see, e.g., Bellazzini et al. 2012; Fabricius et al. 2014; Lardo et al. 2015; Kamann et al. 2018; Ferraro et al. 2018) and astrometric studies (see, e.g., Watkins et al. 2015; Bellini et al. 2017; Bianchini et al. 2018; Vasiliev 2018; Libralato et al. 2018, 2019; Baumgardt et al. 2019; Sollima et al. 2019) have shown evidence of internal kinematical properties more complex than those included in the traditional view of globular clusters: many clusters are characterized by internal rotation and velocity anisotropy, dynamical ingredients neglected in most studies aimed at modeling clusters’ properties and their evolution.

A few investigations have further sharpened the dynamical picture of globular clusters by studying the kinematical properties of the multiple stellar populations within them and provided additional key constraints to understand their formation and dynamical evolution. Richer et al. (2013), Bellini et al. (2015), and Bellini et al. (2018), using HST proper-motion measurements for, respectively, 47 Tucanae and NGC2808, and Milone et al. (2018), using Gaia proper-motion measurements of 47 Tuc, have shown that second-generation (2G) stars (defined as those with enhancements in the abundances of elements such as Na, N, and He) are characterized by a more radially anisotropic velocity dispersion than first-generation (1G) stars.111Although we use the terms first- and second-generation to refer to the different populations, we point out that according to some formation models (see, e.g., Bastian et al. 2013; Gieles et al. 2018) all the populations form at the same time, and in those cases the two populations are more properly referred to as first- and second-population. Cordero et al. (2017) studied the rotational properties of the multiple populations of M13 and found that the extreme 2G population rotates more rapidly than the intermediate 2G. In M5, Lee (2017) found that 2G stars are characterized by a significant rotation while 1G stars show no rotation. An opposite trend was found instead in NGC 6752: for this cluster the analysis presented in Lee (2018) revealed that 1G stars rotate more rapidly than 1G stars. Finally, Dalessandro et al. (2018) found a puzzling difference in the velocity dispersion of 1G and 2G stars in the outer regions of NGC6362.

Although the observational investigation of the kinematic dimension of the multiple population phenomenon is still in its very early stages, the kinematical differences emerging from the first observational studies mentioned above clearly show the dynamical richness and complexity of multiple-population clusters and can be exploited to provide key constraints for models of their formation and dynamical evolution.

On the theoretical side, the numerous questions concerning the different aspects of the origin and dynamics of multiple stellar populations are still open and at the center of much attention. A number of numerical investigations have explored different aspects of the dynamics of multiple-population globular clusters, such as the evolution of the structural properties of different populations (see, e.g., Vesperini et al. 2013, 2018; Mastrobuono-Battisti & Perets 2013, 2016; Miholics et al. 2015; Fare et al. 2018), the dynamics of binary stars (Vesperini et al. 2011; Hong et al. 2015, 2016, 2019), and the effects of dynamical evolution on the stellar mass function (Vesperini et al. 2018). A few studies have started to specifically focus on the kinematical properties of multiple populations and attempted to distinguish kinematical features imprinted by the formation process from those which are instead shaped by the cluster’s long-term dynamical processes (see, e.g., Bekki 2010, 2011; Mastrobuono-Battisti & Perets 2013, 2016; Hénault-Brunet et al. 2015; Bellini et al. 2015; Dalessandro et al. 2018; Hong et al. 2019).

The goal of this paper is to explore the fundamental dynamical aspects of the evolution of the kinematical properties of multiple-population clusters and the connection between the evolution of the kinematical and structural properties. Specific attention will be devoted to the interplay between the effects of internal dynamics and those of the external tidal field of the host galaxy, and related implications on the internal kinematics of 1G and 2G stars. Although our $N$-body models are not aimed at a detailed comparison with observational data for specific globular clusters, we also explore how the intrinsic properties of our numerical simulations translates into corresponding projected observables, with a careful assessment of possible degeneracies and limitations in accessing the phase space of nearby star clusters.

The outline of the paper is the following: in Section 2, we describe the method and initial conditions of our $N$-body simulations, in Section 3, we present the results of our study, Section 4 includes a discussion of our results, and in Section 5 we state our conclusions.

2 Method and Initial Conditions

We have created a suite of direct $N$-body simulations using NBODY6 with GPU acceleration (Aarseth 2003; Nitadori & Aarseth 2012), and ran them on the Big Red II cluster at Indiana University.

The star clusters in our simulations move on circular orbits about the centre of their host galaxy. The potential of the host galaxy is modeled as that of a point-mass, and the equations of motions are solved in a co-rotating reference frame centred on the cluster (Heggie & Hut 2003), rotating with angular speed, $\Omega$, equal to that of the angular speed of the cluster’s orbital motion. The orbital plane is set up to be parallel to $x$-$y$ plane, and the cluster’s orbital angular velocity vector is parallel to the $z$-axis. All initial conditions begin with $N=$ 32 768 equal-mass stars, and stars are removed from the $N$-body system once they move beyond a distance equal to two times the cluster’s Jacobi radius, $r_{\rm J}$.

Our systems start with an equal amount of second-generation (2G) and first-generation (1G) stars. Since our investigation is focused only on the long-term evolution of multiple-population star clusters, we do not study here the early stages of a cluster’s evolution which might have driven the number ratio of 2G to 1G stars toward values close to those observed in Galactic clusters (see, e.g., Decressin et al. 2007; D’Ercole et al. 2008).

Following the prediction of a number of multiple-population formation scenarios we start with a 2G initially more compact than the 1G; our systems are characterized by internal rotation, and as predicted by the hydrodynamical simulations of Bekki (2010, 2011), we assume that the 2G is initially rotating more rapidly than the 1G population.

To set up our initial conditions, we produce the relevant rotating components from $N$-body realizations of the models by Varri & Bertin (2012). These distribution function-based equilibria are characterized by differential rotation, such that the rotation curve rises approximately linearly with radius in the cluster’s inner regions, reaches a peak in the intermediate parts, and then decreases and vanishes at the edge of the system. These configurations are oblate and flattened along a direction parallel to the rotation axis, and their rotational and structural properties are determined by four dimensionless parameters, $W_{0},\hat{\omega},\bar{b}$, and $c$, which define the central concentration, the rotation strength, and the shape of the rotation curve in the intermediate and outer regions of the models, respectively (see Varri & Bertin 2012, for further details on the properties of these models). The specific parameters we have adopted for our initial conditions are the following: for the 1G system, $W_{0}=3$, $\hat{\omega}=0.1$, $\bar{b}=0.1$, and $c=1$; for the 2G system $W_{0}=7$, $\hat{\omega}=0.3$, $\bar{b}=0.1$, and $c=1$. The two models have been scaled so that the initial ratio of the 1G to the 2G half-mass radius is equal to 5. The choice of the values attributed to the dimensionless parameters has been informed by a number of previous phenomenological studies in which such a family of models has been used to interpret the structural and kinematic properties of Galactic globular clusters (e.g., see Bianchini et al. 2013; Kacharov et al. 2014; Bellini et al. 2017). Similarly, the value of the ratio of the half-mass radii assumed here has been based on the current empirical knowledge of the spatial distribution of multiple populations in Galactic globular clusters, as based on existing photometric studies (e.g., Lardo et al. 2011). In all cases, the values adopted here should be considered only as representative, as a global modeling exercise (based on multi-component distribution function-based equilibria) of multiple-population star clusters has not been conducted yet in the current literature.

Similar to Tiongco et al. (2018), we present the analysis of our results using an inertial reference frame. In our analysis, we will refer to the global angular velocity vector, $\boldsymbol{\omega}$, calculated from an entire system of particles, and we define this vector to be the rotation axis of the system. We define the rotation axis’s direction in space using the angles $\phi_{\omega}$ and $\theta_{\omega}$: $\phi_{\omega}$ is the angle between the projection of $\boldsymbol{\omega}$ in the $x$-$y$ plane (the orbital plane) and the $x$-axis, measured in a counter-clockwise direction from the $x$-axis, and $\theta_{\omega}$ is the angle between $\boldsymbol{\omega}$ and the $z$-axis, measured starting from the $z$-axis (the axis pointing outwards from the orbital plane).

Table 1 lists the simulations in this study and the initial orientation of the rotation axis of the cluster, $\theta_{\rm\omega,i}$ (all models have their initial rotation axis in the $x$-$z$ plane). The parts of the cluster that intrinsically have zero rotation (the outermost regions) are not assumed to be tidally locked, i.e., the outermost regions initially have zero rotation in the inertial reference frame.

We initialize each $N$-body system such that the ratio of the cluster’s truncation radius is equal to the Jacobi radius, making its ratio of the half-mass radius to Jacobi radius ($r_{\rm h}/r_{\rm J}$, i.e., the cluster’s initial filling factor) equal to 0.106. We evolve each simulation until about 20–25% of the cluster’s initial mass is left.

3 Results

3.1 Rotation

3.1.1 Case of the rotation axis parallel to the orbital angular velocity vector

We begin the presentation of our results by considering first the simplest case among those under consideration (see Table 1). Therefore, we focus our attention on the evolution of the rotational properties of model Theta0, in which we assume that the internal and orbital angular velocity vectors are initially aligned (e.g., see Ernst et al. 2007; Kim et al. 2008). Later, we will discuss the evolution of the radial distribution of the angular momentum in the more general case, where the internal and orbital angular velocity vectors are not aligned. We will then move onto the discussion of the evolution of the degree of anisotropy of the velocity dispersion tensor in our $N$-body models.

Fig. 1 shows the evolution of the rotational velocity profiles of the two populations for the Theta0 model. We adopt a line of sight parallel to the rotation axis of the system (the $z$-axis) and the rotational velocity profile (as a function of the cylindrical radius) is calculated from the tangential components of the velocity vectors of the particles, as projected on the plane perpendicular to the chosen line of sight. The profiles are constructed by combining 5 snapshots around the desired time (except for time $t/t_{\rm rh,i}=0$, where only 1 snapshot is used), taking into consideration all particles within the Jacobi radius and binning them in cylindrical shells with heights encompassing the entire cluster.

In agreement with previous studies (see e.g., Einsel & Spurzem 1999; Ernst et al. 2007; Kim et al. 2008; Tiongco et al. 2017), we find that the overall cluster’s rotation decreases due to two-body relaxation redistributing angular momentum from the inner regions to the outer regions and the loss of angular momentum carried away by escaping stars.

The effects of two-body relaxation erase the initial differences between the spatial and kinematical properties of the two populations. The profiles shown in Fig. 1 illustrate how the rotational velocity profiles of the two populations evolve and become increasingly similar. Specifically, the 2G population’s peak rotational velocity decreases and its location moves to higher fractions of $R_{\rm h}$. Meanwhile, the changes in the rotation curve of 1G stars are more subtle and the location of its peak velocity has only a moderate shift toward the inner regions. We point out that, although over most of the cluster’s extension the 2G population is characterized by a larger rotational velocity than the 1G population, as 2G stars diffuse out of the central core and start populating the cluster’s outermost regions, their outer rotational velocity may be slower than that of the 1G population. This transition in the difference between the 1G and the 2G rotational velocities is therefore a kinematic manifestation and a consequence of the initial structural differences of the two populations.

After the rotational profiles of the two populations become indistinguishable, they share the same evolution towards smaller values of the rotational velocity, eventually converging towards a solid-body rotation, in agreement with what found in Tiongco et al. (2016b, 2017). As shown in Tiongco et al. (2016b), the combined effects due to the external tidal field and the preferential loss of prograde orbiting stars result in an internal rotation only partially synchronized with the external orbital rotation with an internal angular velocity $\omega\simeq 0.5\Omega$ (see also Claydon et al. 2017). In Fig. 1, we also show with a dashed line the solid-body rotation corresponding to $\omega=0.5\Omega$, to illustrate how the rotation curve is the combination of the intrinsic rotation curve (that is gradually being erased) and the growing solid-body rotation due to the effects of the external tidal field on the cluster’s kinematics (though the value of this slope can not be disentangled from the composite rotation curve until the final stages of dynamical evolution). The final panel of Fig. 1 shows that there is still some intrinsic rotation left in the cluster by this point in the simulation, not yet converging to $\omega\simeq 0.5\Omega$.

Further details of the evolution of the rotational properties of the 1G and the 2G populations are presented in Fig. 2; this figure shows the evolution of the difference between 1G’s and 2G’s rotational velocity measured at different radii from the cluster centre, and at the peak rotational velocity of each population. This figure illustrates how the difference between the 1G and the 2G rotational velocity evolves as a function of three different measures of dynamical evolution within the cluster: (i) time normalized to the initial half-mass relaxation time, (ii) the fraction of the initial mass left in the cluster, and (iii) the ratio of 3D half-mass radius of 1G to the 3D half-mass radius of 2G, $\rm r_{h,1G}/\rm r_{h,2G}$.

The picture emerging from this figure further illustrates how initial differences in the 1G and 2G rotational properties are gradually erased during a cluster’s evolution: the inner regions of the cluster are those where the local relaxation time is shorter and where memory of the initial rotational differences is lost more rapidly. The right panel of Fig. 2 shows the evolution of the rotational difference as a function of $\rm r_{h,1G}/\rm r_{h,2G}$; this ratio is a global measure of the spatial mixing of the two populations and allows us to illustrate the connection between the spatial and kinematical mixing of the two populations. It is interesting to point out that the mixing in the rotational profiles can be reached before complete spatial mixing. We illustrate this point further in Fig. 3, by showing the evolution of the projected number density profiles ($\Sigma(R)$) featuring the same times as Fig. 1 (however, it should be noted that the residual spatial radial gradient remaining when the two populations are kinematically mixed is very weak and likely difficult to observe).

3.1.2 Case of the rotation axis at a generic angle to the orbital angular velocity vector

Building on the understanding developed in the simpler case described in the previous subsection, here we follow the evolution of the kinematical properties of the $N$-body models in which we have relaxed the assumption concerning the initial alignment of the internal and orbital angular velocities. We now consider the case of a stellar system in which the initial internal rotation axis points in a generic direction with respect to the one defined by the orbital angular velocity vector. In Tiongco et al. (2018), we already showed that tidally perturbed rotating stellar systems characterized by an initial rotation axis oriented in a generic direction can develop a number of complex kinematic features, such as precession and nutation of the cluster’s rotation axis, a radial variation in the rotation axis’ orientation within the cluster, and counter-rotation between different parts of the cluster. As a more realistic extension of our physical set-up, we now explore how the initial differences in the structural properties of the 1G and the 2G populations affect the extent and development of our previously found features for the two stellar populations within the cluster model.

We first present in Fig. 4 the evolution of the rotation curves for each population of model Theta45. This figure demonstrates the importance of projection effects in the characterization of the rotation curve of the stellar system, as observed from different directions adopted for the line of sight. In the top rows, we consider a line of sight along the $z$-axis, which is the most favourable direction to detect the effects of (partial) tidal synchronization. Similar to Fig. 1, we denote with a dashed line the behaviour corresponding to $\omega=0.5\Omega$. Since model Theta45 begins with a smaller value of $v_{\rm rot}$ along the $z$-axis than model Theta0, we can see the Theta45 model converging to the expected partial synchronization at $t/t_{\rm rh,i}=25$, as shown in the final panel. In the bottom rows of Fig. 4, we consider a line of sight along the $x$-axis. Such a direction is indeed the least favourable one to characterize the rotational properties in the late evolution; unsurprisingly, we do not see the effects of (partial) tidal synchronization, therefore the rotation curve goes to zero in the outermost regions. The precession of the rotation axis causes the value of $v_{\rm rot}$ along the $x$-axis to become negative around $t/t_{\rm rh,i}=10$. Eventually, all signatures of rotation manifested along this line of sight are erased, as shown in the final panel.

Since precession can cause the direction of the rotation axis to vary greatly throughout the dynamical evolution of the tidally perturbed stellar system, we also show a more general characterization of the time evolution of a cluster’s rotational properties, no longer restricted to a specific line of sight. In Fig. 5, we show the evolution of the components of global angular velocity vector $\boldsymbol{\omega}$ for each population in all of our models. To obtain $\boldsymbol{\omega}$, we calculate the moment of inertia tensor, I, and the angular momentum vector, L, and solve the equation $\textbf{{L}}=\textbf{{I}}\boldsymbol{\omega}$. The study of the vector $\boldsymbol{\omega}$ allows us to focus on the cluster’s actual rotational properties, including the direction of rotation at any given time. We will further discuss the manifestation of the complex kinematical features of our clusters along different lines of sight and the associated projection effects in Section 4.

Similar to the results of Tiongco et al. (2018), the large-scale oscillations taking place over several relaxation times in the $x$ and $y$ components indicate that the rotation axis of each population is precessing. Such features are most prominent for model Theta45. Also present are smaller oscillations (taking place over fractions of relaxation times), indicating nutation of the rotation axis. The precession originates from the torque generated by the tidal field acting on the rotating, oblate cluster, which is initially oriented in a generic direction with respect to the symmetry introduced by the host galaxy’s tidal field. In Fig. 6, we further illustrate the precession and nutation of the rotation axis by plotting the evolution of the $x$ and $y$ components of $\boldsymbol{\omega}$ for model Theta45.

Fig. 5 also illustrates the gradual loss of angular momentum from each cluster model. As a result, the strength of rotation in all components decreases and the precession oscillations are damped. In the late stages of the long-term dynamical evolution of both models Theta45 and Theta90, the system shows a small residual rotation in the $z$-component. As discussed in the previous subsection, this feature corresponds to the solid-body rotation induced by the cluster’s interaction with the external tidal field (we wish to note that model Theta0 is still preserving some degree of its initial intrinsic rotation).

After all previous phases of investigation, we are now in the position to appreciate the evolution of the kinematical differences of the two populations: Fig. 5 also shows the time evolution of the components of the vector $\boldsymbol{\omega}$ for the 1G and the 2G populations, separately. A key result emerging from the $N$-body models that are oriented in a generic direction is that during the long term evolution of the system, the rotation axes of the 1G and 2G populations become misaligned. The difference in the direction of the 1G’s and 2G’s rotation axes is demonstrated in Fig. 7, where we show for models Theta45 and Theta90 the evolution of $\theta_{\omega}$ of each population.

The initial divergence in the directions of the rotation axes of each population (as seen in Fig. 7) is due to the 1G population being initially more extended and tidally filling; this implies that the interaction with the external tidal field and its effect in driving the cluster’s rotation towards a partially synchronized rotation around the direction perpendicular to the cluster’s orbital plane (the $z$ direction) acts more efficiently on the 1G population. As the two populations mix spatially and kinematically (at $t/t_{\rm rh,i}\approx$ 15–20), the rotation axes of each population become aligned again and eventually converge to a condition of partially synchronized rotation about the $z$-axis (and thus closer to $\theta_{\omega}=0$).

As we have shown in Tiongco et al. (2018), the non-trivial interplay between several internal dynamical processes and the effects of the external tidal field can lead to a number of complex kinematical features, including a radial variation in the orientation of the rotation axis. We show in Fig. 8 that both populations in the cluster evolve to develop a radial variation of the direction of the rotation axis (i.e., a variation of $\theta_{\omega}$ with radius); interestingly, multiple-population clusters are characterized by an additional complexity as the two populations do not follow the same radial profile of $\theta_{\omega}$. More specifically, the 1G’s rotation axis tends to point in a direction closer to that of the cluster’s orbital angular velocity vector (i.e. closer to 0 degrees) across a broader range of radii. As already discussed for Fig. 7, this is a consequence of the fact that the more extended spatial distribution of 1G stars causes this population to be more significantly affected by the external tidal field.

3.2 Velocity dispersion anisotropy

We now focus our attention on the evolution of the velocity anisotropy in our $N$-body models. Fig. 9 shows the evolution of the velocity anisotropy profiles of multiple populations in models Theta0 and Theta45. The initial conditions sampled from the Varri & Bertin (2012) equilibria already possess some initial intrinsic velocity anisotropy that is strongest when observing the cluster along a line of sight parallel to the rotation axis (see Varri & Bertin 2012, for further details). We first wish to focus on characterising the degree of evolutionary anisotropy (i.e., the one resulting from long-term, relaxation-driven effects). For this reason, we choose a line of sight perpendicular to the rotation axis for the analysis of model Theta0. On the other hand, for model Theta45 we adopt a generic direction as the line of sight, to highlight the features associated with both evolutionary and initial intrinsic velocity anisotropy. We define the anisotropy parameter as the ratio of the tangential velocity dispersion to the radial velocity dispersion, $\sigma_{\rm T}/\sigma_{\rm R}$., both measured from projected velocities on the plane perpendicular to the line of sight.

Drawing parallels to the results of Tiongco et al. (2016a), we find that the evolution of the anisotropy profile of the 1G population is similar to that of a tidally filling system, while the evolution of the anisotropy of the more compact 2G subsystem is similar to that of a tidally underfilling case. Specifically, as the 2G population expands due to two-body relaxation, it populates the outer halo with 2G stars on primarily radial orbits; this is clearly seen in Fig. 9 as the 2G’s evolution is characterised by a growing radial anisotropy in the cluster’s outer regions. As shown in Fig. 10, where we compare the ratios of the tangential velocity dispersion profiles for each population, $\sigma_{\rm T,2G}/\sigma_{\rm T,1G}$ and also the ratios of the radial velocity dispersion profiles, $\sigma_{\rm R,2G}/\sigma_{\rm R,1G}$, the larger radial anisotropy of the 2G population is due to the fact that this population has a smaller tangential velocity dispersion than the 1G. This behaviour is qualitatively consistent with that reported by Richer et al. (2013); Bellini et al. (2015); Milone et al. (2018), as resulting from state-of-the-art astrometric studies of several nearby Galactic globular clusters conducted with HST and Gaia DR2. The less concentrated and more tidally filling 1G population never develops a strong radial anisotropy in the velocity space. As in the previous case, this behaviour is in agreement with the findings of Tiongco et al. (2016a) for tidally filling systems.

As the system evolves and preferentially loses stars on more radial orbits (see e.g., Takahashi & Lee 2000; Baumgardt & Makino 2003; Tiongco et al. 2016a), the overall radial anisotropy of the system decreases due to relaxation-driven star escape. Moreover, the magnitude of the strongest radial anisotropy decreases as the Jacobi radius moves inwards (due to mass loss), stripping away the radially anisotropic layers and revealing the more isotropic layers closer to the centre of the cluster (see Giersz & Heggie 1997). The loss of preferentially radial orbits also creates some tangential anisotropy in the outermost regions.

Similar to the case of the rotational velocity profiles, the two population’s anisotropy profiles evolve to become increasingly similar as the two populations mix, after which the evolution proceeds similarly to that of a single population system. It is interesting to point out, however, that the two population may evolve to have a similar rotational profile while still be characterized by different anisotropy profiles (compare, for example, the anisotropy and rotation profiles for model Theta0: at $t=15t_{\rm rh,i}$ in Figs 1 and 9).

4 Discussion

As discussed in the introduction, a growing number of observational studies have revealed that some Galactic globular clusters have retained some memory of the kinematical differences expected in scenarios predicting that 2G stars form more concentrated and more rapidly rotating than 1G stars. The observational characterization of the kinematics of multiple populations, however, has been so far limited to a small set of cases and further studies will be essential to formulate further empirical constraints and a more complete picture of the dynamical dimension of the multiple population phenomenon.

On the theoretical side, the results of our direct $N$-body simulations have shown that the rich interplay between the dynamics of multiple-population formation, the dynamical processes affecting the cluster’s long-term evolution and the cluster’s interaction with the host galaxy’s tidal field can lead to a number of complex features and non-trivial differences between the kinematics of 1G and 2G populations. The awareness of these complexities has motivated the decision to focus our investigation exclusively on the interplay of only three fundamental physical ingredients: internal angular momentum, external tidal perturbations, and two-body relaxation. With such an illustrative goal in mind, the idealised $N$-body simulations presented in this study are appropriate to capture only the basic evolutionary processes at play. A detailed comparison with kinematic and photometric data of Galactic globular clusters would certainly require a careful exploration of a number of additional dimensions of the problem, both from the physical and numerical perspective, ranging from the effect of stellar evolution to the scaling with the number of particles. On the purely dynamical side, the presence of a mass spectrum would introduce much greater complexity to the definition of the problem (especially regarding the implications on the process of transport of the angular momentum), which would require a dedicated investigation. Finally, a number of degrees of freedom associated with the design of the initial conditions (such as the initial amount of angular momentum or the strength of the external tidal field, as specified here by the initial filling factor and the orbital angular velocity) would also deserve a separate exploration, by means of a suite of dedicated $N$-body simulations.

Although our $N$-body simulations are still idealised and only aimed at exploring the fundamental dynamical processes taking place in multiple-population clusters, in this section we provide a few examples illustrating how the rich kinematical properties of the stellar systems we have explored in this study would (or would not) be revealed by data typically collected in current observational studies. It is important to emphasise how the kinematic complexities emerging from the idealised star cluster models presented here and the possible challenges in their observational identification call for caution in the interpretation of observational studies based only on a partial characterisation of a cluster kinematics. Indeed, such caution further highlights the importance of building a complete picture based on the combination of line-of-sight velocities and proper motion measurements.

The lack of detection of kinematical differences can be due to the fact that in dynamically old clusters have effectively lost memory of all the fingerprints of the formation epoch and of those features usually produced during dynamical evolution. However, projection effects combined with the availability of only partial kinematical information may also limit our ability to identify existing kinematical differences or, possibly, reveal them but in a much weaker form.

Fig. 11 shows rotation curves of model Theta45 at $t=1.5t_{\rm rh,i}$, as calculated along three different lines of sight. The top panels show the rotation curves calculated as in Section 3, equivalent to using proper motion data, and the bottom panels show the rotation curves calculated equivalent to using radial/line-of-sight velocities. The figure demonstrates that observing the cluster along different lines of sight will show different strengths in the rotation curve depending on how close the line of sight is to being parallel (for proper motions) or perpendicular (for radial velocities) to the actual rotation axis. While the top panels show an obvious difference in the rotation curves of the two populations, the differences are less apparent using radial velocities. Additionally, the rotation curves again show that one population doesn’t necessarily rotate faster than the other in all regions of the cluster, i.e., the profiles show that 1G rotates faster than 2G in the outermost regions. This representative example highlights the need for observations of the internal kinematics of multiple populations in globular clusters across a broad radial range, ideally throughout the entire cluster extension.

When calculating the rotation curve using radial velocities, the position angle (PA) of the rotation axis is chosen as the one that maximises the amplitude of $\langle\Delta v_{\rm r}\rangle$ (see e.g., Bellazzini et al. 2012). However, in performing such an assessment, the size and radial distribution of the sample of stars under consideration can have a significant impact. Often, all of the stars are taken into account, but previous observational and theoretical studies have shown that position angle of the rotation axis can vary as a function of the radial distance from the cluster centre (see e.g., Gebhardt et al. 2000; Bianchini et al. 2013; Boberg et al. 2017; Tiongco et al. 2018). Thus, in the case of multiple populations co-existing in a star cluster, if the populations are not yet spatially mixed, each population may also have a different position angle of its rotation axis; this effect has already been partly noted in M13 by Cordero et al. (2017). In Fig. 12, we show the time evolution of the individual position angle of the rotation axis of each population, along three different lines of sight (this figure corresponds to Fig. 7). We find that the differences in PA are small enough so as not to make a significant difference in the rotation curves in each population. On the basis of this example, we conclude that, although possible misalignements of the PAs of the two populations appear to play only a secondary role in the characterisation of their corresponding kinematic observables, it is an aspect that, when possible, should always be empirically assessed, especially when interpreting observations of kinematical differences or lack thereof.

5 Conclusions

We have presented the results of a set of direct $N$-body simulations aimed at understanding the fundamental aspects of the evolution of the internal kinematics of multiple stellar populations in tidally perturbed star cluster models. We have followed the long-term dynamical evolution of multiple population clusters which, as suggested by hydrodynamical simulations of multiple-population formation, are composed of an extended, slower rotating first-generation (1G) population and a second-generation (2G) subsystem more spatially concentrated and more rapidly rotating than the 1G system in which it is embedded.

We focused our attention on the evolution of the degree of internal differential rotation and anisotropy in velocity space of the two populations, and we have explored their connection with the dynamical evolution of the cluster’s structural properties and the effects of the host galaxy’s tidal field.

We have studied the evolution of a tidally perturbed, rotating stellar system in which the orientation of the rotation axis is perpendicular to the plane of the cluster’s orbit around the host galaxy. We then considered the more general case of systems for which the initial rotation axis has a generic direction relative to orbital plane. In the latter case, as shown in Tiongco et al. (2018), a star cluster model can develop a number of complex kinematical properties resulting from the evolutionary interplay between relaxation-driven processes and external influences by the tidal field.

In all cases, we find that during a cluster’s long-term evolution, its overall internal rotation decreases due to the redistribution and loss of angular momentum, and that the rotation curves of the two populations become increasingly similar, eventually converging to a common rotational profile. The rotation curves of the two populations converge first in the inner regions of the cluster and subsequently in the outer regions. Before converging to a common behaviour, 2G stars are rotating more rapidly than the 1G component in the cluster’s inner regions; however, in the cluster’s outermost regions, where the 2G stars gradually move to as they diffuse from the inner regions, the 2G population is typically rotating more slowly than the 1G population.

Both populations eventually evolve towards a condition of internal rotation which is partially synchronized with the orbital motion, such that the internal angular velocity equals to half of the orbital angular velocity, as found in previous studies (Tiongco et al. 2016b, 2017; Claydon et al. 2017).

The results of the long-term dynamical evolution of multiple-population clusters starting with a rotation axis in a generic orientation relative the cluster’s orbital angular velocity vector show a number of complex kinematical features. As described in Tiongco et al. (2018), the orientation of cluster’s global rotation axis varies with time due to precession and nutation oscillations induced by the cluster’s interaction with the host galaxy’s tidal field. The rotation axis direction may also vary with radius within the cluster: the outer regions are more strongly affected by the external tidal field and their kinematical properties evolve to acquire a rotation around a direction closer to the direction of the cluster’s orbital angular velocity vector. The inner regions, on the other hand, are less affected by the external tidal field and tend to maintain the direction of the initial intrinsic rotation. As a consequence of the initial differences in the spatial distributions of the two populations, the impact of the tidal effects on the internal kinematical properties is stronger for the less spatially concentrated 1G stars, and 1G and 2G stars develop a misalignment in their global angular velocity vectors and different radial profiles of the direction of angular velocity vectors.

We also find significant differences in the velocity dispersion anisotropy of the two populations. The 2G system is initially compact, tidally underfilling, and develops a radial velocity anisotropy as 2G stars diffuse towards the cluster’s outer regions. The 1G system, on the other hand, is initially tidally filling and never develops a significant radial anisotropy. Differences in the anisotropy profiles of 1G and 2G can therefore be the kinematical manifestation of their different initial spatial properties. Similar to the evolution of rotation profiles, the two populations’ anisotropy profiles eventually become indistinguishable as the two populations mix.

Finally, we have discussed the possible observational limitations that can affect the ability to detect the kinematical differences and complex features that we have found in our $N$-body simulations. Projection effects may weaken or completely hide the signatures of the expected primordial or dynamically-induced kinematical differences. Moreover, the identification of many of the key fingerprints of the formation and dynamical evolution requires the ability to obtain the variation of the 2G and 1G kinematical properties as functions of the radial distance from the cluster centre. These limitations should be taken into account, and caution should be used in the interpretation of observational data.

The kinematical complexity emerging from our star cluster models further emphasizes the need and importance of observational studies aimed at building a complete kinematical picture of the multiple population phenomenon, as based on a combination of radial and proper motion measurements, ideally covering the broadest possible range of distances from the cluster centre.

Acknowledgements

We wish to thank the referee for a careful reading and review of the manuscript. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at IU is also supported in part by Lilly Endowment, Inc. MT acknowledges financial support from Indiana University’s President’s Diversity Dissertation Fellowship. ALV from a JSPS International Fellowship and Grant-in-Aid (KAKENHI-18F18787), a Marie Sklodowska-Curie Fellowship (MSCA-IF-EF-RI NESSY 658088), and the Institute for Astronomy at the University of Edinburgh.