© The Authors 2026

1 Introduction

The globular cluster (GC) system of the Milky Way (MW) includes more than 160 observable objects. One of the most complete and recent information about them was provided by the catalogue of Baumgardt & Vasiliev (2021)¹, which contains the structural parameters for 168 objects identified as GCs. The masses of most GCs lie in the range of ~10⁴–10⁶ M_⊙, and their half-mass radii (r_hm) fall within 2–20 pc. Since the ages of GCs are comparable to the age of the Galaxy itself, they are definitely one of the best indicators in its dynamical evolution, including the formation of various structures and the whole assembly history (Antonini et al. 2012; Kruijssen et al. 2020). Some of the GCs have a pronounced heterogeneous stellar population, that is, multiple stellar populations (Bica et al. 2016). This might indicate not only the self-enrichment processes, but also the intense interactions with a surrounding environment. One popular explanation for the presence of multiple stellar populations in Galactic GCs (see Bastian & Lardo 2018; Gratton et al. 2019, and references therein) is the external origin of chemically distinct gas or stellar material, which may also result from close interactions caused by potential mergers of individual clusters (Ferraro et al. 2009; Origlia et al. 2011; Khoperskov et al. 2018; Mastrobuono-Battisti et al. 2019). Direct GC close encounters or even collisions might seem extremely rare. However, according to several studies, such interactions might be fairly probable for GCs in the central MW region. For instance, Miocchi et al. (2006) modelled the evolution of massive and compact GCs in a triaxial gravitational potential and showed that although the impact on internal structure and orbital evolution from even face-on GC collisions is negligible, GC merging events into nuclear star cluster-like systems are potentially possible due to dynamical friction in the very inner regions of the Galaxy (Capuzzo-Dolcetta & Miocchi 2008a,b; Kuvatova et al. 2026). According to Khoperskov et al. (2018) and Mastrobuono-Battisti et al. (2019), in a system of about 100 massive (~10⁷ M_⊙) GCs, a few physical collisions can occur for each billion years, along with several close passages even with possible mass exchange. Their simulations included cases of GC mergers into a single system within 0.5 Gyr. In the long term, the resulting clusters lose mass and become more compact. When the merging clusters have similar dynamical properties, their stellar populations mix rapidly; otherwise, clusters with differing densities produce a system with multiple populations and distinct density profiles.

In this context, numerical studies of GC populations indicate that modelling clusters in a joint system leads to structural differences compared to their isolated evolution, primarily through tidal interactions and repeated close passages (Khoperskov et al. 2018). However, the strength of this effect is highly sensitive to the spatial configuration of the system and is expected to be most significant in dense environments or during the early evolutionary stages (Miocchi et al. 2006). Combes et al. (1999) further argued that clusters undergoing substantial mass loss can develop elongated shapes by precessing around the z-axis. While this does not constitute direct evidence that simultaneous modelling of multiple GCs necessarily leads to deviations from spherical symmetry, it suggests a possible link: if interactions within a GC system enhance the mass-loss efficiency, they might indirectly promote morphological transformations of this type.

In our previous studies (Ishchenko et al. 2023a,b, hereafter, Paper I, and Paper III), the orbital modelling of GCs backward in time for 10 Gyr showed that some GCs complete a significant number of revolutions around the Galactic centre (e.g. ~620 for Terzan 4 and ~540 for NGC 6624). Several such clusters are associated with regions such as the central bulge and disk. The long-term presence of several objects within a relatively small Galactic centre volume can naturally lead to a certain number of mutual encounters. In our study in Paper III, we numerically modelled 147 GCs in five selected time-variable MW-like potentials (TVPs) extracted from the cosmological database IllustrisTNG-100 (Nelson et al. 2018) with the aim to conduct a statistical analysis of the likelihood of close encounters and possible collisions among GCs. It was found that most of these events occurred in the central region of the Galactic disk, within ~4 kpc of the centre, with two ring-like zones with an increased encounter density observed at radii of about ~1 and ~2 kpc of the Galactic centre (see Fig. 8 in Paper III). Across all modelled potentials, the rate of close passages (≤50 pc) was estimated as ~10 events every billion years. Moreover, actual collision events were detected for 22 GC pairs with a very high likelihood. Some of the most active participants in these events are Terzan 4, Terzan 2, NGC 6624, NGC 6440, Terzan 5, and Liller 1, all of which belong to the MW central GC subsystem and are located near the Galactic centre (see the detailed Table 3 in Paper III). For a more detailed analysis of collisions and close passages, we here selected the Terzan cluster family: Terzan 2 and Terzan 4 as the most active collisional pair in all our five MW-like potentials. Terzan 5 was added to the list because of its high probability of close passages with Terzan 4 and also with Terzan 2. Terzan 5 is also one of the most massive clusters known today, with ~10⁶ M_⊙. These GCs were directly modelled as a full N-body system with up to several million particles each and with a detailed stellar evolution.

