© The Authors 2025

1 Introduction

Most globular clusters (GCs) show star-by-star variations in element abundances that point to distinct stellar populations. In general, observations show a first population (1P) with element abundances that are consistent with those of field stars having the same metallicity, and a second (2P) or more populations that present anomalous abundances in light elements (such as Na, O, Al, Mg, C, N, and He) while, in most cases, not showing any significant difference in heavy elements (see, Gratton et al. 2012, 2019; Bastian & Lardo 2018; Milone & Marino 2022, for reviews discussing evidence of multiple populations in spectroscopic and photometric studies)¹.

The multiple stellar populations in GCs do not show only chemical differences but also differences in their spatial and kinematic properties: 2P stars are found to be spatially more centrally concentrated and have internal kinematics characterised by a more rapid rotation and a radially anisotropic velocity distribution, whereas the 1P population is more isotropic or tangentially anisotropic. The strength of these dynamical differences decreases with the cluster’s dynamical age (see, e.g., Dalessandro et al. 2019, 2024; Libralato et al. 2023; Cordoni et al. 2020, 2025, and references therein for studies of the spatial and kinematic differences of multiple stellar populations in GCs, including examples of dynamically older clusters where differences between the 1P and the 2P dynamical properties have been erased by GC dynamical processes). The dynamical differences revealed by observations can provide some insight into the initial dynamical state of the populations emerging from the formation process, but investigations of the dynamics of multiple populations are necessary to establish a closer and more quantitative link between initial and present-day spatial and kinematic properties, since evolutionary processes can lead to different degrees of dynamical mixing during a cluster’s evolution.

Various formation scenarios have been proposed for multiple populations in GCs (see, e.g., Bastian & Lardo 2018; Gratton et al. 2019; Milone & Marino 2022, for a few reviews). While there is no consensus on many aspects of the formation of multiple populations, most models agree that the 2P would form with a more centrally concentrated spatial distribution than the 1P (see, e.g., D’Ercole et al. 2008; Bekki 2010, 2011; Gieles et al. 2018; Calura et al. 2019; Wang et al. 2020; Lacchin et al. 2022). In the last decade, several works have studied the long-term evolution of GCs with multiple populations through simulations focusing on the effects of dynamical processes on the evolution of the structural and kinematic properties of multiple populations (see, e.g., Bekki 2010; Mastrobuono-Battisti & Perets 2013; Vesperini et al. 2013, 2018; Vesperini et al. 2021; Hénault-Brunet et al. 2015; Tiongco et al. 2019; Sollima 2021; Lacchin et al. 2022; Hypki et al. 2025). Some investigations have also explored the possibility that some of the clusters with more complex chemical properties, characterised by a significant spread in iron, may result from the mergers of clusters with different iron abundances (Amaro-Seoane et al. 2013; Gavagnin et al. 2016; Khoperskov et al. 2018, see, e.g.,); see also Ishchenko et al. (2023) for a study of the possible collision history of Galactic GCs.

In this work, we follow the long-term evolution of GCs with multiple populations to characterise the dynamical mixing process in the energy-angular momentum phase space, E − L, and to connect the phase-space mixing process with observational tracers for the spatial and kinematic differences between the populations. We also explore the strength of the signatures of these dynamical differences over the dynamical history of the GCs and the possibility of establishing a link between the initial and the present-day dynamical properties of multiple populations.

The outline of the paper is as follows: Section 2 describes our simulations and the initial conditions of our models. In Sect. 3, we explore the mixing process through the behaviour of the stellar distribution in the energy and angular momentum phase space of each population. In Sect. 4, we analyse three observable properties that characterise multiple populations’ spatial and kinematic differences, and we compare them with the mixing process in phase space. Lastly, in Sect. 5, we discuss and summarise our results.

2 Initial conditions and simulations of globular clusters with multiple stellar populations

We follow the long-term evolution of seven simulated GCs with two distinct stellar populations. All models were evolved using the MOCCA code (Giersz et al. 2013; Hypki & Giersz 2013), which is based on the Henon’s Monte Carlo methods (Hénon 1971b,a) and includes prescriptions for stellar evolution (through the SSE and BSE codes Hurley et al. 2000, 2002), close stellar encounters (through the FEWBODY code Fregeau et al. 2004), and the dynamical formation of binaries. MOCCA accounts for the effects of a host galaxy tidal field by including a tidal radius truncation in the cluster (Hypki & Giersz 2013; Giersz et al. 2013, for further details on the MOCCA code). The time variation of the tidal fields associated with eccentric orbits or the tidal shocks due to, for example, Galactic disc crossings, are not included, and the study of these effects requires N-body simulations that properly account for a realistic modelling of the external tidal field (see, e.g., Cai et al. 2016; Berczik et al. 2025).

Five of the seven models presented here were also used in Livernois et al. (2024) to study energy equipartition in clusters with multiple populations. Here, we briefly describe the initial conditions of our models and refer to Livernois et al. (2024) for further details. The choice of initial conditions generally follows the results of simulations modelling the formation of multiple populations, finding that 2P stars form in a high-concentration system in the innermost regions of the 1P system (see, e.g., D’Ercole et al. 2008; Bekki 2010; Calura et al. 2019). All simulations start with the 1P modelled as a spherical system following the density profile of a King (1966) model with a central dimensionless potential of W₀ = 5, and the 2P modelled as a King model with W₀ = 7. All models start with a tidally filling 1P population extending to the cluster’s tidal radius, while the 2P population is concentrated in the innermost regions of the 1P system with an initial ratio of the 1P to the 2P half-mass radii equal to 20 or 10 (see Table 1). For all the simulations, we adopt an initial mass ratio between the populations of M_1P/M_2P = 4. We evolve clusters with 2 × 10⁶, 1 × 10⁶, and 0.5 × 10⁶ stars. All models have an initial tidal radius corresponding to that expected for clusters with the given mass at a Galactocentric distance of 4 kpc in a logarithmic potential with a circular velocity of 220 km/s. We also explore the case for a weaker tidal field (corresponding to a galactocentric distance of 8 kpc for the model with 0.5 × 10⁶ stars; N05M-wf in Table 1).

The initial conditions adopted follow the general trend of different possible formation scenarios predicting a centrally concentrated 2P population and are able to produce general properties at 12 Gyr consistent with those observed in GCs (see, Vesperini et al. 2021; Livernois et al. 2024; Hypki et al. 2022, and discussion therein). All simulations start with stellar masses of 1P and 2P stars following a Kroupa initial mass function (Kroupa 2001) between 0.1 and 100 M_⊙ and without primordial binaries.

The goal of this work is to study the dynamics behind multiple population mixing in phase space. The initial conditions explored allowed us to address some specific points concerning the mixing process, specifically:

• (a)

  Dynamical age and relaxation time: the different models we have simulated reach different dynamical ages and degrees of mixing at the end of the simulations, with the dynamical age increasing as initial masses decrease. This allows us to characterise the phase space properties of systems at different stages of their mixing dynamical history.