Terzan 2 is a relatively metal-rich GC ([Fe/H] = −0.86), while Terzan 4 is moderately metal poor ([Fe/H] = −1.38). Both belong to the inner bulge population (Schiavon et al. 2024) and exhibit evidence of multiple stellar populations (Kunder et al. 2021). However, the most prominent example of a chemically complex stellar system is Terzan 5. This GC shows a wide spread in iron abundance (~1 dex) and a significant stellar age spread of approximately 5–8 Gyr. Its metallicity distribution is clearly multi-modal (Massari et al. 2014), with three distinct iron peaks: an extremely high-metallicity population (29%, [Fe/H] ≈ 0.25), a moderately metal-rich population (62%, [Fe/H] ≈ −0.3), and a relatively metal-poor population (6%, [Fe/H] ≈ −0.8). This complexity indicates a highly inhomogeneous stellar population.

Several studies have questioned whether Terzan 5 is a genuine GC at all (Origlia et al. 2011; Bailin & von Klar 2022). Its orbit and chemical properties clearly suggest an in situ origin within the Galactic bulge, however (Origlia et al. 2025). The observed population complexity might partially be attributed to self-enrichment processes and the formation of subsequent stellar generations (Ferraro et al. 2009; Valcarce & Catelan 2011). However, the predicted mass fraction of second-generation stars is insufficient to explain the observed peculiar features. Alternative scenarios invoke a more complex dynamical history. For example, the presence of multiple populations could result from the merger of two originally distinct stellar systems (Gavagnin et al. 2016; Khoperskov et al. 2018; Mastrobuono-Battisti et al. 2019; Pfeffer et al. 2021).

Based on our previous findings in Paper III, we carried out a more detailed dynamical analysis here, including a separate and common orbital-evolution modelling of the Terzan cluster family. We analysed the individual and collective mass loss of the clusters and as the energy evolution of their mutual gravitational interaction.

In accordance with this, the paper is organised as follows. In Sect. 2, we introduce our N-body reconstruction of the chosen GCs. In Sect. 3, we estimate the mass loss of the GCs during their evolution. Section 4 analyses the GCs collision probability on cosmological timescales. In Sect. 5, we present additional individual simulations to investigate the mutual gravitational effect of GCs. Finally, in Sect. 6 we summarise our results and draw the main conclusions.

2 Full N-body reconstruction of Terzan 2, 4, and 5

2.1 Time-variable potential implementation

To enhance the physical realism of our simulations, we employed data from the cosmological magnetohydrodynamical IllustrisTNG-100 simulation (Nelson et al. 2018). One of the key advantages of this cosmological simulation database is that the time-dependent gravitational potentials for a large number of different individual galaxies can be extracted.

From the entire set of simulated galaxies in TNG, we selected the MW-like object with ID 411321, which at zero redshift (z = 0) exhibits properties that most closely match the observed parameters of today’s MW. In particular, the selection was based on the similarity in total mass, as well as the spatial scales of the disk and halo components. The sampling and fitting procedures for this selected potential and for other MW-like potentials are discussed in detail in Table 1 and Fig. 2 of Paper I. The code routines for sampling and fitting the selected potentials are also publicly available on GitHub².

For the full N-body integration of the GCs, we used the popular parallel high-order code φ–GPU³, which employs a fourth-order Hermite integration scheme with hierarchical individual block time-steps. This direct body code has proven itself well in modelling GCs in our previous papers (Ishchenko et al. 2024; Weis et al. 2025; Kuvatova et al. 2026). An external potential 411321 was especially incorporated into the code. In addition, a stellar evolution library was implemented (Banerjee et al. 2020; Kamlah et al. 2022a,b), enabling the treatment of the main stellar evolution processes during the dynamical modelling.

Fig. 1 Simultaneous orbital evolution of the Terzan clusters in the TVP 411321. The orbital evolution is present in three planes, (X–Y), (X–Z), and (R–Z), where R is the planar Galactocentric radius. The total integration time is an 8 Gyr look-back time. The black marks show the current positions, and empty marks show the positions at −8 Gyr.