• (b)

  Primordial velocity anisotropy: all models except the N1M-ra model start with an isotropic velocity distribution. The N1M-ra model starts with a primordial radial anisotropy following an Osipkov–Merrit velocity anisotropy profile (Osipkov 1979; Merritt 1985). This anisotropy profile defines the kinematics of the whole cluster without distinction between populations. The centrally concentrated 2P stars populate the central regions of the clusters where the velocity distribution is close to isotropic, while the regions characterised by the stronger radial anisotropy are initially dominated by the 1P stars.

• (c)

  Different concentrations of the 2P: the five models in common with Livernois et al. (2024) have an initial half-mass radius ratio between the populations of r_h,1P/r_h,2P = 20. We included two new models with an initial concentration of r_h,1P/r_h,2P = 10 (named ’r10’), which made the 2P less concentrated.

Table 1 summarises the principal differences in initial conditions for the model analysed in this work. We saved each model state at every 100 Myr in order to have a continuous description of the time evolution of the properties we analyse here.

For the subsequent analysis, we measured a cluster’s dynamical age with the ratio of the cluster’s physical age to its half-mass relaxation time calculated at that age, n_h = t/t_rh. We use the Spitzer half-mass relaxation time (Spitzer 1987): $$t_{rh}=\frac{6.5\times {10}^2\mathrm{M}\mathrm{y}\mathrm{r}}{\ln (\lambda N)}{\left(\frac{M}{{10}^5{M}_{\odot }}\right)}^{1/2}\left(\frac{M_{\odot }}{⟨m⟩}\right){\left(\frac{r_h}{\mathrm{p}\mathrm{c}}\right)}^{3/2},$$(1)

where λ = 0.02 (Giersz & Heggie 1996), N is the total number of stars in the cluster, M the total mass of the cluster, ⟨m⟩ is the mean mass, and r_h is the half-mass radius.

In Tables A.1 and A.2, we report the properties of all our models at 7 Gyr and 12 Gyr, respectively. The final properties are generally consistent with those of many Galactic GCs. However, we strongly emphasise that these models are not meant to describe any specific observed GCs, nor do the initial conditions explored provide full coverage of the parameter space of possible initial and final conditions. The models are meant to provide a starting point for understanding the mixing process of multiple populations, as well as the current and new tools we are introducing to quantify the degree of mixing.

Fig. 1 Phase-space distribution of stars in the 1Ps and 2Ps for the models N2M and N05M at different times. At t = 0 Gyr, both clusters show a similar relative distribution of stars in the 1Ps (red) and 2Ps (blue), due to their similar half-mass radius and mass ratios. As the clusters evolve, the 2P expands towards the regions initially dominated by the 1P. By t = 12 Gyr, the N2M model is still partially mixed, while the N05M model has a significant degree of mixing. The insets show the marginalised distribution of energy and angular momentum, likewise highlighting the evolution of dynamical mixing. In all panels, the dashed lines follow the maximum allowed angular momentum at a given energy, defined by the angular momentum of the circular orbit of the corresponding energy.

3 Signatures of phase-space mixing

We studied the degree of dynamical mixing by tracing the time evolution of the distribution function (DF) for each stellar population. The simulations we used are spherically symmetric, and therefore, we expected the DF to be only a function of the energy (E) and total angular momentum (L = |L|). We identified the DF of each population as F_1P(E, L) for the initially extended 1P population, and F_2P(E, L) for the initially compact 2P population. As the cluster evolved, both populations, which were initially in different regions of the phase space of energy and angular momentum, diffused in phase space and eventually populated the same region. Once this occurred, we considered the populations to be fully mixed.

We followed the mixing of the stellar populations in each simulation by measuring the instantaneous energy and angular momentum per unit mass for each star every 100 Myr. We measured these values by taking the 3D positions and velocities of each star along with the potential energy given by the global mass distribution in the cluster: $$E_i=0.5\times (v_{r,i}^2+v_{\theta ,i}^2+v_{\varphi ,i}^2)+\mathrm{\Phi }\left(r_i\right),$$(2) $$L_i=r_i\times \sqrt{v_{\varphi ,i}^2+v_{\theta ,i}^2}.$$(3)

To compare simulations with different sizes and masses, we normalised the energy by Φ_h = 0.5GM_tot/r_h and the angular momentum by $L_h=\sqrt{0.5G{M}_{\text{tot}}{r}_h}$. Φ_h and L_h corresponds to the potential of a point mass of half of the total mass of the cluster at the distance of r_h, and the angular momentum of the circular orbit at r_h in the same point mass potential, respectively.

Figure 1 shows the DF of each stellar population for the simulated clusters with the lowest (N2M) and the highest (N05M) degree of mixing at three different epochs: 0 Gyr, 7 Gyr, and 12 Gyr. As described in Sect. 2, both clusters have different initial conditions but share the same relative ratios for each population’s half-mass radius and total masses, giving them similar normalised distributions in the phase space. Initially (t = 0), the 1P (red contours) is extended (larger values of L) and less concentrated (shallower total energy). Moreover, it occupies the upper regions of the phase space. The 2P (blue contours) is more concentrated and limited to the deepest part of the cluster’s potential. The inset panels show the initial difference between the configurations of both populations in the marginalised distributions of angular momentum and total energy. As the clusters evolve, both populations start to mix in phase space. The 2P ‘diffuses’ to larger energies and higher angular momentum values, a process also driven by the reduction of the potential caused by the cluster’s mass loss due to stellar evolution and stellar escape. The increase of the angular momentum values also goes hand in hand with the more extended orbits. At 7 Gyr (middle panel), we observe that model N2M still shows a significant difference between the two populations, while model N05M exhibits a higher degree of mixing, particularly in the energy distribution. The N05M model is also dynamically older with a dynamical age of n_h(7 Gyr) = 4.9 compared to model N2M with n_h(7 Gyr) = 1.0. At 12 Gyr, the populations in model N2M are still in the process of mixing. By contrast, both populations appear thoroughly mixed for model N05M, as the distributions of stars only show minor differences.

We constructed a grid on the energy-angular momentum space to quantify the differences in the distribution of both stellar populations. In each bin, we took the squared difference of the normalised distributions and then summed over all bins using $$S^2=\sum _i^{N_L}\sum _k^{N_E}{(\rho {(L_i,E_k)}_{2P}-\rho {(L_i,E_k)}_{1P})}^2\times \left(N_L{N}_E\right),$$(4)

where ρ(L_i, E_k)_1P and ρ(L_i, E_k)_2P are the number of stars in each (i, k) bin normalised by the total number of stars in each population to make the distributions comparable. While the initial value of S² might depend on the initial separation of the populations, as the populations mix, S² tends to zero.

Figure 2 shows the evolution of S² for all models in our sample as a function of time (panel a) and their dynamical age n_h (panel b). At 12 Gyr, the models show clear diversity in the degree of mixing for the clusters. The model N2M shows the least mixing with S² = 9.99, while the model N05M shows the most mixing at 12 Gyr with S² = 0.52. In between these two extremes, we find two groups with similar S²: first, the models N2M-r10 (S² = 4.07), N1M (S² = 3.62), and N05M-wf (S² = 3.59); second, the models N1M-r10 (S² = 1.77) and N1M-ra (S² = 1.61). All models follow a more similar evolutionary path of S² with respect to the dynamical age n_h (panel b). One substantial difference arises for the N2M model, where S² increases around n_h = 1, triggered by the cluster’s core collapse, which we discuss later in this section. The models N2M-r10 and N1M-r10 also follow a common evolutionary path; they begin with a less concentrated 2P and are initially more mixed. This translates to a lower initial value of S². The N2M-r10 and N1M-r10 models follow a similar but parallel evolutionary path to the N2M, N1M and N05M models. The N05M and N1M-ra models start with comparable values of S² as the N2M, N1M, and N05M-wf models, but the populations mix more rapidly.

3.1 Core collapse

Four models in our sample undergo core collapse within the 13 Gyr of evolution (see Fig. B.3). Among these, we see in Fig. 2 that the N2M model shows multiple peaks in S². While not as strong, the N1M, N05M, and N05M-wf models also show a peak in S² (panel a in Fig. 2). These peaks are consistent with the times of the first core collapse. Models N1M-ra, N2M-r10, and N1M-r10 do not undergo core collapse within the 13 Gyr of evolution.

We now discuss how the small peaks in S² found in some of the simulations are the fingerprints of the evolution of the innermost regions at the time of core collapse, and how they relate to the degree of spatial mixing of the 1P and 2P populations in the cluster’s central regions at that time. The N2M model undergoes core collapse at t = 7.06 Gyr (n_h(7.06 Gyr) = 0.97). At this time, the two populations do not show a significant mixing in phase space where the 2P dominates the cluster’s centre with a fraction of the total number of stars belonging to the 2P population within the 1Per cent Lagrangian radius equal to 0.87, while the global 2P fraction is equal to 0.52. Panel a in Fig. 3 shows the phase space at ~500 Myr before the first core collapse, showing that the populations are not yet mixed. On the other hand, the N05M model is considerably more mixed before it undergoes core collapse at t = 8.62 Gyr (n_h(8.62 Gyr) = 7.46), as shown in panel b of Fig. 3. The N05M model has a central 2P fraction of f_2P(r_1%) = 0.63 and a global fraction of f_2P = 0.6 at the time of the core collapse.

The difference in the central concentration of each population drives the observed small peaks in S² during the core collapse. Panel c in Fig. 3 shows the value of S² calculated using only the angular momentum distribution (i.e. marginalising over the energy space), S²(L), and the value of S² calculated using only the energy distribution, S²(E). In both simulations, the time evolution of S²(L) does not show any feature associated with core collapse, whereas S²(E) for the N2M simulation is characterised by the same peaks associated with the core collapse and the subsequent core oscillations already revealed by the evolution of S² in Fig. 2. As the core starts contracting in systems like the N2M models where 2P stars are the dominant population in the inner regions, more 2P stars follow the core contraction than 1P stars. This affects the potential energy distribution of the two populations and produces the observed peaks in S² and S²(E). The strength of the peaks is determined by their degree of mixing: the further the system is from mixing, the more prominent are the peaks in S² and S²(E).

Fig. 2 Evolution of the S2 parameter for all models. Panel (a) shows the evolution of S2 with respect to the physical time in the simulations. All models show different degrees of mixing at most times. At t = 12 Gyr (vertical line), model N2M is the least mixed, while model N05M is the most mixed in phase space. Panel (b) shows the evolution of S2 as a function of the number of half-mass relaxation times. The symbols show the value of S2 at 12 Gyr and the respective number of relaxation times.

Fig. 3 Effects of partial mixing during core collapse. Panels a and b: phase-space distribution of the models N2M and N05M at ~500 Myr before core collapse. The N2M model is far from fully mixed, while the N05M model has a significant degree of mixing. Panel c: evolution of the S2 parameter centred at the time of core collapse for the models N2M and N05M (tcc = 7.06 Gyr and tcc = 8.62 Gyr). The markers indicate the corresponding time of phase-space distributions in panels a and b. Paneld shows the value of the S2 parameter only for the energy component of the phase space. The peaks in S2(E) show the effects of core collapse and subsequent core oscillations for the model N2M, where the 2P dominates the inner regions. While for N05M model, the two populations are almost completely mixed at the time of core collapse, and the evolution of S2(E) does not show any feature associated with the core collapse.

3.2 Initially anisotropic model

The N1M and N1M-ra models start with the same initial conditions except for the initial radial anisotropy profile. The presence of primordial radial anisotropy impacts the overall evolution of the cluster. The N1M model has lost more mass and stars than the N1M-ra model. Furthermore, the N1M model undergoes core collapse at t = 7.54 Gyr, while the model N1M-ra continues to evolve towards core collapse at 13 Gyr. The delay of core collapse due to primordial velocity anisotropy has been thoroughly described by Breen et al. (2017) and Pavlík & Vesperini (2021) in the case of single-population clusters.

Figure 4 shows the S² parameter, as in Fig. 2, for the N1M and N1M-ra models, as well as the marginalised distribution differences, S²(L) and S²(E). The N1M-ra model develops more dynamical mixing than the N1M model within their common range of dynamical ages. Two effects drive the difference in mixing as indicated by the S² parameter. The initial velocity anisotropy in the N1M-ra model makes the angular momentum distribution for the 1P more compact and similar in extent to the 2P (see Fig. B.2). As the cluster evolves, the 2P develops radial anisotropy while the 1P slowly loses its anisotropy, becoming isotropic at about n_h = 3 (see Fig. 8 and discussion therein). On the other hand, the N1M model evolves quickly to have a radially anisotropic 2P and a tangentially anisotropic 1P, adding to the differences in the angular momentum distribution. The second effect is the delayed core collapse, which allows the two populations in the N1M-ra model to spatially mix at a steady pace. In contrast, in the N1M, and as discussed in Sect. 3.1, the core collapse delays the energy mixing if one of the two populations dominates the core.

Fig. 4 Evolution of the S2 parameter (a), as well as S2 measured over the angular momentum (S2(L), panel b) and energy (S2(E), panel c) spaces for the N1M and N1M-ra models. Primordial radial anisotropy increases the population mixing rate while keeping the cluster dynamically younger. The model N1M-ra does not undergo core collapse within its 13 Gyr of evolution.

4 Observable population mixing parameters derived from simulations

We discussed in Sect. 3 how differences in the phase-space distribution, traced by the S² parameter, can help quantify the degree of dynamical mixing of the two populations in a GC. However, the S² parameter can be measured only in numerical simulations. In this section, we explore three observable tracers for multi-population mixing, and we compare them with the S² parameter.