2.2 Initial models for Terzan 2, 4, and 5

As mentioned in Sect. 1, we selected three Terzan clusters for our study, namely Ter2, Ter4, and Ter5. These clusters exhibit a high probability of mutual interactions over their dynamical orbital histories. To reconstruct the evolutionary history of the mass and spatial parameters of the selected clusters at least approximately, we adopted the following initial models for each cluster individually.

We assumed that the stellar mass distribution in our models spanned the range 0.08–100 M_⊙, following the Kroupa initial mass function (Kroupa 2001). We also assumed for each cluster its own metallicity, which affects the stellar mass distribution and subsequent evolution. The models were initially brought into dynamical equilibrium assuming the King models (King 1966), which are characterised by the half-mass radius r_hm and the central concentration parameter W₀. These three parameters (initial M, r_hm, and W₀) were selected so that after 8 billion years of stellar and orbital evolution, the cluster system would have an M and r_hm comparable to observational data (see Sect. 3).

Each particle in the simulation represented ten real stars of the same type, which allowed us to fit the cluster models more efficiently, particularly in terms of total mass and half-mass radius, to match their observed present-day values. This approximation in the particle number does not affect the system two-body relaxation time significantly. A detailed discussion of the effect of this reduction factor can be found in Panamarev et al. (2019) and in Ishchenko et al. (2024). We refer to this setup as the common 3xGC run, as all three Terzan clusters are modelled simultaneously.

Table 1 summarises these physical parameters of the clusters at T = −8 Gyr. The initial positions and velocities for GC centres of mass used in the simulations correspond to their calculated stages 8 Gyr ago, obtained by backward integration of each cluster (treated as a point mass; these results were presented in detail in our Paper I). The orbital evolution of the selected Terzan clusters in 411321 time-variable external potential is shown in Fig. 1. We also indicate the position for each GCs at −8 Gyr as empty marks in the plots. The effect of dynamical friction on the GC orbital evolution in the context of the point-mass approximation is discussed in Appendix A.

The numerical experiments were aimed to obtain well-fitted individual models for each GC by varying key parameters, including the total mass, half-mass radius, and King concentration parameter. The ref. runs, in which full N-body simulations of each individual Terzan clusters were carried out separately, were conducted on the computer cluster of the Fesenkov Astrophysical Institute in Almaty, Kazakhstan and as on the LOTR cluster at the Main Astronomical Observatory of the National Academy of Sciences of Ukraine. Each run used a pair of NVIDIA RTX 4080 or RTX 4090 GPUs, achieving an average sustained performance of ≈21.2 Tflops per card. The common 3xGC run was executed on the JUWELS Booster supercomputer at the Jülich Supercomputing Centre under the madnuc project, using four nodes, each with a 4×NVIDIA A100 GPU. The total modelling time took ≈5 months for the largest common 3xGC run.

Fig. 2 Evolution of Mtid and rhm for the Terzan clusters due to mass loss based on common 3xGC run (coloured lines) and ref. run (yellow lines). The dotted lines show the estimated current mass, rhm, and the grey zone shows the ±2σ observational uncertainty based on the catalogue by Baumgardt & Vasiliev (2021) for today.

3 Mass loss

The long-term mass evolution of a GC is shaped by internal and external galactic dynamical processes. One significant source of cluster mass loss (especially in the first few hundred million years of evolution) is high-mass stellar evolution. As a result of the high-mass stars stellar evolution, we also observe the kick processes of black holes and neutron stars from the cluster. Clusters on eccentric orbits can experience stronger tidal forces, especially near their pericentre, which also enhances the stellar stripping and accelerates the subsequent mass loss (Spitzer 1987; Webb et al. 2014). The evolving galactic potential can further amplify these effects by driving changes in orbital energy and angular momentum over time. The stripped stars often form extended tidal tails that trace the cluster orbit and even provide observable signatures of its dynamical history (Habibi et al. 2013; Park et al. 2018; Arnold & Baumgardt 2025). This is particularly relevant for central MW region clusters such as Ter5, Ter4, and Ter2, which follow compact eccentric orbits with pericentric distances of ≲0.5 kpc. Their proximity to the Galactic centre and repeated passages through these dense regions make them highly susceptible to tidal stripping, which offers valuable insight into mass loss and the formation of tidal debris.