Stars considered for the analysis presented in this section are limited to those brighter than two magnitudes below the turn-off magnitude. This limit encompasses the range of current observational data for multiple populations in GCs (see, e.g., Libralato et al. 2022; Martens et al. 2023; Cordoni et al. 2020, 2025). For each model, and at each time, we calculated the projected 2D radius (R) and the radial and tangential proper motions (v_pmR, v_pmT) for fifty different lines of sight.

The analysis of these ‘observables’ is not meant to be directly compared with observations but rather to illustrate how these measurements can trace the degree of dynamical mixing of the stellar populations and how they are linked to the phase-space mixing discussed in the previous sections. Therefore, we did not include kinematic errors or completeness effects, and we assumed that the two populations were clearly distinguishable with no 2P (1P) star misclassified as belonging to the 1P (2P) population.

4.1 Spatial differences

To quantify the differences in concentration and spatial distribution between the two populations, we used the A⁺parameter (Alessandrini et al. 2016; Lanzoni et al. 2016). The A⁺parameter measures the area difference between the cumulative stellar count of two samples of stellar tracers within a specific radial range [usually to measure the differences in the spatial distributions of populations with different masses, see, e.g., Alessandrini et al. (2016); Lanzoni et al. (2016) for a study using A⁺to measure the segregation of blue-stragglers relative to main-sequence and red giant stars, and Weatherford et al. (2018); Della Croce et al. (2024) for a study using A⁺to quantify the degree of mass segregation of more massive main sequence stars as well as the link between A⁺and the black hole content in GCs]. In our case, we followed a similar formulation for multiple populations as the one introduced by Dalessandro et al. (2018), $$A^{+}\left(R_{max}\right)={\int }_0^1{\varphi }_{1P}\left(R^{\mathrm{\prime }}\right)-{\varphi }_{2P}\left(R^{\mathrm{\prime }}\right)d{R}^{\mathrm{\prime }},$$(5)

where ϕ(R′) is the cumulative radial distribution, R′ = R/R_max, and R_max is the maximum distance within which A⁺is calculated. In this work, we adopted R_max = 2R_h. We also normalised the cumulative radial profile so that ϕ(R_max) = 1.

Figure 5 shows the evolution of the $A_2^{+}$ parameter for all models as a function of time (panel a) and n_h (panel b). If we look at a specific time, such as at 12 Gyr, we find that our models develop different degrees of spatial mixing: models N2M and N2M-r10 show the highest difference in concentration (more negative $A_2^{+}$), while model N05M exhibits the most spatial mixing of the models ($A_2^{+}$ ~0). On the other hand, if we look at the dynamical age n_h, we find that the picture changes significantly. After n_h ~ 1, all models converge, and differences between the $A_2^{+}$ parameter for the various models are less significant.

It is interesting to note that the initial concentration differences between the r10 models and some of the r20 models mostly disappear by a dynamical age n_h ~ 1–2. Moreover, in all the models, the initial spatial differences between 1Ps and 2Ps significantly decrease during the clusters’ early evolutionary phases. Even in dynamically young clusters (for example in clusters with n_h ~ 1–3), these spatial differences are significantly smaller than those in the initial conditions imprinted by the formation process. These results show that the strength of present-day differences cannot be used to infer the extent of primordial differences and that even small present-day differences in dynamically young clusters may actually be the evolutionary outcome of multiple populations forming with strong dynamical differences [see also Tiongco et al. (2019); Vesperini et al. (2021); Sollima (2021) and the discussion in Cadelano et al. (2024)].

We do not find any clear signature of core collapse in the evolution of A⁺in Fig. 5. However, we can see some interesting behaviour if we look at the individual areas of the cumulative stellar distribution, particularly in the cluster centre (R < 0.5 R_h). In Fig. 6, we show the area of the cumulative stellar distribution within 0.5R_h for the models N2M (panel a) and N05M (panel b), following the comparison in Fig. 3. As we mentioned before in Sect. 3.1, the two populations in the N05M spatially mix before the cluster undergoes core collapse. At core collapse, the two populations increase their concentration at a similar rate. On the other hand, the model N2M is not spatially mixed by the time of core collapse, and 2P has a higher concentration than 1P, as shown by the larger value for the area of ϕ(R < 0.5R_h).

Fig. 5 $A_2^{+}$ Parameter evolution for all GCs as a function of time (a) and number of relaxation times (b). All models follow a similar evolutionary path as the two populations spatially mix. The panel a scatter reduces when we trace the cluster dynamical ages. After one relaxation time, the initial difference between the rh ratio of 20 (N2M, N1M, N1M-ra, N05M, and N05M-wf) and rh ratio of 10 models (N2M-r10 and N1M-r10) is indistinguishable.

4.2 Velocity anisotropy differences

In this section, we focus on the differences between the velocity anisotropy profiles of 1P and 2P stars. Several observational investigations of the velocity anisotropy profile of multiple populations have shown a 2P with varying degrees of radial anisotropy and a 1P that is isotropic or slightly tangential (see Bellini et al. 2015; Cordoni et al. 2020, 2025; Libralato et al. 2023; Dalessandro et al. 2024). Multiple studies investigating the evolution of the radial profiles of the velocity dispersion and anisotropy have shown that the observed velocity anisotropy profiles are consistent with those found in systems that assume initial conditions characterised by a tidally filling 1P system and a tidally underfilling 2P subsystem [see, e.g., Vesperini et al. (2021); Sollima (2021); Livernois et al. (2024); see also, e.g., Tiongco et al. (2016) and references therein for a similar evolution of the anisotropy in tidally filling and tidally underfilling single-population clusters].

With the goal of quantifying these differences, we adopted a new parameter that provides a global measure of the difference between the velocity anisotropy of the two populations by means of the area between the best-fitted velocity anisotropy models for 1P and 2P. We first calculated the velocity dispersion profiles for the radial and tangential components (projected on the sky) from the selected sample of stars. Then, we fitted a parametric model for the radial velocity dispersion and the velocity anisotropy. For the radial velocity dispersion, we used a fourth-order polynomial following the parameterisation in Watkins et al. (2022). For the tangential component, we used $${\sigma }_{\tan }^2={\sigma }_{\mathrm{r}\mathrm{a}\mathrm{d}}^2(1-{\beta }_{2\mathrm{D}}),$$(6)

where we define the projected velocity anisotropy profile, ${\beta }_{2\mathrm{D}}=1-{\sigma }_{\tan }^2/{\sigma }_{\text{rad\hspace{0.17em}}}^2$, with the parameterisation: $${\beta }_{2\mathrm{D}}(R)=\{\begin{array}{ll}\frac{{\beta }_i{R}^2}{(R_{\mathrm{a}}^2+R^2)}(1-\frac{R}{R_{\mathrm{t}}})& \text{\hspace{0.17em}for\hspace{0.17em}}R<R_{\mathrm{t}}\\ 0& \text{\hspace{0.17em}for\hspace{0.17em}}R\ge R_{\mathrm{t}}.\end{array}$$(7)

This parameterisation is based on the Osipkov-Merrit velocity anisotropy profile (Osipkov 1979; Merritt 1985) with the anisotropy radius R_a, but allows for tangential anisotropy (β_i < 0) and a maximum radial anisotropy smaller than max(β_2D(R)) < 1. We also included a truncation term, (1 − R/R_t), so that the anisotropy profile could once again be isotropic (β_2D = 0) at the truncation radius R_t. The truncation term allows us to properly characterise our simulations’ velocity anisotropy profile at large radii. For the purpose of this analysis, we fixed the anisotropy truncation radius at the respective tidal radius of each simulated cluster’s snapshot. We then calculated the anisotropy area A_β for each population by first normalising the radius with the tidal radius and then integrating over β_2D(R).

Panels a and b in Fig. 7 show the best-fit velocity anisotropy profiles for the N2M model at two different times (t = 4 Gyr and t = 12 Gyr, respectively). At 4 Gyr, the 2P shows radial anisotropy with an anisotropy area of A_β,2P = 0.282, while the 1P is close to being isotropic with A_β,1P = 0.011. At 12 Gyr the 2P becomes slightly less radially anisotropic with A_β,2P = 0.235, and the 1P is now tangentially anisotropic with A_β,1P = −0.132. In panel (c), we show the time evolution of the anisotropy areas for model N2M. Both populations began isotropic with A_β = 0; during the first 100 Myr of evolution, both populations developed radial anisotropy (with a stronger anisotropy for the 2P). The development of radial anisotropy follows early mass loss due to stellar evolution, followed by a shallowing of the potential and an expansion of the cluster. The difference in strength of radial anisotropy is a consequence of the initial extension of the populations and the tidally filling and underfilling conditions for the 1Ps and 2Ps, respectively. As the cluster evolves, the 2P steadily loses radial anisotropy, while the 1P becomes isotropic and then develops tangential anisotropy.

Panel d in Fig. 7 shows the anisotropy area difference between the populations for the N2M model. This difference is defined as $$\mathrm{\Delta }{A}_{\beta }=A_{\beta ,2\mathrm{P}}-A_{\beta ,1\mathrm{P}}.$$(8)

If both populations are close to isotropic (A_β,2P ~ A_β,1P ~ 0) or have similar anisotropy profiles (A_β,2P ~ A_β,1P ≠ 0), then ΔA_β is close to zero, we consider the populations to be kinematically mixed. On the other hand, as ΔA_β becomes different from zero, the populations are further from being mixed. This could be the case when 2P has a strong radial anisotropy, and 1P is isotropic (panel a in Fig. 7), or the two populations have opposite velocity anisotropy (as in panel b in Fig. 7, where 2P is radially anisotropic and 1P is tangentially anisotropic). We see that, for the N2M model, the anisotropy area difference remains relatively constant, which is consistent with it being the furthest from mixing at the end of the simulations.

We show the evolution of the anisotropy area for each population and the anisotropy difference for all models in Fig. 8. Our choice of initial conditions affects the populations differently. In panel a, the 1P for all models that are initially isotropic follows the same path: starting at A_β = 0, developing a very mild radial anisotropy, becoming isotropic, and then developing tangential anisotropy. Only the N05M model returns to being close to isotropic at the end of the simulation. On the other hand, the N1M-ra model’s 1P begins with primordial radial anisotropy, which translates to an initial A_β = 0.4. The cluster’s 1P proceeds to lose radial anisotropy and becomes isotropic at n_h = 3. The evolution of the velocity anisotropy in the 1P is driven by its tidally filling condition and is consistent with tidally filling single population simulations (see, e.g., Tiongco et al. 2016). In panel b, the 2P shows two main paths of evolution driven by its initial concentration: as the 2P expands, clusters with an initially more concentrated 2P develop a more substantial radial anisotropy (A_β ~ 0.4) than those with an initial concentration of r10 (A_β ~ 0.3). As for the N1M-ra model, we note that, given the primordial velocity anisotropy profile of this model, the 2P stars populate mainly the inner isotropic regions, while the radial anisotropy mainly affects the 1P populations extending in the outer anisotropic regions.

In panel c of Fig. 8, we show the area difference for each model. Following the discussion for panels a and b, the imprint of the initial conditions is still present on the area difference. The initially isotropic models with r20 concentration develop a ΔA_β ~ 0.35 after the first 100 Myr, while the models with r10 develop a ΔA_β ~ 0.25. We see that for the two concentration groups, ΔA_β slowly decreases as the clusters become dynamically older, showing that differences in velocity anisotropy persist over several relaxation times. Only after n_h = 6 does the N05M model show a continuous decline in ΔA_β, approaching complete kinematic mixing. On the other hand, the N1M-ra model has a ΔA_β = 0.05 after the first 100 Myr, starting in a ‘close to kinematic mixing’ state. This is because the 2P develops radial anisotropy, while the 1P was already radially anisotropic. As the cluster evolves dynamically and the 1P becomes isotropic, the ΔA_β value of the N1M-ra model increases and becomes similar to those with an initial concentration of r10. By the end of the simulation, as the 1P of the N1M-ra model becomes tangentially anisotropic, its ΔA_β value becomes consistent with the r20 models.

In Dalessandro et al. (2024), we estimated the anisotropy area (referred to as α_β in that paper) and the anisotropy area difference (referred to as $({\alpha }_{\beta }^{SP}-{\alpha }_{\beta }^{FP})$ in that paper) for a sample of Galactic GCs based on their observed kinematics, using the same definitions described here. We presented the N1M, N1M-ra, and N05M models to compare simulations with observations. Interestingly, we find a consistent trend between the observations and simulations in the evolution of these quantities. While further comparisons of the observed GCs and the simulations are beyond the scope of our current sample of simulations, we stress here the usefulness of the anisotropy area difference to trace the mixing process.

Fig. 6 Area of the cumulative stellar distribution within 0.5Rh (ϕ(<0.5Rh)) for the two populations in models: (a) N2M and (b) N05M. As in Fig. 3, we centre the time evolution at the respective time of the first core collapse. The two populations in model N2M have not achieved complete spatial mixing and ϕ2P(R < 0.5Rh) > ϕ1P(R < 0.5Rh). The model N05M has already attained a significant spatial mixing and ϕ2P(R < 0.5Rh) ~ ϕ1P(R < 0.5Rh). Furthermore, for model N05M, both populations show the same rate of increased concentration at core collapse, which explains the lack of signatures in S2(E) (see Fig. 3).

Fig. 7 Velocity anisotropy area (Aβ) and area difference (ΔAβ) for the N2M model. Panels a and b show the velocity anisotropy for the two populations (red upper triangles, PO1P, and downwards blue triangles, PO2P) in the N2M model at 4 Gyr and 12 Gyr, respectively. The hatched regions represent the respective anisotropy area for each population from the best-fit velocity dispersion and anisotropy parameterisation (see main text). Panel c shows the evolution of the anisotropy areas for the two populations over time; as the 2P becomes less radially anisotropic, the 1P becomes tangentially anisotropic and keeps a constant difference with the 2P, as shown in panel d, where we show the evolution of the anisotropy area difference. The vertical lines in panels c and d represent the times for which we show the velocity anisotropy profiles in panels a and b.

Fig. 8 Anisotropy area (Aβ) and area difference (ΔAβ) as described in Fig. 7, for all models in our sample with respect to their dynamical ages. For better visualisation, we separated nh into a linear section (nh ≤ 6) and a logarithmic section (nh > 6).

4.3 Projected angular momentum differences

In Sect. 3, we showed that differences in the marginalised angular momentum can also trace kinematic mixing (see Figs. 1, 3, 4, and B.3). Furthermore, differences in the velocity anisotropy imply differences in the radial and tangential velocity dispersion components, which can manifest in the angular momentum distribution. As initially suggested in Sect. 3.2, stars on more radial orbits have smaller angular momentum at a given radius.

In this section, we explore the connection between the projected angular momentum L_z = v_T × R (where v_T is the tangential velocity on the plane of the sky, and R is the projected distance from the cluster’s centre), and the velocity anisotropy. In a non-rotating cluster, the L_z distribution at any given radius is centred at zero, and, therefore, we focused on the dispersion of the L_z distribution, σ_Lz, as an observational tracer of the 3D angular momentum. As in Sect. 3, we normalised by the half-mass angular momentum L_h of a circular orbit in an equivalent point-mass potential.

In order to provide some initial guidance on the link between σ_Lz and the cluster’s velocity anisotropy, we show in Fig. 9 the radial profile of σ_Lz for four single-population King models with different radial anisotropy profiles built with the Limepy code (Gieles & Zocchi 2015, 2017). The value of σ_Lz at a given radius in the cluster’s outer regions decreases for more radially anisotropic models. This decrease in σ_Lz for more radially anisotropic velocity distributions is consistent with the expected smaller fraction of orbits with larger angular momentum (corresponding to a narrower distribution of L_z for more anisotropic models).

Following the correlations observed in the single-population King models, the differences in velocity anisotropy between the stellar populations discussed in Sect. 4.2 should also appear as differences in σ_Lz. Fig. 10 shows the σ_Lz profile for the two populations of the N2M and N05M models at 12 Gyr (panels a and b). The N2M model, which exhibits the largest velocity anisotropy area difference, also displays distinct σ_Lz profiles, and as expected, σ_Lz (1P) ≥ σ_Lz (2P). The N05M model shows almost no differences between the σ_Lz profiles, consistent with being close to complete mixing. Considering how the differences in σ_Lz appear more significant outside of the half-light radius, we follow the time evolution of the difference between the 1P and 2P values of σ_Lz at r = 2R_h for all epochs in our simulations. Panel c shows the time evolution of Δσ_Lz for models N2M and N05M. This plot shows that, as the cluster evolves and mixes, Δσ_Lz approaches zero and provides another tracer to measure the kinematic mixing of the two populations.

Figure 11 shows the evolution of Δσ_Lz for all the models in our sample in terms of physical time and dynamical age. The models N2M and N 05 M show the most and the least differences in σ_Lz between the populations. The models with initial r10 concentration begin with smaller values of Δσ_Lz than the r20 models, and after 1 Gyr, the models with initial r10 concentration overlap with the N1M-ra model. The latter also moves to a larger difference in σ_Lz, consistent with 1P evolving from a radially anisotropic state to isotropic and then tangentially anisotropic. As we examine the dynamical age, we see that the evolutionary paths change. Instead of exhibiting a set of parallel tracks, all the models converge into a single path. The models N1M-ra, N1M-r10, and N2M-r10 join at a single path at approximately n_h = 0.1 early in the dynamical evolution of the clusters. As we measure the differences at 2R_h, we see degeneration between the effects of a weaker radially anisotropic 2P (N1M-r10 and N2M-r10) and two populations with different degrees of radial anisotropy (N1M-ra, the 1P becomes isotropic at n_h = 3, see panel a in Fig. 8). The rest of the models follow a single path and begin to overlap at n_h = 4. While most of our models do not get past this point, we expect them to follow the path of the N05M model.

Fig. 9 Projected angular momentum dispersion σLz profiles for King models with different velocity anisotropy profiles (panel a) and their respective velocity anisotropy profiles (panel b). As the King models become more radially anisotropic (smaller values of ra/rh), the projected angular momentum dispersion decreases. The insets on panel a show the distribution of the projected angular momentum at r ~ 2Rh for the models with the most and the least radial anisotropy. The grey area in panel a shows the range of radial distances of stars used to calculate the histograms in the panel. We normalised σLz by the half-mass angular momentum Lh, corresponding to the angular momentum of a circular orbit of an equivalent point-mass potential (see Sect. 3).

Fig. 10 Projected angular momentum dispersion σLz profiles for both populations in the N2M model (panel a) and the N05M model (panel b) at 12 Gyr. The vertical grey solid lines in panels a and b represent the difference in σLz between the populations at R = 2Rh. Panel c shows the time evolution of σLz difference at 2Rh; the closer this value is to zero, the closer the two populations are to kinematically mixing.

Fig. 11 Evolution of the projected angular momentum difference at 2Rh for all models in our sample as described in panel c of Fig. 10. In panel a, this is with respect to the physical time of the simulation, and in panel b this is with respect to their dynamical age. Once we consider the dynamical age, all models converge into the same evolutionary path, independently of their initial conditions.

4.4 Connection with the S² parameter

We explored three different observational parameters to quantify the degree of mixing of the multiple populations. All three trace different aspects of the mixing process: $A_2^{+}$ follows spatial mixing, whereas ΔA_β and Δσ_Lz trace differences in the velocity distribution anisotropy of the two populations and follow kinematic mixing. In Sect. 3, we defined the S² parameter as a way to describe overall phase-space mixing. Now, we compare the S² parameter – which cannot be observed directly – with the three observable tracers of mixing.

Figure 12 shows the three projected quantities we analysed as a function of the S² parameter, with increasing mixing, i.e. decreasing S² towards the right. We find the exact overall behaviour when comparing each parameter with the dynamical time. As the clusters become more mixed (smaller value of S²), the observational parameters approach zero. However, we also distinguish additional features due to the cluster’s different evolutionary state. Panel a of Fig. 12 shows the effects of core collapse and the subsequent divergence into different evolutionary paths. This divergence is due to the degree of spatial mixing in the core at core collapse. The clusters that do not undergo core collapse do not show a temporary increase in S² (corresponding to the bump shown in Fig. 2). In the case of the velocity anisotropy area difference in panel b, we see similar behaviour as in Fig. 8. Lastly, for the difference in the projected angular momentum in panel c, all the models begin to converge onto the same path after attaining a mixing level of S² ≃ 9. It is important to note, however, that at this value of S² ≃ 9, differences between models persist when examining the velocity anisotropy area difference.

Fig. 12 Comparison of the three projected parameters that quantify the degree of mixing of the multiple populations with the S2 parameter from the phase-space mixing. The S2 parameter decreases to the right following the increasing phase-space mixing of the two populations, as stated at the top of the panel a.

5 Discussion and summary

In this work, we analysed the long-term evolution of simulations with multiple populations to characterise the dynamical mixing process. Our simulations follow the evolution of models with an initially concentrated and underfilling 2P in an extended, more massive, and tidally filling 1P. We explored different initial conditions to examine differences in the concentration of the 2P, the primordial velocity anisotropy of the cluster, and the strength of the tidal field.

Multiple studies have shown that the different stellar populations in GCs have a large variety of spatial and kinematic mixing as traced by differences in concentration (see, e.g., Dalessandro et al. 2019; Leitinger et al. 2023; Cadelano et al. 2024), rotation (see, e.g., Cordero et al. 2017; Cordoni et al. 2020; Dalessandro et al. 2021, 2024; Martens et al. 2023; Leitinger et al. 2025), and velocity anisotropy (see, e.g., Richer et al. 2013; Bellini et al. 2015; Cordoni et al. 2020, 2025; Libralato et al. 2023; Dalessandro et al. 2024). To trace the evolution of the mixing processes, we followed the differences in the phase-space distribution of each population in our models. We defined the S² parameter as the sum of the energy and angular momentum distribution square differences over the whole phase space. The S² parameter allows us to build a comprehensive dynamical picture of the mixing process, including both the spatial and kinematic components of the populations mixing. We find that (as expected) mixing increases with dynamical age (see Fig. 2). We also show that core collapse and primordial velocity anisotropy can alter the mixing rate of the cluster, and core collapse may even cause a temporary slight inversion of the phase-space mixing (see Figs. 3 and 4). As shown in Fig. 2, models can achieve similar values of S² at different physical and dynamical ages, either due to variations in their mixing rates or initial conditions. Our models reach different degrees of dynamical mixing at 12 Gyr; among all the models presented, the N05M and N2M exhibit, respectively, the most and the least mixing.

We followed the evolution of three “observable” parameters that can trace the degree of spatial and kinematic mixing and found that they display the same general trends as in previous observational studies (see, e.g., Dalessandro et al. 2019, 2024; Cadelano et al. 2024), where the differences in the spatial and kinematic properties of multiple populations decrease with dynamical age. We also find that in our models, the strength of the initial 1P and 2P differences already decreases in the clusters’ early evolutionary phases when the clusters are still dynamically young (n_h ~ 1–2) (see also Vesperini et al. 2021; Sollima 2021; Cadelano et al. 2024). Present-day differences still contain important information about the 1P and 2P dynamical properties imprinted at the time of cluster formation. However, it is important to note that the observed present-day strength of these differences can be much smaller than the initial ones, even in dynamically young clusters, and that different initial conditions may lead to similar present-day properties. Furthermore, and as a cautionary note, our data do not consider limitations typically affecting observational analyses, such as incompleteness or uncertainties in population membership, which can add stochasticity to the evolutionary paths and reduce the detected spatial differences shown in Fig. 5.

We quantified the difference in the velocity anisotropy profiles of each population with the ΔA_β parameter, which measures the area difference between the velocity anisotropy profiles. This parameter is the most lasting kinematic signature for mixing, as the models still keep differences after multiple relaxation times. The overall trend is driven by a radially anisotropic 2P and a tangentially anisotropic 1P. For all the initially isotropic models, the 1P remains isotropic in the inner regions and develops tangential anisotropy in the outer regions, while 2P evolves towards a radially anisotropic velocity distribution with the strength of radial anisotropy depending on the initial concentration of 2P: initially less concentrated 2P spatial distributions results into a shallower radial velocity anisotropy (see Fig. 8). The model with an initial radially anisotropic velocity distribution illustrates the complexity of the interplay between the effects due to different spatial and velocity distributions. In the initially anisotropic model we studied (N1M-ra with r_h,1P/r_h2,P = 20), the difference between 1P and 2P anisotropy, after n_h ~ 1, evolves similarly to that of an initially isotropic model but with smaller initial spatial differences (N1M-r10). Additional models are necessary to further explore these aspects of the clusters’ dynamical mixing.

Finally, we studied the evolution of σ_Lz, the dispersion of the angular momentum calculated from the tangential component of the proper motion velocity on the plane of the sky and the projected radial distance from the cluster’s centre. We show that this parameter is linked to the degree of radial anisotropy, and that the differences between the values of this parameter for 1P and 2P, Δσ_Lz, trace the differences in the degree of anisotropy of their velocity distributions. More importantly, this parameter allows for a more direct way to quantify differences in the velocity distribution of the two populations and provides another observational measure of the kinematics of multiple populations, as well as the degree of kinematic mixing. We explored its evolutionary path (see Figs. 10 and 11) as well as its co-evolution with global phase-space mixing measured by S² (see Fig. 12). For dynamical ages n_h ≳ 1–2, the evolutionary paths of Δσ_Lz for models with different initial conditions approach a common universal track. As pointed above, dynamical differences observed in old clusters provide key information on the differences between the 1P and 2P spatial and kinematic properties set by formation processes and their evolution. However, even for dynamically young clusters, it is necessary to exercise caution in interpreting the specific strength of present-day differences as a close indication of those imprinted at the time of multiple population formation.

Models extending the investigation presented here to a broader set of initial conditions are necessary to further explore the evolutionary path of the parameters introduced in this paper, as well as their link to the clusters’ initial conditions. Including other dynamical indicators, such as differences in the degree of energy equipartition and internal rotation, will provide additional key steps towards a comprehensive picture of the dynamical mixing of multiple populations in GCs.

Acknowledgements

EV and FA acknowledge support from STScI grant AR-16157. EV acknowledges also support from the John and A-Lan Reynolds Faculty Research Fund. ED acknowledges financial support from the Fulbright Visiting Scholar program 2023. ED is also grateful for the warm hospitality of Indiana University, where part of this work was performed. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute. FA also acknowledges the following Python packages used during the research presented here: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), emcee (Foreman-Mackey et al. 2013) and Astropy (Astropy Collaboration 2022), while all figures were made using Matplotlib (Hunter 2007).