In Fig. 2, we present the evolution of the mass within a tidal radius and of the half-mass radii for Terzan clusters in more detail. The cluster tidal radius (also known as Jacobi or King radius) was calculated based on the numerical iteration of the M_tid and r_tid values according to the King equation, $$r_{\mathrm{t}\mathrm{i}\mathrm{d}}={\left[\frac{G\cdot M_{\mathrm{t}\mathrm{i}\mathrm{d}}}{4{\mathrm{\Omega }}^2-{\kappa }^2}\right]}^{1/3},$$(1)

where G, M_tid, Ω, and κ are the gravitational constant, the tidal cluster mass, the circular and the epicyclic frequencies, respectively, of a near-circular orbit in the Galaxy at the GC current location (for more details, see King 1962 and Ernst et al. 2011).

The grey bands in Fig. 2 represent the ±2σ uncertainties of the present-day tidal mass and half-mass radius estimates for the selected GCs, taken from the Baumgardt & Vasiliev (2021) catalogue. For Ter2 and Ter4, the model uncertainties for the final tidal mass are ≈29%, whereas for Ter5, the uncertainty is only ~7%. This indicates that in our constructed initial models, the mass of Ter5 is determined most precisely. On the other hand, the masses of Ter2 and Ter4 are less well constrained.

The tidal stripping experienced by the GCs is also expected to lead to the formation of extended tidal tails. These tidal features are clearly present in our simulations. However, their detection in observations remains challenging, as their typical volume densities are very low (≈0.01 M_⊙/pc³), significantly below the average stellar mass density in the central Galaxy regions (Gregersen 2009). As a result, their identification likely requires full 6D phase-space information, which may become available with future data releases from the Gaia space mission.

The plots show that Ter 5 lost about 45% of its initial mass. For Ter2 and Ter4, the mass-loss evolution was quite similar, but by the present day, they already lost almost 75% of their initial masses. To summarise, Ter5 loses no more than half of its initial mass, in contrast to the other two clusters, which undergo significantly stronger mass loss. Ter5 clearly exhibits the smallest fluctuations in its parameters throughout the evolution and provides the best match to its present-day properties. In contrast, the values of the parameters for Ter2 and Ter4 at the end of the simulation are systematically higher than the observed values and even lie outside their respective ±2σ uncertainty ranges. In addition, the initial models we assumed for the Terzan cluster family members are probably good enough for the purpose of studying their long-term mutual gravitational interactions.

As initial conditions for these ref. runs for each GC, we used exactly the same particle data as listed in Table 1. Using these reference runs, we were able to examine the individual differences in the dynamical evolution of each Terzan GC. In Fig. 2, we present the dynamical mass-loss evolution for each cluster for the ref. runs and the common 3xGC run. The plots show that in general, the particle and mass loss, together with the half-mass radius evolution, are similar for all three systems, with some differences, especially for the lower-mass Ter2 and Ter4. In the case of the more massive Ter5, the common 3xGC run shows practically the same mass-loss evolution as in the Ter5 ref. run.

To understand this mass-loss behaviour in different GCs on the cosmological timescale, we compared the mutual gravitational effect of the clusters in the common 3xGC run. The most massive Ter5 clearly has the most significant gravitational perturbation acting on the other two Terzan systems. The effect of the low-mass Ter2 and Ter4 systems on Ter5 is significantly weaker. This effect is clear in the mass-loss and half-mass radius plots in Fig. 2. This plot shows the combined effect of global gravitational perturbations, which act constantly, together with collisions (close passages) between individual Terzan systems. The effect of the global gravitational perturbation is visible almost immediately, after a few hundred million years from the beginning of the simulation, in the higher mass-loss rate compared with the corresponding ref. runs. By the end of the simulation, the mass loss for the individual ref. runs and for the common 3xGC model becomes very similar. This effect is naturally stronger for Ter2 and Ter4.

We refer to Appendix C for details about the internal dynamic distribution of massive star remnants for Terzan clusters after an evolution of 8 Gyr. We present the analysis as from common 3xGC run and as ref. run in this appendix.

4 Terzan 2, 4, and 5 collision tracking on a cosmological timescale

Throughout the simulation, we monitored the possible close approaches between GCs at fixed distances shorter than dR = 100 pc. In Table 2 we list the GC pairs involved in these events, along with the time of the encounter. We also present the minimum separation and the relative velocity between their centres of mass together with the corresponding values of their tidal radii at that moment.