Appendix A Simulation properties at 7 Gyr and 12 Gyr

Appendix B Additional figures

We include in this appendix three additional figures to complement Figures 1 and 2. Figures B.1 and B.2 show the phase-space distribution for all models at t = 0, 4, 8, 12 Gyr, while Fig. B.3 shows the S² parameter (as in Fig. 2) and the corresponding marginalised version for angular momentum S²(L) and energy S²(E).

Fig. B.1 Phase-space distribution evolution for models N2M, N1M, N05M and N05M-wf at t = 0 Gyr, t = 4 Gyr, t = 7 Gyr and t = 12 Gyr. For each model, we include the corresponding dynamical age traced by the number of half-mass relaxation times nh.

Fig. B.2 As in Fig. B.1, but for models N1M-ra, N2M-r10 and N1M-r10.

Fig. B.3 Time evolution of the S2 parameter for all models as in panel (a) in Fig. 2, but separated and also including the marginalisation over angular momentum S2(L) and energy S2(E). We have marked a vertical line at the time of the first core collapse for all models that undergo core collapse within the 13 Gyr of evolution.

References

All Tables

All Figures

Fig. 1 Phase-space distribution of stars in the 1Ps and 2Ps for the models N2M and N05M at different times. At t = 0 Gyr, both clusters show a similar relative distribution of stars in the 1Ps (red) and 2Ps (blue), due to their similar half-mass radius and mass ratios. As the clusters evolve, the 2P expands towards the regions initially dominated by the 1P. By t = 12 Gyr, the N2M model is still partially mixed, while the N05M model has a significant degree of mixing. The insets show the marginalised distribution of energy and angular momentum, likewise highlighting the evolution of dynamical mixing. In all panels, the dashed lines follow the maximum allowed angular momentum at a given energy, defined by the angular momentum of the circular orbit of the corresponding energy.