The most significant of the recorded encounters are four close collisions between Ter2 and Ter4 and two collisions between Ter4 and Ter5. In these cases, the sum of the tidal cluster radii exceeds their minimum separation. A particularly notable event is the encounter between Ter2 and Ter4, during which the clusters approached to within 10 pc, while their tidal radii were 12 and 17 pc, respectively. In this case, however, the relative fly-by velocity between clusters is also one of the highest for this pair (~320 km/s). As expected, Ter4 remains the most dynamically active cluster in terms of close encounters and interactions with the other two systems.

In Fig. 3, we present the evolution of the semi-major axis, a, of each cluster over a simulation of 8 Gyr and mark the collision events (according to Table 2). The encounters are distributed relatively uniformly in time and in semi-major axis, indicating no clear dependence of the close interaction probability on orbital phases.

The right panel of Fig. 4 shows the collision pairs from Table 2 in a 3D space projection, with a colour-coding by time. The coloured spheres represent the doubled tidal radii of the two clusters at the moment of the encounter. The left panel displays the position of cluster pairs that approached each other within 200 pc. A sharp increase in the number of such pairs is evident for dR = 200 pc. In the X − Y projection, they form a ring-like structure around the Galactic centre. This effect was already observed in Paper III (see Fig. 8).

Fig. 3 Evolution of the GC semi-major axes during the simulation with marked collision events in the common 3xGC run.

Fig. 4 Collisions of the Terzan clusters in 3D projections for different dR conditions. Left: dR = 200 pc. Right: dR = 100 pc.

Fig. 5 Global potential energy evolution of the cluster systems for the common 3xGC run and the individual ref. runs. Left: full-time evolution of the models up to 8 Gyr. Right: same energy evolution with a very detailed output near the close passage marked A in Table 2.

5 Gravitational potential perturbations caused by the mutual global effect of Terzan 2, 4, and 5

Based on our previous analysis of the common 3xGC run (three GCs modelled together), we already detected some common effects of three separate GCs in the time-evolved external gravitational field of the Galaxy. For a detailed analysis of this common gravitational effect of the GCs on each other, we decided to carry out a new set of separate N-body runs in which we modelled the dynamical evolution of each Terzan GC individually in the same external MW-like potential.

To quantitatively check this effect, we analysed the behaviour of the energy outputs for the ref. runs and the common 3xGC runin high detail. We present these results in the left panel of Fig. 5. This plot shows the total potential energy of the systems, namely the sum of particle self-gravity and their energy in the external gravitational field. In the common 3xGC run results (black line), the energy perturbation level is significant because the three GCs collectively affect each other gravitationally. This perturbation is present throughout the entire simulation time. To obtain a more detailed view of this cyclic perturbation, we present the very detailed time output of these energies during one of the Ter2 − Ter4 close-passage events (marked A in Table 2) in the right panel of Fig. 5. In addition to the higher potential energy level in the case of the common 3xGC run, we also detect low-level short-timescale periodic perturbations of ~200 Myr that reflect the mutual interactions between the three GCs.

To demonstrate the collective behaviour of the GC system in the case of the common 3xGC run, we also analysed the ellipsoidal shape of the GC main stellar body in detail. To determine the semi-major axis, we used the well-tested algorithm triax. This method is described in detail in Appendix C in Regály et al. (2023). The triax method is based on the well-known Jacobi eigenvalue matrix algorithm⁴. As a result of the method, we obtained the ratios of the semi-major axis b/a and c/a and as a final parameter: the evolution of the triaxiality, defined as T = (a² − b²)/(a² − c²).

Figure 6 shows the time evolution of the cluster triax as the T parameter for common 3xGC run and ref. runs for highly detailed outputs for the A case of closed passage. In the top panel, we present the time evolution of the mutual distances between the pairs of clusters, which is calculated as $${\mathrm{D}\mathrm{R}}_{\mathrm{D}\mathrm{C}}=\sqrt{{(X_{\mathrm{G}\mathrm{C},\mathrm{i}}-X_{\mathrm{G}\mathrm{C},\mathrm{j}})}^2+{(Y_{\mathrm{G}\mathrm{C},\mathrm{i}}-Y_{\mathrm{G}\mathrm{C},\mathrm{j}})}^2+{(Z_{\mathrm{G}\mathrm{C},\mathrm{i}}-Z_{\mathrm{G}\mathrm{C},\mathrm{j}})}^2}.$$

In the middle panel of Fig. 6 we present the triaxial parameter definition for the common 3xGC run. In the bottom plot, we present the same triax evolution for ref. runs. The figures clearly show the differences between the two types of runs. The ref. runs evolve smoothly for the selected time interval, without significant changes in the triax parameter at one level for Ter2 in 0.78, for Ter4 in 0.74, and for Ter5 in 0.60. For the common 3xGC run, the behaviour is completely different: there are strong and short fluctuations of T caused by the changing mutual separation between the clusters. The overall variability of T is very strongly correlated with the DR_DC in the top panel. In Fig. B.1, we also provide the evolution b/a and c/a axes separately for the common 3xGC run and ref. runs.

Fig. 6 Time evolution of the mutual distance between GC pairs and the axial ratios for the common 3xGC run and ref. runs. The middle and bottom panels show the triaxiality. The top and middle panels correspond to the common 3xGC run, and the bottom panel shows the ref. runs for the Terzan clusters.

6 Conclusions

The evolution of our GC subsystem remains one of the most actively discussed topics in contemporary Galactic astrophysics. GCs provide crucial insights into the early stages of Galaxy formation, including past accretion events and the subsequent dynamical evolution of the system.

After performing a set of high-resolution dynamical N-body simulations of the three Galactic centre GCs Ter2, Ter4, and Ter5, we analysed the individual and mutual evolution of these systems in detail and identified a list of pairwise collisions (close encounters; see Table 2). These collisions persisted during the whole modelling time of 8 Gyr. In a few cases, the separation between the cluster pairs became smaller than the sum of their r_tid. During these events, the relative velocities of the clusters were quite high, ~200–300 km/s. These high velocities prevent the possible formation of gravitationally bound Terzan cluster pairs.

An interesting finding of our simulations is the slightly different mass-loss behaviour. The individual single Terzan full N-body simulations (ref. runs) show a significantly lower initial mass loss (particle loss) during the first few hundred million years compared to the common 3xGC model, in which all three clusters were simulated simultaneously (see Fig. 2). At the same time, the potential energy differences between the common 3xGC and ref. runs were significant (see Fig. 5).

The collective behaviour of the GC system in the common 3xGC run was analysed, revealing a clear effect of mutual interactions on the evolution of the cluster shapes. The application of the triax algorithm demonstrated strong correlations between the axis ratios and the mutual distances of the cluster pairs, especially for Ter2, which evolved from nearly spherical to distinctly prolate (b/a and c/a ~ 0.75–0.80), and for Ter5, which showed similar but less pronounced changes (see Fig. 6). In contrast, Ter4 remained relatively stable. In the ref. run, however, the clusters preserved almost spherical shapes (with axis ratios ≳0.95), highlighting that the triaxial deformations observed in the common 3xGC case are driven by mutual gravitational interactions among the modelled Terzan clusters.

The high-mass stellar remnant normalised cumulative distributions for these two groups of runs also showed a more compact number distribution result for the common 3xGC run today. This effect is especially clearly manifested for the low-mass systems Ter2 and Ter4 (see Appendix C).

We conclude based on these findings that the full group of GCs within this zone of a few kiloparsec should preferably be included in the modelling during the dynamical evolution of Galactic central-region clusters. This complex modelling is especially important in the case of low-mass systems because their high-mass neighbours, for example, Ter5 or Liller 1, which have current masses of more than 10⁶ M_⊙, can affect their internal dynamical and shape behaviour significantly (see Figs. 2 and 6).

Deviations in the GC shapes from spherical symmetry are commonly attributed either to internal mechanisms, such as rotation, or to external effects, primarily tidal forces exerted by the surrounding environment. In support of the latter, previous studies have reported systematic elongation of clusters located near the Galactic bulge, with their major axes preferentially aligned towards the Galactic centre (Chen & Chen 2010). This suggests that the external tidal field can play a significant role in shaping the GC morphology in dense environments. However, observational constraints on the shapes of our chosen GCs remain absent in the literature. Therefore, while the effect of neighbouring massive clusters and the bulge potential is physically plausible, its quantitative effect on the morphology of individual systems is still not well constrained.