In the text

Fig. 2 Evolution of the S2 parameter for all models. Panel (a) shows the evolution of S2 with respect to the physical time in the simulations. All models show different degrees of mixing at most times. At t = 12 Gyr (vertical line), model N2M is the least mixed, while model N05M is the most mixed in phase space. Panel (b) shows the evolution of S2 as a function of the number of half-mass relaxation times. The symbols show the value of S2 at 12 Gyr and the respective number of relaxation times.

In the text

Fig. 3 Effects of partial mixing during core collapse. Panels a and b: phase-space distribution of the models N2M and N05M at ~500 Myr before core collapse. The N2M model is far from fully mixed, while the N05M model has a significant degree of mixing. Panel c: evolution of the S2 parameter centred at the time of core collapse for the models N2M and N05M (tcc = 7.06 Gyr and tcc = 8.62 Gyr). The markers indicate the corresponding time of phase-space distributions in panels a and b. Paneld shows the value of the S2 parameter only for the energy component of the phase space. The peaks in S2(E) show the effects of core collapse and subsequent core oscillations for the model N2M, where the 2P dominates the inner regions. While for N05M model, the two populations are almost completely mixed at the time of core collapse, and the evolution of S2(E) does not show any feature associated with the core collapse.

In the text

	Fig. 4 Evolution of the S2 parameter (a), as well as S2 measured over the angular momentum (S2(L), panel b) and energy (S2(E), panel c) spaces for the N1M and N1M-ra models. Primordial radial anisotropy increases the population mixing rate while keeping the cluster dynamically younger. The model N1M-ra does not undergo core collapse within its 13 Gyr of evolution.