In our current modelling, we restricted the simulated clusters to the list of systems that survived the early turbulent dynamical history of the Galactic centre formation and evolution. These surviving systems are clearly less affected by global common gravitational interactions than the GC systems that were already dissolved during the early Galaxy evolution on cosmological timescales. The first attempt at systematic direct N-body modelling of such already dissolved GC systems in a realistic time-evolving MW potential was presented in our previous paper, Kuvatova et al. (2026). In the follow-up studies, we plan to extend this combined evolutionary approach to the global system of dissolved and surviving Galactic centre GCs. The question of individual GC survivability is a complex multi-parameter problem. According to our current understanding, it requires not only a description of the internal dynamics of the GC, but also a correct account of interactions with other nearby systems.

Acknowledgements

The authors thank the anonymous referee for a very constructive report and suggestions that helped significantly improve the quality of the manuscript. The work was carried out with the support of the Ministry of Science and Higher Education of the Republic of Kazakhstan within the framework of Project No. AP26100669, “Galactic Archaeology of Gaia-Enceladus Globular Clusters as Indicators of the Formation History of the Milky Way Disk”. MI and PB thank the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams, grant No. PAN.BFB.S.BWZ.329.022.2023. We also gratefully acknowledge the Polish high-performance computing infrastructure PLGrid (HPC Center: ACK Cyfronet AGH) for providing computer facilities and support within computational grant No. PLG/2026/019243.

References

Appendix A Orbital evolution of Terzan 2, 4, and 5 in the TVP with dynamical friction

Here, we include additional computations aimed at quantifying the impact of dynamical friction (DF) on the determination of GCs’ orbital parameters. In an astrophysical framework, DF describes the gradual deceleration experienced by a massive object moving through a background of field stars, caused by the cumulative effect of gravitational interactions with surrounding particles. This phenomenon was first identified by Chandrasekhar and von Neumann in their seminal papers (Chandrasekhar & von Neumann 1942, 1943). Subsequently, Chandrasekhar provided a more rigorous and quantitative formulation of the theory in later works (Chandrasekhar 1943a,b). The resulting acceleration due to DF acting on a GC can be expressed as (Binney & Tremaine 2008): $$\frac{\mathrm{d}{\mathbf{V}}_{\mathrm{G}\mathrm{C}}}{\mathrm{d}\mathrm{t}}=-\frac{4\pi G^2\rho M_{\mathrm{G}\mathrm{C}}}{V_{\mathrm{G}\mathrm{C}}^3}\chi \cdot \ln \mathrm{\Lambda }\cdot {\mathbf{V}}_{\mathrm{G}\mathrm{C}}\text{\hspace{0.17em}with\hspace{0.17em}}\chi =\frac{\rho (<V_{\mathrm{G}\mathrm{C}})}{\rho }.$$(A.1)

In general, the function χ and the Coulomb logarithm Λ are functions of the velocity of the massive body, as well as of the physical characteristics of the surrounding medium. A detailed discussion of the adopted dynamical friction formalism is provided in Just et al. (2011). Motivated by the results presented in that study, we assume constant values of ln Λ = 5 and a χ = 0.5 throughout our calculations.

In Fig. A.1, we present the evolution of the orbits of our three GCs over 8 Gyr of integration time for models calculated with a dynamical friction algorithm. In the body of the main paper, this figure can be compared with Fig. 1. In Fig. A.2, we also show the evolution of the Terzan clusters’ orbital parameters a and ecc with and without DF. In our opinion, applying dynamical friction in its present form does not significantly change the Terzan cluster orbits on the global phase-space diagrams.

Fig. A.1 Orbital evolution of the Terzan clusters integrated simultaneously in the TVP external potential 411321, including DF. The trajectories are shown in the (X − Y),(X − Z), and (R − Z) planes, where R is the planar Galactocentric radius. The total integration time corresponds to a look-back time of 8 Gyr. Black markers indicate the current positions, while non-filled markers indicate the positions at −8 Gyr.

Fig. A.2 Evolution of a and ecc for the Terzan clusters integrated simultaneously in the TVP external potential 411321. Left: without DF. Right: with DF.

Appendix B Evolution of b/a and c/a axial ratios

In Fig. B.1, we show the evolution of the axial ratios, demonstrating strong correlations between the b/a and c/a ratios for the common 3xGC run and the ref. runs. The systems evolve from nearly spherical to distinctly prolate configurations (b/a and c/a ~0.75–0.80), especially in the case of Ter2. Ter5 exhibits similar but less pronounced changes. In contrast, Ter4 remains relatively stable. In the ref. run, however, the clusters retain nearly spherical shapes (axis ratios ≳ 0.95), highlighting that the triaxial deformations observed in the common 3xGC case are driven by mutual gravitational interactions among the Terzan clusters. It should be noted that the evolution of the b/a and c/a axes for Ter4 requires further modelling and research.