In the text

Fig. 5 $A_2^{+}$ Parameter evolution for all GCs as a function of time (a) and number of relaxation times (b). All models follow a similar evolutionary path as the two populations spatially mix. The panel a scatter reduces when we trace the cluster dynamical ages. After one relaxation time, the initial difference between the rh ratio of 20 (N2M, N1M, N1M-ra, N05M, and N05M-wf) and rh ratio of 10 models (N2M-r10 and N1M-r10) is indistinguishable.

In the text

Fig. 6 Area of the cumulative stellar distribution within 0.5Rh (ϕ(<0.5Rh)) for the two populations in models: (a) N2M and (b) N05M. As in Fig. 3, we centre the time evolution at the respective time of the first core collapse. The two populations in model N2M have not achieved complete spatial mixing and ϕ2P(R < 0.5Rh) > ϕ1P(R < 0.5Rh). The model N05M has already attained a significant spatial mixing and ϕ2P(R < 0.5Rh) ~ ϕ1P(R < 0.5Rh). Furthermore, for model N05M, both populations show the same rate of increased concentration at core collapse, which explains the lack of signatures in S2(E) (see Fig. 3).

In the text

Fig. 7 Velocity anisotropy area (Aβ) and area difference (ΔAβ) for the N2M model. Panels a and b show the velocity anisotropy for the two populations (red upper triangles, PO1P, and downwards blue triangles, PO2P) in the N2M model at 4 Gyr and 12 Gyr, respectively. The hatched regions represent the respective anisotropy area for each population from the best-fit velocity dispersion and anisotropy parameterisation (see main text). Panel c shows the evolution of the anisotropy areas for the two populations over time; as the 2P becomes less radially anisotropic, the 1P becomes tangentially anisotropic and keeps a constant difference with the 2P, as shown in panel d, where we show the evolution of the anisotropy area difference. The vertical lines in panels c and d represent the times for which we show the velocity anisotropy profiles in panels a and b.

In the text

	Fig. 8 Anisotropy area (Aβ) and area difference (ΔAβ) as described in Fig. 7, for all models in our sample with respect to their dynamical ages. For better visualisation, we separated nh into a linear section (nh ≤ 6) and a logarithmic section (nh > 6).
In the text

Fig. 9 Projected angular momentum dispersion σLz profiles for King models with different velocity anisotropy profiles (panel a) and their respective velocity anisotropy profiles (panel b). As the King models become more radially anisotropic (smaller values of ra/rh), the projected angular momentum dispersion decreases. The insets on panel a show the distribution of the projected angular momentum at r ~ 2Rh for the models with the most and the least radial anisotropy. The grey area in panel a shows the range of radial distances of stars used to calculate the histograms in the panel. We normalised σLz by the half-mass angular momentum Lh, corresponding to the angular momentum of a circular orbit of an equivalent point-mass potential (see Sect. 3).

In the text

Fig. 10 Projected angular momentum dispersion σLz profiles for both populations in the N2M model (panel a) and the N05M model (panel b) at 12 Gyr. The vertical grey solid lines in panels a and b represent the difference in σLz between the populations at R = 2Rh. Panel c shows the time evolution of σLz difference at 2Rh; the closer this value is to zero, the closer the two populations are to kinematically mixing.

In the text

	Fig. 11 Evolution of the projected angular momentum difference at 2Rh for all models in our sample as described in panel c of Fig. 10. In panel a, this is with respect to the physical time of the simulation, and in panel b this is with respect to their dynamical age. Once we consider the dynamical age, all models converge into the same evolutionary path, independently of their initial conditions.
In the text

	Fig. 12 Comparison of the three projected parameters that quantify the degree of mixing of the multiple populations with the S2 parameter from the phase-space mixing. The S2 parameter decreases to the right following the increasing phase-space mixing of the two populations, as stated at the top of the panel a.
In the text

	Fig. B.1 Phase-space distribution evolution for models N2M, N1M, N05M and N05M-wf at t = 0 Gyr, t = 4 Gyr, t = 7 Gyr and t = 12 Gyr. For each model, we include the corresponding dynamical age traced by the number of half-mass relaxation times nh.
In the text

	Fig. B.2 As in Fig. B.1, but for models N1M-ra, N2M-r10 and N1M-r10.
In the text

	Fig. B.3 Time evolution of the S2 parameter for all models as in panel (a) in Fig. 2, but separated and also including the marginalisation over angular momentum S2(L) and energy S2(E). We have marked a vertical line at the time of the first core collapse for all models that undergo core collapse within the 13 Gyr of evolution.
In the text