Fig. B.1 Time evolution of the axial ratios, b/a and c/a, for the common 3xGC run and the ref. runs.

Appendix C High-mass stellar remnant distribution for the inner Terzan cluster regions

We focused in particular on compact stellar remnants located in the central region (inner 30 pc). By employing detailed stellar evolution modules within our N-body framework, we traced the internal kinematics and long-term orbital behaviour of remnants, primarily white dwarfs (WDs), neutron stars (NSs), and black holes (BHs; Kamlah et al. 2022a). The most dynamically impactful stage in the lives of massive stars occurs during their terminal supernova explosions. Such energetic events generally impart fallback-scaled natal kicks with random orientations to the newly formed remnants (Banerjee et al. 2020).

In Fig. C.1, we show the present-day normalised cumulative number distributions of stellar remnants as a function of distance from the density centres of the GCs for the common 3xGC run and the ref. runs. In the case of the most massive GC, Ter5, we find almost no differences between these runs: the normalised cumulative distributions are nearly identical. We interpret this as a result of the very weak mutual gravitational influence of Ter2 and Ter4 on the internal structure of Ter5. In contrast, the central cumulative remnant distributions of Ter2 and Ter4 appear more compact in the common 3xGC run than in the ref. runs. This interesting behaviour of the high-mass remnant distribution can be understood in terms of the stronger gravitational influence of the massive Ter5 cluster on its lower-mass neighbours, Ter2 and Ter4. Due to periodic gravitational perturbations and stripping induced by Ter5, the internal distribution of remnants within the other Terzan clusters becomes slightly more concentrated, particularly for the WD population.

Fig. C.1 Present-day distribution of high-mass stellar remnants, WDs, NSs, and BHs, for both types of runs. Red, orange, and yellow colours represent the ref. runs. Light blue, violet, and black represent the common 3xGC run.

All Tables

All Figures

	Fig. 1 Simultaneous orbital evolution of the Terzan clusters in the TVP 411321. The orbital evolution is present in three planes, (X–Y), (X–Z), and (R–Z), where R is the planar Galactocentric radius. The total integration time is an 8 Gyr look-back time. The black marks show the current positions, and empty marks show the positions at −8 Gyr.
In the text

	Fig. 2 Evolution of Mtid and rhm for the Terzan clusters due to mass loss based on common 3xGC run (coloured lines) and ref. run (yellow lines). The dotted lines show the estimated current mass, rhm, and the grey zone shows the ±2σ observational uncertainty based on the catalogue by Baumgardt & Vasiliev (2021) for today.
In the text

	Fig. 3 Evolution of the GC semi-major axes during the simulation with marked collision events in the common 3xGC run.
In the text

	Fig. 4 Collisions of the Terzan clusters in 3D projections for different dR conditions. Left: dR = 200 pc. Right: dR = 100 pc.
In the text

	Fig. 5 Global potential energy evolution of the cluster systems for the common 3xGC run and the individual ref. runs. Left: full-time evolution of the models up to 8 Gyr. Right: same energy evolution with a very detailed output near the close passage marked A in Table 2.
In the text

	Fig. 6 Time evolution of the mutual distance between GC pairs and the axial ratios for the common 3xGC run and ref. runs. The middle and bottom panels show the triaxiality. The top and middle panels correspond to the common 3xGC run, and the bottom panel shows the ref. runs for the Terzan clusters.
In the text

Fig. A.1 Orbital evolution of the Terzan clusters integrated simultaneously in the TVP external potential 411321, including DF. The trajectories are shown in the (X − Y),(X − Z), and (R − Z) planes, where R is the planar Galactocentric radius. The total integration time corresponds to a look-back time of 8 Gyr. Black markers indicate the current positions, while non-filled markers indicate the positions at −8 Gyr.

In the text

	Fig. A.2 Evolution of a and ecc for the Terzan clusters integrated simultaneously in the TVP external potential 411321. Left: without DF. Right: with DF.
In the text

	Fig. B.1 Time evolution of the axial ratios, b/a and c/a, for the common 3xGC run and the ref. runs.
In the text

	Fig. C.1 Present-day distribution of high-mass stellar remnants, WDs, NSs, and BHs, for both types of runs. Red, orange, and yellow colours represent the ref. runs. Light blue, violet, and black represent the common 3xGC run.
In the text