# Milky Way globular clusters on cosmological timescales.
## I. Evolution of the orbital parameters in time-varying potentials

Maryna Ishchenko, Margaryta Sobolenko, Peter Berczik, Sergey Khoperskov, Chingis Omarov, Olexander Sobodar, Maxim Makukov

<figure><img src="image_1.png" /></figure>

► To cite this version:

Maryna Ishchenko, Margaryta Sobolenko, Peter Berczik, Sergey Khoperskov, Chingis Omarov, et al. Milky Way globular clusters on cosmological timescales. I. Evolution of the orbital parameters in time-varying potentials. Astronomy & Astrophysics - A&A, 2023, 673, 10.1051/0004-6361/202245117. insu-04848741

HAL Id: insu-04848741

https://insu.hal.science/insu-04848741v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Milky Way globular clusters on cosmological timescales

## I. Evolution of the orbital parameters in time-varying potentials

Maryna Ishchenko<fnref n="1" />,<fnref n="2" />, Margaryta Sobolenko<fnref n="1" />, Peter Berczik<fnref n="3" />,<fnref n="4" />,<fnref n="1" />,<fnref n="2" />, Sergey Khoperskov<fnref n="5" />,<fnref n="6" />, Chingis Omarov<fnref n="2" />, Alexander Sobodar<fnref n="1" />, and Maxim Makukov<fnref n="2" />

Received 1 October 2022 / Accepted 5 April 2023

## ABSTRACT

*Context.* Recent observational data show that the Milky Way galaxy contains about 170 globular clusters. A fraction of them is likely formed in dwarf galaxies that were accreted onto the Milky Way in the past, while the remaining clusters were formed in situ. Therefore, the different parameters, including the orbits, of the globular clusters are a valuable tool for studying the Milky Way evolution. However, because the evolution of the 3D mass distribution of the Milky Way is poorly constrained, the orbits of the clusters are usually calculated in static potentials.

*Aims.* We study the evolution of the globular clusters in several external potentials, where we aim to quantify the effects of the evolving galaxy potential on the orbits of the globular clusters.

*Methods.* For the orbit calculation, we used five Milky Way-like potentials from the IllustrisTNG-100 simulation. The orbits of 159 globular clusters were integrated using the high-order N-body parallel dynamic code ϕ-GPU, with initial conditions obtained from the recent Gaia Data Release 3 catalogues.

*Results.* We provide a classification of the globular cluster orbits according to their 3D shapes and association with different components of the Milky Way (disk, halo, and bulge). We also found that the energy – angular momentum of the globular clusters in the external potentials is roughly similarly distributed at the present time. However, neither total energy nor total angular momentum of the globular clusters are conserved due to time-varying nature of the potentials. In some extreme cases, the total energy can change up to 40% (18 objects) over the last 5 Gyr of evolution. We found that the in situ formed globular clusters are less affected by the evolution of the TNG potentials than clusters that were likely formed ex situ. Therefore, our results suggest that time-varying potentials significantly affect the orbits of the globular clusters, thus making them vital for understanding the formation of the Milky Way.

**Key words.** globular clusters: general – Galaxy: kinematics and dynamics – Galaxy: evolution – methods: numerical

## 1. Introduction

According to the ΛCDM model, the Milky Way (MW) globular clusters (GCs) are the first stellar associations that formed in the early Universe as gravitationally bound systems (Gratton et al. 2019). As a result, their typical ages are >10–12 Gyr (Marín-Franch et al. 2009; VandenBerg et al. 2013; Valcin et al. 2020), and the current masses are ≥10⁵ M☉ (Harris et al. 2013; Kharchenko et al. 2013; Baumgardt et al. 2019; Baumgardt & Vasiliev 2021). Many GCs survived over long galactic timescales, and they thus became relics from the earliest era of galaxy formation (Garro et al. 2022). Recent observations reveal a picture according to which the MW halo contains a large number of GCs (Harris 2010; Vasiliev & Baumgardt 2021), stellar streams (Ibata et al. 2021), and satellite galaxies (McConnachie 2012; McConnachie et al. 2021). Most of the discovered and studied streams in the halo are remnants of either GCs or dwarf galaxies that were destroyed in

the process of the merging history of our Galaxy (Mateu 2023; Ibata et al. 2021; Ferrone et al. 2023). They underwent complete destruction (or mass loss) under the influence of the complex time-dependent gravitational field of our Galaxy. The study of the kinematics and chemical abundances of the GCs that were captured by the MW is a very good tool for studying the history of our Galaxy (Xiang & Rix 2022), and also to search for possible progenitors for these GCs (Massari et al. 2019; Malhan et al. 2022).

The publication of the Gaia (ESA) catalogues (Gaia Collaboration 2021, 2018) made full phase-space information available for almost the entire population of the GCs, which enables an investigation of their orbits. Several recent works have studied the orbits of the GCs in the MW potential under the assumption that it axisymmetric and does not evolve in time (Massari et al. 2019; Myeong et al. 2019; Malhan et al. 2022). These two assumptions result in the conservation of the total energy and of the z-component of the

<footnote n="1">¹ Main Astronomical Observatory, National Academy of Sciences of Ukraine, 27 Akademika Zabolotnoho St, 03143 Kyiv, Ukraine e-mail: marina@mao.kiev.ua</footnote>

<footnote n="2">² Fesenkov Astrophysical Institute, 050020 Almaty, Kazakhstan</footnote>

<footnote n="3">³ Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, Mönchhofstrasse 12-14, 69120 Heidelberg, Germany</footnote>

<footnote n="4">⁴ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), MTA Centre of Excellence, Konkoly Thege Miklós út 15-17, 1121 Budapest, Hungary</footnote>

<footnote n="5">⁵ Leibniz Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany</footnote>

<footnote n="6">⁶ GEPI, Observatoire de Paris, Université PSL, CNRS, 5 place Jules Janssen, 92190 Meudon, France</footnote>
angular momentum ($L_z$) of the GCs. This has been used to group some of the GCs and link these kinematically coherent groups with their galaxy progenitors (either MW or accreted satellite). However, these assumptions are not fully correct because the MW galaxy contains the bar and boxy-peanut bulge in its centre (Dwek et al. 1995; Wegg & Gerhard 2013; Ness & Lang 2016). In the outer parts, the axisymmetry of the MW potential is thought to be broken by massive orbiting satellites (Laporte et al. 2018, 2020), including the large and small Magellanic Cloud systems (Gómez et al. 2015; Petersen & Peñarrubia 2021; Conroy et al. 2021). In the current picture, most of the mass of the disk is assembled from the smooth accretion of gas combined with the accretion along cold filaments (Katz & Gunn 1991; Birnboim & Dekel 2003; Kereš et al. 2005; Agertz et al. 2009). Although the mass evolution of the MW is still uncertain, different models suggest that it rapidly increases at early times (until ≈8 Gyr), followed by a slower evolution (Snaith et al. 2014), which agrees with the mass growth of MW-like galaxies that is constrained via halo abundance matching (van Dokkum et al. 2013; Fattahi et al. 2019). While some models suggest that the mass growth of MW-type galaxies affects the coherence of accreted and in situ stellar populations in the kinematic space (e.g. in $E - L_z$ Grand et al. 2019; Panithanpaisal et al. 2021; Khoperskov et al. 2023; Pagnini et al. 2023), the impact on the motions of the GCs remains unclear.

The main feature of our study of the MW GC subsystem phase-space distributions is the time-variable external potential. This feature significantly differs in our investigation from other similar works (see e.g. Massari et al. 2019; Vasiliev 2019; Bajkova & Bobylev 2021; Bajkova et al. 2021; Malhan et al. 2022; Armstrong et al. 2021; Haghi et al. 2015), in which the authors used a fixed or analytically time-evolving MW potential. In a series of articles about the evolution of GCs, and in particular in this article, we aim to investigate how the evolution of the mass and the spatial scales of both stellar and dark matter components of the MW affect the orbital parameters of the GCs. Because the evolution of the MW parameters is poorly constrained, we analyse the dynamics of the GCs in the MW-like evolving potentials as directly obtained from IllustrisTNG cosmological simulations (Pillepich et al. 2018; Nelson et al. 2018).

The paper is organised as follows. In Sect. 2.1 we present the GCs selection with the error distribution from the catalogue of the Gaia data release 3 (DR3). In Sects. 2.2 and 2.3 we introduce the time-dependent gravitational potentials of the MW-type galaxies that were selected from IllustrisTNG-100 and briefly describe the numerical integrator. In Sect. 3 we present the individual orbits of the GCs in the evolving TNG potentials. In Sect. 4 we describe the GC phase-space evolution. In Sect. 5 we summarise our results.

## 2. Methods

### 2.1. Globular cluster sample

We used the recent catalogues of GCs released by Vasiliev & Baumgardt (2021) and Baumgardt & Vasiliev (2021), which contain information about more than 160 individual objects. The catalogues include the complete 6D phase-space information required for the initial conditions of our simulations: right ascension (RA), declination (Dec), heliocentric distance ($D_\odot$), proper motion in right ascension (PMRA), proper motion in declination (PMDEC), and the radial velocity (RV). Because of the dynamical simulation of the GCs, we also used the self-gravity of the

clusters together with the Galaxy external potential. We excluded the Mercer 5 object because there is no mass information for this GC in the catalogues listed above.

Before the orbital integration, we analysed the GCs Gaia measurement error influence. This error might change the orbital shape during the orbital integration (see Bajkova & Bobylev 2021). We analysed the errors in PMRA, PMDEC, RV, and $D_\odot$ for 159 GCs from the catalogue above. In Fig. 1 we present the distribution of the relative errors of these parameters. The vast majority (up to 95%) of the errors are well below 5%, but some GCs (Crater and Lae 3) have quite a high value in proper motion. Analysing the errors in the heliocentric distance, we found only 18 objects with relative errors larger than 5% and only three GCs with errors larger than 10% (namely 2MASS-GC01, VVV-CL001, and Pal 10). The influence of the measurement errors on the orbital integration is presented in Sect. 3.1.

To transform the positions and velocities into the Cartesian galactocentric rest frame (see Johnson & Soderblom 1987; Bovy 2011, for the coordinate transformation equations), we assumed a galactocentric distance of the Sun $R_\odot = 8.178$ kpc (Gravity Collaboration 2019; Reid & Brunthaler 2004), a height above the Galactic plane $Z_\odot = 20.8$ pc (Bennett & Bovy 2019), and the velocity of the local standard of rest (LSR) $V_{LSR} = 234.737$ km s$^{-1}$ (Bovy et al. 2012; Drimmel & Poggio 2018). Accordingly, the Sun is located at $X_\odot = -8178$ pc, $Y_\odot = 0$ pc and $Z_\odot = 20.8$ pc in our Cartesian galactocentric coordinate system. For the peculiar velocity of the Sun, we used the following values with respect to the LSR: $U_\odot = 11.1$ km s$^{-1}$, $V_\odot = 12.24$ km s$^{-1}$, and $W_\odot = 7.25$ km s$^{-1}$ (Schönrich et al. 2010).

### 2.2. Time-variable gravitational potentials from the Illustris simulation

Because the time evolution of the MW potential is poorly constrained from observations, we selected the MW-like galaxies from the most recent Illustris cosmological numerical simulations to integrate the GC orbits. For our purpose, we used the publicly available snapshot data from the IllustrisTNG-100 (Pillepich et al. 2018; Nelson et al. 2018, 2019; Springel et al. 2018; Marinacci et al. 2018; Naiman et al. 2018). These simulations successfully reproduce the main scaling relations and evolution of galaxy sizes (Genel et al. 2018), including the formation of realistic disk galaxies (Pillepich et al. 2019), the gas-phase mass-metallicity relation (Torrey et al. 2019), and the range of the observed kinematic properties of MW-type galaxies (Lovell et al. 2018). The Illustris simulations were used to investigate various processes of galactic evolution, including the gas-stripping phenomena (Yun et al. 2019; Łokas 2020), the central black hole activity (Habouzit et al. 2019; Donnari et al. 2021), and star formation quenching (Genel et al. 2018; Davies et al. 2020). These results ensure that the investigation of the properties of different types of galaxies from the Illustris cosmological simulations allows us to study the GC dynamics in an environment similar to the MW.

The IllustrisTNG-100 is characterised by a simulation box of approximately 100 Mpc$^3$. In a box of this size, a sufficient number of the MW-mass disk galaxies can be resolved (in our sample, 54 objects$^1$) with a mass resolution of $7.5 \times 10^6 M_\odot$ for dark matter and $1.4 \times 10^6 M_\odot$ for the baryonic particles. For our analysis, we identified the MW-like galaxy candidates in the Illustris simulations with at least $10^5$ dark matter particles and at least

<footnote>$^1$ MW-like potentials are presented on the web page of the project https://bit.ly/3b0lafw</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Distribution of the relative errors (from left to right): radial velocity (eRV), proper motions in right ascension (ePMRA), and declination (ePMDEC), and in heliocentric distance $eD_⊙$. Blue histograms represent the error distributions. The following GCs are not shown because the error values are too high: Crater, Lae 3, NGC 6760, BH 261, NGC 6553, and Pal 13, except in the graph in the bottom right panel.</figcaption></figure>

<table><caption>**Table 1.** Parameters of the time-varying potentials selected from the IllustrisTNG-100 simulation at redshift zero.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>#411321</th><th>#441327</th><th>#451323</th><th>#462077</th><th>#474170</th><th>Milky way</th></tr></thead><tbody><tr><td>Disk mass, M<sub>d</sub></td><td>10<sup>10</sup> M<sub>☉</sub></td><td>7.110</td><td>7.970</td><td>7.670</td><td>7.758</td><td>5.825</td><td>6.788</td></tr><tr><td>Halo mass, M<sub>h</sub></td><td>10<sup>12</sup> M<sub>☉</sub></td><td>1.190</td><td>1.020</td><td>1.024</td><td>1.028</td><td>0.898</td><td>1.000</td></tr><tr><td>Disk scale length, a<sub>d</sub></td><td>1 kpc</td><td>2.073</td><td>2.630</td><td>2.630</td><td>1.859</td><td>1.738</td><td>3.410</td></tr><tr><td>Disk scale height, b<sub>d</sub></td><td>1 kpc</td><td>1.126</td><td>1.356</td><td>1.258</td><td>1.359</td><td>1.359</td><td>0.320</td></tr><tr><td>Halo scale height, b<sub>h</sub></td><td>10 kpc</td><td>2.848</td><td>1.981</td><td>2.035</td><td>2.356</td><td>1.858</td><td>2.770</td></tr></tbody></table>

**Notes.** The last column shows the parameters of the corresponding MW components according to Bennett et al. (2022) at present.

10³ baryonic particles (stars and gas) at redshift zero. From the simulated galaxies in the IllustrisTNG-100, we selected our five galaxies (411321, 441327, 451323, 462077, and 474170) that reproduce MW parameters well at redshift zero (disk and halo masses with the spatial scales) best at present (see Table 1).

Due to the limited particle number of our selected TNG-TVPs, we cannot resolve a separate bulge component. This is mainly caused by the limited mass resolution of the IllustrisTNG-100, that is, we simply do not have enough particles to determine this component as a separate bulge in our galaxy mass model.

We also used a circular velocity value at the solar distance (≈8 kpc) in the model as an additional parameter to select the TNG galaxies that represent the MW-type systems best. This value indicates the position of the Sun at present. According to the age and chemical compositions of the stars in the solar neighbourhood, we know that over the past few billion years, there were no large changes in the radial motion of the masses. This means that the circular velocity at the distance of the Sun in the Galactic disk should remain approximately constant during the last few billion years near the $V_☉$ ≈ 235 km s⁻¹ (Mardini et al. 2020). In Appendix A, we present the circular velocity evolution for the selected five TNG external potentials. In the same

figures, we also present the time evolution of the selected TNG-TVP circular velocity radial distribution. These figures were obtained directly from the mass distribution of the corresponding IllustrisTNG-100 simulation snapshots for three different redshifts $z = 0, -5$ Gyr ($z = 0.48$), and $-10$ Gyr ($z = 1.74$). Because the galaxy mass decreases as a function of lookback time, the circular velocity decrease at all galactocentric radii is clearly visible.

To obtain the spatial scales of the disks and dark matter haloes, we decomposed the mass distribution using the Miyamoto-Nagai (MN) $Φ_d(R, z)$ (Miyamoto & Nagai 1975) and Navarro-Frenk-White (NFW) $Φ_h(R, z)$ (Navarro et al. 1997) potentials,

$$ \begin{aligned} \Phi_{\text{tot}} &= \Phi_d(R, z) + \Phi_h(R, z) \\ &= -\frac{GM_d}{\sqrt{R^2 + (a_d + \sqrt{z^2 + b_d^2})^2}} - \frac{GM_h \ln(1 + \frac{\sqrt{R^2+z^2}}{b_h})}{\sqrt{R^2+z^2}}, \end{aligned} \quad (1) $$

where $R = \sqrt{x^2 + y^2}$ is the planar galactocentric radius, $z$ is the distance above the plane of the disk, $G$ is the gravitational constant, $a_d$ is the disk scale length, $b_{d,h}$ are the disk and halo scale
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Evolution of halo and disk masses, and their characteristic scales for (from top to bottom): 411321, 441327, 451323, 462077, and 474170 TNG-TVPs in time. The halo mass $M_h$, disk mass $M_d$, NFW halo scale parameter $b_h$, and the MN disk scale parameters $a_d$ and $b_d$ are presented from left to right. Solid green lines with dots show the parameters recovered from the IllustrisTNG-100 data. Dotted and dash-dotted blue lines correspond to the values after the interpolation and smoothing with the 1 Myr time step that was used in the orbital integration.</figcaption></figure>
heights, respectively, and $M_d$ and $M_h = 4\pi\rho_0 b_h^3$ ($\rho_0$ is the central mass density of the halo) are the masses of the disk and halo, respectively.

In our case, all these mass distribution parameters are func-
tions of lookback time. We performed this fit for each available
snapshot in each of our selected galaxies. The code that we used
to find these parameters is also publicly available at GitHub<fnref n="2" />.

For each snapshot available in the IllustrisTNG database, we
decomposed the mass distribution of each selected galaxy. Next,
we interpolated the variation in masses and spatial scales of dark
matter and stellar disk among the snapshots with a time step of
1 Myr. Finally, to avoid some short-timescale perturbations that
are likely linked to the mergers or errors in the galaxy identifica-
tion, we additionally smoothed the evolution of the main galaxy
parameters. These procedure results are shown for our five time-
variable potentials (TNG-TVPs) in Fig. 2.

## 2.3. Integration procedure

For the GC orbit integration, we used the high-order parallel dynamical *N*-body code *φ*-GPU<fnref n="3" /> (Berczik et al. 2011, 2013). The code is based on the fourth-order Hermite integration scheme with hierarchical individual block time steps. To calculate the force that acts on the particles during the integration, we combined the particle self-gravity with the time-variable external potential described in the previous subsection. To calculate the self-gravity (as a first-order simplification), we applied the current mass of the individual GCs. We applied the integration time-step parameter $\eta = 0.01$ (Makino & Aarseth 1992), which gives the required integration accuracy for our investigation. As a test case, we ran a typical simulation with the fixed (i.e. for *t* = 0) values of the current external potential 411321 with the most recent halo and disk masses and scale lengths. During this test case, the total relative energy drift ($\Delta E/E_{t=0}$) over a 10 Gyr backward integration was below $\sim 2.5 \times 10^{-13}$. For the production runs, we embedded our 159 GC into the selected five external TNG-TVPs as a point masses (set the masses of GC as for today). For each GC point, we set the central values from Baumgardt & Vasiliev (2021) and Vasiliev & Baumgardt (2021) for the position and velocity.

# 3. Characteristics of the individual globular cluster orbital evolution in five time-varying potentials

## 3.1. Influence of the measurement errors

Prior to the analysis of the orbital types of the GCs, we inves-
tigated the impact of the initial condition uncertainties on the
initial velocity values and structure of the orbits. For each GC
in our sample, we generated an additional ten initial conditions
in which the initial velocities take the normal distribution of the
errors for PMRA, PMDEC, RV (blue circles), and $D_\odot$ (green
circles) from Baumgardt & Vasiliev (2021)<fnref n="4" /> into account; see
Fig. 3. The two randomisations have a similar effect on the
GC initial velocities. In general, the errors for the velocities
are even slightly dominated by the errors due to the distance
determination. The effect from the velocity determination errors

is even slightly stronger for a few objects (six) with indexes
below 30.

In Fig. 4 we illustrate the impact of the initial conditions
uncertainties on the orbital evolution of two selected GCs:
Crater (belongs to the 5% with the largest velocity uncertainties:
211%, 100% and 0.4% for the proper motions and radial veloc-
ity) as for the NGC 6121 (with one of the lowest values: 0.1%,
0.12%, and 2%). These two GCs are marked by crosses in Fig. 3.
The velocity uncertainties significantly change the orbital char-
acteristics in the case of Crater. In particular, the extension of the
orbit along the principal axes, especially in the disk plane, is very
different in various GCs realisations. At the same time, the verti-
cal motions are less affected. In contrast, for NGC 6121, whose
parameters are well measured, the orbits are less affected by the
velocity uncertainties, as shown in the bottom row of Fig. 4.
Thus, since the velocity uncertainties for most of our GCs are
similar to that of NGC 6121, we assume that the orbits of the
majority of the GCs are robust to the variations of the initial
condition.

Figure 4 shows that the influence of the relative errors of the heliocentric distances on the orbital shape of our GC sample is indeed less significant than the other relative errors. As an example, we present the orbits of two selected GC: Carter (largest heliocentric distance) and NGC 6121 (smallest heliocentric distance). The effect from the relative distance error is significantly smaller than the effect from the proper motions and radial velocity errors that is presented in the same plot. Next, we focus our attention only on the GC-point orbits obtained using the central values for the positions and velocities from Baumgardt & Vasiliev (2021) and Vasiliev & Baumgardt (2021).

## 3.2. Evolution of the globular cluster orbital elements

Assuming the cosmologically motivated time variable poten-
tial, we automatically took the possible influence of the Milky
Way satellite galaxies on the GC subsystem into account. If we
compare the individual GC orbital elements (semi-major axis a
and eccentricity *ecc*) variation over the whole integration time
(sometimes more than hundreds or thousands of GC orbital reve-
lations). In Fig. 5 we present the orbital element relative changes
in five TNG-TVPs. In Fig. 6 we show the time evolution of a and
ecc. The strong influence on the GC orbits has been discussed
in detail in Garrow et al. (2020) and Boldrini & Bovy (2022).
Our Figs. 5 and 6 are directly comparable with Figs. 1–3 from
Garrow et al. (2020). The relative changes in a and ecc lie in
a similar range. In our case, the a evolution is even larger by a
factor of two. We also note the very similar distribution of the
relative quantities for the different TNG-TVP potentials.

As an illustration of the time evolution of the orbital elements of all GCs, we show the data for 411321 TNG-TVP in Fig. 6. The separation between the inner and outer GCs is very clear. The inner GCs ($a \lesssim 3$ kpc) have more regular and larger eccentricity changes during the evolution. The outer GCs ($a \gtrsim 3$ kpc) have much smaller eccentricity changes during the whole backward integration time. This Fig. 6 can be directly compared with Fig. 3 from Garrow et al. (2020). Our potential keeps the individual GCs bound, even with large a during the whole integration time of 10 Gyr.

For example, Garrow et al. (2020) and Boldrini & Bovy (2022) added the main MW satellites during their GC orbit investigations as an additional gravitational perturbation in the fixed MW potential approximation. In contrast in our case, the IllustrisTNG-100 mass-assembling history of our theoretical

<footnote n="2">² The ORIENT: https://github.com/Mohammad-Mardini/The-ORIENT</footnote>

<footnote n="3">³ *N*-body code φ-GPU: https://github.com/berczik/phi-GPU-mole</footnote>

<footnote n="4">⁴ Error values for PMRA, PMDEC, RV, and $D_\odot$ from https://people.smp.uq.edu.au/HolgerBaumgardt/globular/orbits_table.txt</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Influence of the error randomisation on the initial velocity components in the galactocentric Cartesian reference frame. The distance determination is shown by green circles and the velocity determination by blue circles. The red cross marks Crater, and the black cross indicates NGC 6121.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Changes in the orbits caused by the uncertainties in the radial velocity and proper motions (see footnote 4) for Crater and NGC 6121. The presented orbits were integrated for a 10 Gyr lookback time. The different colours with transparency represent the orbits for five random realisations of the initial conditions in 411321 TNG-TVP.</figcaption></figure>

MW-like galaxies (see Fig. 2) provide us with the same or even a larger scale of orbital perturbations for individual GCs during their time evolution.

In Fig. 7 we present the evolution of the main orbital GC parameters: apocenters and pericenters (for a similar plot see Bajkova et al. 2021, see Fig. 9). On the X-axis, we plot the apo and peri values of each GC in the static (fixed) potential, that is, we use the initial values of the 411321 IllustrisTNG-100 potential. The general behaviour of these quantities in our Fig. 7 and in

the work of Bajkova et al. (2021) is quite similar. Up to roughly 5 Gyr lookback time, no significant changes in these quantities are visible. After approximately 8 Gyr, however, the differences (especially in apocenter values) become more significant (up to ~25%).

In Fig. 8 we present the evolution of the apocenters and pericenters during the full 10 Gyr lookback-time integration in 411321 TNG-TVP. For the visualisation, we selected several GCs with different heliocentric distances. In our representative
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Relative orbital semi-major *Δa/a* (left) and eccentricity *Δe/e* (right) changes during the GC orbital evolution for all five TNG-TVPs.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Evolution of the GC orbital semi-major and eccentricity during the whole backward integration time in the case of 411321 TNG-TVP. The inner GCs ($a \le 3$ kpc) have more regular and larger eccentricity changes during the evolution. The outer GCs ($a > 3$ kpc) have much smaller eccentricity changes during the whole backward integration time.</figcaption></figure>

sample, Liller and Terzan 4 are the closest GCs to the Galactic centre.

Figure 8 shows that Liller and Terzan 4 are the closest GCs with a clearly decreasing apocenter with minimum values at 8 Gyr with further increase. NGC 6717 slowly declines with the minimum at the same time. However, the pericenter evolution is less clear. This behaviour can be understood by taking the stronger dependence of the apocenter value (compared to the pericenter) on the mass and size changes of our TNG-TVPs into account. At the same time, the GC 5927 apo- and pericenter changes are synchronised and have a minimum value at 3 Gyr and after it starts to increase. We present the evolution of the

apocenters and pericenters for all GCs on the website<fnref n="5" /> for all five TNG-TVPs and for each GC separately.

We also present in Fig. 9 the comparison of apocenters and pericenters of the evolution for selected GCs in the five TNG-TVPs simultaneously. The pericenter evolution in all the five TNG-TVPs looks very similar. The apocenter evolution shows some differences for different potentials, which we can understand as a stronger influence of the time-variable masses and sizes of MW-like potentials on the orbits.

### 3.3. Types of orbits

We calculated the orbits of 159 GCs at a lookback time for 10 Gyr in each of our TNG-TVP potentials. The orbital evolution for our GC sample is presented online<fnref n="5" />. The visualisation was carried out in three projections of the Cartesian galactocentric rest frame: (X, Y), (X, Z), and (R, Z), where R is the planar galactocentric radius. The orbits change in different potentials, but their general shape remains similar.

The shape filled by an orbit in an axisymmetric potential might be classified as a short-axis tube, thin tube, and so on, depending on specifics of the potential (e.g. oblate or prolate) and the particular combination of the integrals of motion the orbit possesses (Carpintero & Aguilar 1998; Merritt 1999). For our purposes, after a visual analysis of 159 orbits (in each of our TNG-TVP potentials), we divided them into four main types (categories):

- Tube orbit (TB): One hundred and ten (69%) orbits. The orbit is mainly in the (X, Y) Galactic plane with a hole in the centre. The orbit in the (X, Z) plane has a boxy shape and a trapezoidal shape in the (R, Z) plane.

- Perpendicular tube orbit (PT): Eight (5%) orbits. The orbital plane is close to the meridional, also contains a hole in the centre, and has a trapezoidal shape in the (R, Z) plane.

- Long radial orbit (LR): Nineteen (12%) orbits. These are long-period, near-hyperbolic orbits that are characterised by a high eccentricity without a prominent hole in the centre.

- Irregular orbit (IR): Twenty-two (14%) orbits. The orbit cannot be classified as one of those described above.

<footnote n="5">⁵ We also present the orbits for all 159 GCs in TNG-TVP potentials on the web site of the project https://bit.ly/3b0lafw.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Evolution of the apocenters (left panel) and pericenters (right panel) for the 159 GCs in the 411321 TNG-TVP. The colour-code corresponds to the lookback time.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Evolution of the apocenters (upper panel) and pericenters (bottom panel) for the (from left to right) Liller, Terzan 4, NGC 6717, and NGC 5927 GCs. The presented orbits were integrated for 10 Gyr lookback time. Each point represents the individual orbital passages of the GC in the 411321 TNG-TVP.</figcaption></figure>

We list the orbit types for each GC in Table C.1. In Fig. B.1 we show some typical examples of the orbits we used in the classification above.

### 3.4. Association with the regions of the Galaxy

We assumed that the Galaxy has several spatial components: the bulge, the thin disk, the thick disk, and the main halo. Each GC orbit can be associated predominantly with one of them (based on the distance criteria). This associations weakly

depend on orbits shape. We analysed the orbits of all 159 GCs in each potential to sort them by ownership to the different regions of the Galaxy in the past billion years. For this classification, we used the following characteristic scales according to Bland-Hawthorn & Gerhard (2016) for the orbital apocenter and pericenter of each GC:

- Bulge (BL). The scale height in Z directions is ~0.2 kpc. The scale length in the *R* plane is ~0.7 kpc. With these scales, we determined the boundary dimensions for the bulge as ~2.1 × ~0.6 kpc. We associated nine GCs.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Evolution of the apocenters (upper panel) and pericenters (bottom panel) for (from left to right) Liller, Terzan 4, NGC 6717, and NGC 5927 GCs in five TNG-TVPs, where red shows 411321, blue shows 441327, orange shows 451323, green shows 462077, and black shows 474170, respectively. The presented orbits were integrated for 10 Gyr lookback time.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Distribution of the GCs by their orbit types in four Galactic components: bulge (BL), thin disk (TN), thick disk (TH), and halo (HL). The colour-coding shows different orbit types: tube orbit (TB), perpendicular tube orbit (PT), long radial orbit (LR), and irregular orbit (IR).</figcaption></figure>

– Thin disk (TN). The scale height in Z directions is ~0.3 kpc. The scale length in the *R* plane is ~2.6 kpc. With these scales, we determined the boundary dimensions for the thin disk as ~7.8 × ~0.9 kpc. We associated ten GCs.

– Thick disk (TH). The scale height in Z directions is ~0.9 kpc. The scale length in the *R* plane is ~2 kpc. With these scales, we determined the boundary dimensions for the thick disk as ~6 × ~2.7 kpc. We associated 47 GCs.

– Halo (HL). More than ~2.7 kpc in Z direction and more than ~7.8 kpc in the *R* plane. We associated 94 GCs.

Figure 10 shows that most GCs at the present stage (last billion year) belong to the halo and thick disk, with 59% and 30%, respectively. The bulge and thin disk contain just 5% and 6% GCs. In the distant past, however, the probability is high that some of the GCs belonged to other regions of the Galaxy (e.g. Bajkova et al. 2020; Sun et al. 2023). A possible explanation is

that GCs survive better when they are far from the Galactic centre. We list the associations with different Galaxy components for each GC in Table C.1. This classification is also stable for each of the five different TNG-TVP, which makes these potentials quite representative.

The majority of the GCs from the bulge and the thin and thick disk have tube orbits (see Table C.1). They might have undergone extreme encounter(s) with another GC, which changed their orbit. In contrast, halo GCs display all types of orbits. The most common is a tube orbit, together with irregular and long radial orbits. The halo also harbours objects with all GCs that have perpendicular tube orbits.

### 3.5. Energy changes in globular clusters during evolution

Because we considered the dynamics of the MW GCs in the evolving potentials, the integrals of motions (energy and angular momentum) are not conserved over time. In order to quantify this, we calculated the relative energy changes $\Delta E/E$ during the evolution of GCs over 10 Gyr. In Fig. 11 we present the relative energy changes in individual GCs during the orbit integration. In general, the effect of the Galaxy mass changes in our selected TNG-TVPs (Fig. 2). The maximum relative energy changes are ~45–50%. As an example, in Fig. 11 we highlight five GCs in colour: BH 140 and VVV CL001 (red and blue), which are associated with the Galactic disk and likely formed in situ; AM 4, Sagittarius II, and Lae 3, which are associated with the halo and likely were born outside the Galaxy (see the definition of association in Malhan et al. 2022; Fernández-Trincado et al. 2021). All other 154 GCs are plotted in grey. Figure 11 shows that in situ GCs are less strongly affected by the changes in the gravitational potential than the accreted clusters. This trend is especially visible in the last panel of the Fig. 11 for 474170 TNG-TVP.
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Relative energy $\Delta E/E$ changes during the GC orbital evolution in percent. From left to right from top to bottom: 411321, 441327, 451323, 462077, and 474170 TNG-TVPs. The colour-coding is red for BH 140, green for AM 4, blue for VVV CL001, black for Sagittarius II, orange for Lae 3, and grey for all other GCs.</figcaption></figure>

We also assembled a list of GCs whose relative energy changed by more than 40%. For simplicity, we present the list of GCs from 411321 TNG-TVP only.

- $\Delta E/E = 45\%...50\%$: Crater, Pal 3, Sagittarius II (black in Fig. 11) Eridanus, AM 1, and Pal 4;

- $\Delta E/E = 42\%...44\%$: NGC 2419, Pyxis, Pal 14, Lae 3 (orange in Fig. 11), and Whiting 1;

- $\Delta E/E = 40\%...41\%$: NGC 5694, Arp 2, Terzan 8, Terzan 7, NGC 7006, Pal 12, and Pal 13.

All of them belong to the halo (HL), and in most cases, they have either irregular (IR) or perpendicular tube-type (PT) orbits (see Table C.1).

# 4. Globular clusters in phase space

## 4.1. Present-day properties of the globular clusters

We analysed the GC phase-space distributions in the energy- angular momentum coordinates for each of the TNG-TVPs. Below, we consider specific values of energy and angular momentum normalised by the GC mass. In particular, we anal- yse the phase-space distribution of the GCs at present in three phase-space coordinate combinations: total angular momen- tum ($L_{tot}$) versus total energy ($E$), z-component of the angular momentum ($L_z$) versus total energy ($E$), and the perpendicular

component of the angular momentum ($L_{perp}$) versus total energy ($E$), where $L_{tot} = \sqrt{L_x^2 + L_y^2 + L_z^2}$ and $L_{perp} = \sqrt{L_x^2 + L_y^2}$. The orientation of the Z-axis in our Cartesian galactocentric coordi- nates is directed towards the Galactic south pole. By this defini- tion, the angular momentum of the Sun is negative.

Next, we discuss the main relations in the energy-angular momentum space for the GCs of different orbital types. In Fig. 12, the GCs are colour-coded according to the type of orbit in the 411321 potential. GCs with TB orbits have better GCs bounds. They orbit close to the Galactic centre because their total energies are highest (by absolute value). Most of the GCs with TB orbits show from zero to some prograde rotation with the smallest perpendicular angular momentum component. This might indicate their common dynamical origin or possi- ble orbital evolution as a consequence of their similar angular momentum loss timescales.

Globular clusters with IR orbits show a wide range of total energy and angular momentum components and overlap with GCs with other types of orbits. Their origin is probably not the same because their angular momentum covers quite a wide range.

The GCs with LR orbits have a wide range of total energies, but tend to have lower absolute values of the angular momentum. Interestingly, these GCs have a relatively narrow range of $L_{perp}$ with a weak trend for a retrograde rotation.
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Distribution of orbital types at redshift zero for 411321 TNG-TVP. From left to right: total angular momentum vs. energy ($L_{tot}, E$), zth component of the angular momentum vs. total energy ($L_z, E$), and perpendicular component of the angular momentum vs. total energy ($L_{perp}, E$). The colour corresponds to the type of the GC orbit: blue for the tube orbit (TB), orange for the perpendicular tube orbit (PT), green for the long radial orbit (LR), and red for the irregular orbit (IR). Pal 3, Crater, and Sagittarius II are not shown in ($L_{tot}, E$) and ($L_{perp}, E$) phase space, and Sagittarius II is not shown in the ($L_z, E$) phase space either because their values are extremely high.</figcaption></figure>

The smallest group of the GCs with PT are grouped together in all the energy-angular momentum phase coordinates presented in Fig. 12. They have roughly the same total energy, substantial net rotation, and strong motions perpendicular to the disk plane. The high values of $L_{perp}$ might suggest a common origin for them. As described above, the GCs have the same type of orbits in different TNG-TVPs. Therefore, the described picture is quite similar in all external potentials considered in the paper.

Next, we show the association of the GCs with different Galaxy components presented in the energy-angular momentum phase coordinates. In Fig. 13 (top row), we colour-code the GCs associated with bulge (BL, yellow), thin disk (TN, blue), thick disk (TH, green), and halo (HL, red) Galactic regions. For simplicity, we present the GC parameters only in the 411321 TNG-TVP. The top panel in Fig. 13 shows that as expected, the GCs that belong to the bulge (BL) and thick disk (TH) have the highest negative energy values in the entire sample. The bulge GCs have a relatively narrow energy range than the thick-disk GCs, however. The GCs associated with the halo form a system that is less strongly bound with the Galaxy. Their energy and angular momentum ranges are wide.

In Fig. 13 we present a comparison of the GC positions in phase space for the other four 441327, 451323, 462077, and 474170 TNG-TVPs in the bottom panel. The colours (as in the top panel) represent the associations with the different regions of the Galaxy. Different symbols correspond to different potentials. Comparing the GCs in different potentials at the present, we can conclude that there are no significant differences in their values in the angular momentum space. The visible differences in energy values can easily be explained with different parameters of the galactic potentials (see Table 1) at the present. At the same time, our sample of external potentials shows no strong mutual discrepancies.

## 4.2. Globular cluster subsystem phase-space evolution in five time-varying potentials

Using the time-evolving potentials, we can explore the effect of the Galactic parameter time-variation on the GC phase-space distributions in energy and angular momentum. In Figs. 14 and 15 we show the GC system snapshots for different times

(0.0, 2.5, 5.0, 7.5, and 10 Gyr lookback time) for each of our TNG-TVPs.

Figures 14 and 15 show that all our external potentials provide roughly similar phase-space energy values for individual GCs. The difference between potentials becomes more defined during the evolution. In general, we can expect that at the beginning of the backwards integration ($t = 0$ Gyr), the energies of the GCs are more negative (i.e. their bounding energies are higher; magenta) and at the end of the simulation (10 Gyr, red) their energies are more positive. This would mean that GCs are less bounded, which is not the case. Our expectations are fully met by the potential 441327 (Fig. 14, second row). The other TNG-TVP external potentials show a slightly different behaviour. For example, the 451323 (Fig. 14, third row) and 462077 (Fig. 15, first row) external potentials have the highest bounding energies at 5 Gyr (green). We can explain this behaviour by verifying the evolution of the mass and characteristic scales from Fig. 2. For these two MW-like potentials around 5 Gyr, the mass is higher and the scale sizes are smaller (as for the halo and disk).

For the potential 411321 (Fig. 14, first row), the energy evolution shows a more complex behaviour. The energy values for GCs first decrease (approximately 3 Gyr; orange and approximately 5 Gyr; green) and start to rise only after ~7 Gyr (cyan). This complex behaviour might also be understood based on the evolution of the halo and disk masses shown in Fig. 2 for this TNG-TVP.

In the case of the 474170 TNG-TVP (Fig. 15, second row), the halo mass and characteristic scale are almost constant in the beginning, in contrast to the previous external potentials (up to approximately 3 Gyr; magenta-orange).

In general, we conclude that the evolution of the GC phase-space distributions weakly varies in the five selected TNG-TVPs. The differences are clearly caused by the slightly different accretion history, which originally comes from the IllustrisTNG-100 simulation.

## 5. Conclusions

In this paper, we studied the orbital evolution of 159 Milky Way globular clusters up to 10 Gyr lookback time. For the initial positions, velocities, and masses of each GC, we used the values from the catalogues (Vasiliev & Baumgardt 2021), which
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Position of GCs in different regions of the Galaxy at the present time for the 441321 TNG-TVP, top panel. From left to right: total angular momentum, zth component of the angular momentum, and perpendicular component of the angular momentum vs energy. The colour-coding represents the GCs that belong to the different regions in the Galaxy: orange for the bulge (BL), blue for the thin disk (TN), green for the thick disk (TH), and red for the halo (HL). Bottom panel: same values, but for the 441327 (plus), 451323 (cross), 462077 (star), and 474170 (empty square) TNG-TVPs. The GCs Pal 3, Crater, and Sagittarius II are not shown in ($L_{tot}, E$) and ($L_{perp}, E$) phase spaces and Sagittarius II is not shown in ($L_z, E$) phase space because their values are extremely high.</figcaption></figure>

contain full 6D information from Gaia DR3. Analysing the errors and their distribution for the proper motions of right ascension and declination, heliocentric distances, and radial velocities, we found that the vast majority (up to 95%) of the GCs have relative errors below 5%, which allowed us to study the orbits with great precision.

For the orbital integration, we used the high-order parallel dynamical N-body code ϕ-GPU. To mimic the evolution of the MW mass distribution over time, we used the data from the IllustrisTNG-100 cosmological simulation. From the simulated subhaloes, we selected galaxies that best reproduce the current Milky Way parameters. To obtain the spatial scales of the disks and dark matter haloes, we decomposed the mass distribution using the MN and NFW potentials. Moreover, to verify that the chosen TNG potentials are similar to the potential of our Galaxy, we added an additional parameter: the observed circular velocity of the rotation at the solar radius. Finally, we selected the five TNG galaxies 441321, 441327, 451323, 462077, and 474170 as time-variable external potentials for our study.

The main results of this study are listed below.

- We integrated the orbits of 159 MW GCs up to 10 Gyr evolution in five external potentials. According to our results, we selected a set of GCs that currently belong to the four main structural components of the Galaxy. We defined 9 GCs as currently belonging to the bulge (BL), 10 as belonging to the

thin disk (TN), 47 as belonging to the thick disk (TH), and 94 as belonging to the halo (HL).

- We also proposed the orbit classification according to the GC orbital shapes in the X-Y, X-Z, and R-Z coordinate projections. Based on this classification, we found 110 GCs with tube orbits (TB), 8 GCs with perpendicular tube orbits (PT), 19 GCs with long radial orbits (LR), and 22 GCs with irregular orbits (IR).

- We analysed the evolution of the orbital parameters of the GCs in all five external MW-like potentials. The relative individual changes of the semi-major axis and the eccentricities are quite significant (up to 50%). The time evolution of these orbital elements clearly shows the separation between the inner and outer GCs. The inner GCs ($a \le 3$ kpc) have more regular and stronger eccentricity changes during the evolution. The outer GCs ($a > 3$ kpc) have much weaker eccentricity changes during the whole backward-integration time.

- We analysed the evolution of the apocenters and pericenters for GCs. As the result, we can conclude the more prominent dependences of the apocenter on the evolution of the MW-like potentials (masses and sizes). In case of pericenters, no strong variations were found.

- All our external potentials provide similar energy distributions for individual GCs. Comparing the GCs in different potentials at present, we can also conclude that the angular
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Evolution of the total angular momentum vs energy phase space ($L_{tot}$, *E*) (left panels, top to bottom), zth component of the angular momentum vs energy ($L_z$, *E*) (middle, top to bottom), and the perpendicular component of the angular momentum vs energy ($L_{perp}$, *E*) (right, top to bottom) for 411321, 441327, and 451323 TNG-TVPs. The colours represent data at different times in Gyr lookback time: magenta for $T = 0$, orange for $T = 2.5$, green for $T = 5$, cyan for $T = 7.5$, and red for $T = 10$. The GCs Pal 3, Crater, and Sagittarius II are not shown in the ($L_{tot}$, *E*) and ($L_{perp}$, *E*) phase spaces, and Sagittarius II is not shown in ($L_z$, *E*) phase space because their values are extremely high.</figcaption></figure>

momentum space is similar. Some differences were found in the energy distribution, which can be explained by different parameters of the galactic potentials (relative masses and spatial scales of stellar disks and DM halo). However, more differences (in individual GC energies) appeared when we studied the time evolution of the GC systems in the time-varying potentials.

- We individually checked the relative energy changes $\Delta E/E$ during the orbital evolution of GCs up to 10 Gyr. We identified 18 GCs for which the change exceeds 40%.

This strong change in the relative energies of these clusters also suggests that they might have a non-local origin (see also Malhan et al. 2022) (energy changes from 50% to 40% accordingly): Crater, Pal 3, Sagittarius II, Eridanus, AM 1, Pal 4, NGC 2419, Pyxis, Pal 14, Lae 3, Whiting 1, NGC 5694, Arp 2, Terzan 8, Terzan 7, NGC 7006, Pal 12, and Pal 13. We found that the in situ formed GCs are less strongly affected by the evolution of the TNG potentials than the clusters, which were likely formed in situ.
<figure><img src="image_17.png" /><figcaption>Fig. 15. Same as in Fig. 14, but for 462077 and 474170 TNG-TVPs.</figcaption></figure>

Our study clearly shows that in order to constrain the origin of the MW GCs better, it is vital to use time-variable Galactic potentials for the long-term (>3–4 Gyr) analysis of the GC orbital evolution.

**Acknowledgements.** The authors thank the anonymous referee for a very constructive report and suggestions that helped significantly improve the quality of the manuscript. MI, PB, CO and MM acknowledge support within the grant No. AP14869395 of the Science Committee of the Ministry of Science and Higher Education of Kazakhstan (‘Triune model of Galactic centre dynamical evolution on cosmological time scale’). The work of MI was supported by the Grant of the National Academy of Sciences of Ukraine for young scientists. MS acknowledges the support under the Fellowship of the President of Ukraine for young scientists 2022–2024. The work of PB, MI and MS was also supported by the Volkswagen Foundation under the grant No. 97778. The work of PB was supported by the Volkswagen Foundation under the special stipend No. 9B870 (2022). PB and MI also acknowledge support from the National Academy of Sciences of Ukraine under the Main Astronomical Observatory GPU computing cluster project No. 13.2021.MM. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

## References

Agertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64

Armstrong, B. M., Bekki, K., & Ludlow, A. D. 2021, MNRAS, 500, 2937

Bajkova, A. T., & Bobylev, V. V. 2021, Res. Astron. Astrophys., 21, 173

Bajkova, A. T., Carraro, G., Korchagin, V. I., Budanova, N. O., & Bobylev, V. V. 2020, ApJ, 895, 69

Bajkova, A. T., Smirnov, A. A., & Bobylev, V. V. 2021, Astron. Lett., 47, 454

Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957

Baumgardt, H., Hilker, M., Sollija, A., & Bellini, A. 2019, MNRAS, 482, 5138

Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417

Bennett, M., Bovy, J., & Hunt, J. A. S. 2022, ApJ, 927, 131

Berczik, P., Nitadori, K., Zhong, S., et al. 2011, in International Conference on High Performance Computing HPC-UA 2011, 8

Berczik, P., Spurzem, R., & Wang, L. 2013, in Third International Conference on High Performance Computing, HPC-UA 2013, 52

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Boldrini, P., & Bovy, J. 2022, MNRAS, 516, 4560

Bovy, J. 2011, Ph.D. Thesis, New York University, USA

Bovy, J., Allende Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131

Carpintero, D. D., & Aguilar, L. A. 1998, MNRAS, 298, 1

Chemerynska, I. V., Ishchenko, M. V., Sobolenko, M. O., Khoperskov, S. A., & Berczik, P. P. 2022, Adv. Astron. Space Phys., 12, 18

Conroy, C., Naidu, R. P., Garavito-Camargo, N., et al. 2021, Nature, 592, 534

Davies, J. J., Crain, R. A., Oppenheimer, B. D., & Schaye, J. 2020, MNRAS, 491, 4462

Donnari, M., Pillepich, A., Joshi, G. D., et al. 2021, MNRAS, 500, 4004

Drimmel, R., & Poggio, E. 2018, Res. Notes Am. Astron. Soc., 2, 210

Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716

Fattahi, A., Belokurov, V., Deason, A. J., et al. 2019, MNRAS, 484, 4471

Fernández-Trincado, J. G., Minniti, D., Souza, S. O., et al. 2021, ApJ, 908, L42

Ferrone, S., Di Matteo, P., Mastrobuono-Battisti, A., et al. 2023, A&A, 673, A44

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Garro, E. R., Minniti, D., Gómez, M., et al. 2022, A&A, 658, A120

Garrow, T., Webb, J. J., & Bovy, J. 2020, MNRAS, 499, 804

Genel, S., Nelson, D., Pillepich, A., et al. 2018, MNRAS, 474, 3976

Gómez, F. A., Besla, G., Carpintero, D. D., et al. 2015, ApJ, 802, 128

Grand, R. J. J., Deason, A. J., White, S. D. M., et al. 2019, MNRAS, 487, L72

Gratton, R., Bragaglia, A., Carretta, E., et al. 2019, A&A Rv., 27, 8

Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10

Habouzit, M., Genel, S., Somerville, R. S., et al. 2019, MNRAS, 484, 4413

Haghi, H., Zonoozi, A. H., & Taghavi, S. 2015, MNRAS, 450, 2812

Harris, W. E. 2010, ArXiv e-prints [arXiv:1012.3224]
Harris, W. E., Harris, G. L. H., & Alessi, M. 2013, ApJ, 772, 82

Ibata, R., Malhan, K., Martin, N., et al. 2021, ApJ, 914, 123

Ishchenko, M. V., Sobolenko, M. O., Kalambay, M. T., Shukirgaliyev, B. T., & Berczik, P. P. 2021, Rep. Nat. Academy Sci. Republic Kazakhstan, 6, 94

Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864

Katz, N., & Gunn, J. E. 1991, ApJ, 377, 365

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R. D. 2013, A&A, 558, A53

Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202344233

Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286

Laporte, C. F. P., Belokurov, V., Koposov, S. E., Smith, M. C., & Hill, V. 2020, MNRAS, 492, L61

Łokas, E. L. 2020, A&A, 638, A133

Lovell, M. R., Pillepich, A., Genel, S., et al. 2018, MNRAS, 481, 1950

Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141

Malhan, K., Ibata, R. A., Sharma, S., et al. 2022, ApJ, 926, 107

Mardini, M. K., Placco, V. M., Meiron, Y., et al. 2020, ApJ, 903, 88

Marín-Franch, A., Aparicio, A., Piotto, G., et al. 2009, ApJ, 694, 1498

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4

Mateu, C. 2023, MNRAS, 520, 5225

McConnachie, A. W. 2012, AJ, 144, 4

McConnachie, A. W., Higgs, C. R., Thomas, G. F., et al. 2021, MNRAS, 501, 2363

Merritt, D. 1999, PASP, 111, 129

Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533

Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2

Ness, M., & Lang, D. 2016, AJ, 152, 14

Pagnini, G., Di Matteo, P., Khoperskov, S., et al. 2023, A&A, 673, A86

Panithanpaisal, N., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 920, 10

Petersen, M. S., & Peñarrubia, J. 2021, Nat. Astron., 5, 251

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJ, 781, L31

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Sun, G., Wang, Y., Liu, C., et al. 2023, Res. Astron. Astrophys., 23, 015013

Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484, 5587

Valcin, D., Bernal, J. L., Jimenez, R., Verde, L., & Wandelt, B. D. 2020, J. Cosmol. Astropart. Phys., 2020, 002

VandenBerg, D. A., Brogaard, K., Leaman, R., & Casagrande, L. 2013, ApJ, 775, 134

van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35

Vasiliev, E. 2019, MNRAS, 484, 2832

Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978

Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874

Xiang, M., & Rix, H.-W. 2022, Nature, 603, 599

Yun, K., Pillepich, A., Zinger, E., et al. 2019, MNRAS, 483, 1042
**Appendix A: Evolution of the circular velocity at the distance of the Sun for TNG-TVPs.**

<figure><img src="image_18.png" /><figcaption>**Fig. A.1.** Evolution of the circular velocity at the distance of the Sun (R⊙ ~ 8.12 kpc; Gravity Collaboration 2019) in the Galactic disk as a function of backward-integration time in five TNG-TVPs, *first and third panels*. A straight line indicates the current velocity value of the solar rotation (V⊙ ~ 235 km s⁻¹; Mardini et al. 2020). Dotted green lines represent the original data from the IllustrisTNG-100. Dotted blue lines represent the interpolation and smoothing with a time step of 1 Myr. Second and fourth panels: Evolution of the TNG-TVP rotation curves at the initial (at present in red), average (5 Gyr in green), and final moments of time (10 Gyr in blue). *From left to right from top to bottom:* 411321, 441327, 451323, 462077, and 474170.</figcaption></figure>
Appendix B: Orbit types

<figure><img src="image_19.png" /><figcaption>**Fig. B.1.** Examples of different orbit types in three coordinate planes (X, Y), (X, Z), and (R, Z), where R is the planar galactocentric radius. Tube orbit (TB) of NGC 6717, perpendicular tube orbit (PT) of NGC 6715, long radial orbit (LR) of Palomar 14, and irregular orbit (IR) of NGC 6681, *from top to bottom*. The colour shows the evolution up to 10 Gyr. The blue dot is the initial location, and the red dot is the final location.</figcaption></figure>
## Appendix C: Globular clusters from the Gaia DR3 catalogue.

<table><caption>Table C.1. Globular clusters from the Gaia DR3 catalogue.</caption><thead><tr><th>ID</th><th>Name</th><th>TO</th><th>GR</th><th>ID</th><th>Name</th><th>TO</th><th>GR</th><th>ID</th><th>Name</th><th>TO</th><th>GR</th><th>ID</th><th>Name</th><th>TO</th><th>GR</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th><th>(13)</th><th>(14)</th><th>(15)</th><th>(16)</th></tr></thead><tbody><tr><td>1</td><td>NGC 104</td><td>TB</td><td>HL</td><td>41</td><td>NGC 5927</td><td>TB</td><td>TH</td><td>81</td><td>HP 1<sup>c</sup></td><td>IR</td><td>HL</td><td>121</td><td>NGC 6569</td><td>TB</td><td>TH</td></tr><tr><td>2</td><td>NGC 288</td><td>IR</td><td>HL</td><td>42</td><td>NGC 5946</td><td>TB</td><td>TH</td><td>82</td><td>FSR 1758</td><td>TB</td><td>HL</td><td>122</td><td>BH 261</td><td>TB</td><td>TH</td></tr><tr><td>3</td><td>NGC 362</td><td>LR</td><td>HL</td><td>43</td><td>NGC 5986</td><td>TB</td><td>TH</td><td>83</td><td>NGC 6362</td><td>TB</td><td>HL</td><td>123</td><td>NGC 6584</td><td>IR</td><td>HL</td></tr><tr><td>4</td><td>Whiting 1</td><td>PT</td><td>HL</td><td>44</td><td>FSR 1716</td><td>TB</td><td>TH</td><td>84</td><td>Liller 1</td><td>TB</td><td>BL</td><td>124</td><td>NGC 6624</td><td>TB</td><td>TH</td></tr><tr><td>5</td><td>NGC 1261</td><td>IR</td><td>HL</td><td>45</td><td>Pal 14</td><td>LR</td><td>HL</td><td>85</td><td>NGC 6380</td><td>TB</td><td>BL</td><td>125</td><td>NGC 6626</td><td>TB</td><td>TH</td></tr><tr><td>6</td><td>Pal 1</td><td>TB</td><td>HL</td><td>46</td><td>Lynga 7<sup>b</sup></td><td>TB</td><td>TH</td><td>86</td><td>Terzan 1</td><td>TB</td><td>TN</td><td>126</td><td>NGC 6638</td><td>TB</td><td>TH</td></tr><tr><td>7</td><td>AM 1<sup>a</sup></td><td>LR</td><td>HL</td><td>47</td><td>NGC 6093</td><td>TB</td><td>HL</td><td>87</td><td>Ton 2</td><td>TB</td><td>TH</td><td>127</td><td>NGC 6637</td><td>TB</td><td>TH</td></tr><tr><td>8</td><td>Eridanus</td><td>IR</td><td>HL</td><td>48</td><td>NGC 6121</td><td>TB</td><td>TH</td><td>88</td><td>NGC 6388</td><td>TB</td><td>TH</td><td>128</td><td>NGC 6642</td><td>TB</td><td>TN</td></tr><tr><td>9</td><td>Pal 2</td><td>LR</td><td>HL</td><td>49</td><td>NGC 6101</td><td>IR</td><td>HL</td><td>89</td><td>NGC 6402</td><td>TB</td><td>TH</td><td>129</td><td>NGC 6652</td><td>TB</td><td>HL</td></tr><tr><td>10</td><td>NGC 1851</td><td>LR</td><td>HL</td><td>50</td><td>NGC 6144</td><td>TB</td><td>HL</td><td>90</td><td>NGC 6401</td><td>TB</td><td>TH</td><td>130</td><td>NGC 6656</td><td>TB</td><td>HL</td></tr><tr><td>11</td><td>NGC 1904</td><td>LR</td><td>HL</td><td>51</td><td>NGC 6139</td><td>TB</td><td>TH</td><td>91</td><td>NGC 6397</td><td>TB</td><td>HL</td><td>131</td><td>Pal 8</td><td>TB</td><td>TH</td></tr><tr><td>12</td><td>NGC 2298</td><td>IR</td><td>HL</td><td>52</td><td>Terzan 3</td><td>TB</td><td>TH</td><td>92</td><td>Pal 6</td><td>TB</td><td>TH</td><td>132</td><td>NGC 6681</td><td>IR</td><td>HL</td></tr><tr><td>13</td><td>NGC 2419</td><td>IR</td><td>HL</td><td>53</td><td>NGC 6171</td><td>TB</td><td>TH</td><td>93</td><td>NGC 6426</td><td>TB</td><td>HL</td><td>133</td><td>NGC 6712</td><td>TB</td><td>HL</td></tr><tr><td>14</td><td>Pyxis</td><td>IR</td><td>HL</td><td>54</td><td>ESO 452-11</td><td>TB</td><td>TH</td><td>94</td><td>Djorg 1</td><td>TB</td><td>HL</td><td>134</td><td>NGC 6715</td><td>PT</td><td>HL</td></tr><tr><td>15</td><td>NGC 2808</td><td>TB</td><td>HL</td><td>55</td><td>NGC 6205</td><td>TB</td><td>HL</td><td>95</td><td>Terzan 5</td><td>TB</td><td>BL</td><td>135</td><td>NGC 6717</td><td>TB</td><td>TH</td></tr><tr><td>16</td><td>E 3</td><td>TB</td><td>HL</td><td>56</td><td>NGC 6229</td><td>LR</td><td>HL</td><td>96</td><td>NGC 6440</td><td>TB</td><td>BL</td><td>136</td><td>NGC 6723</td><td>TB</td><td>HL</td></tr><tr><td>17</td><td>Pal 3</td><td>IR</td><td>HL</td><td>57</td><td>NGC 6218</td><td>TB</td><td>TH</td><td>97</td><td>NGC 6441</td><td>TB</td><td>TH</td><td>137</td><td>NGC 6749</td><td>TB</td><td>TH</td></tr><tr><td>18</td><td>NGC 3201</td><td>TB</td><td>HL</td><td>58</td><td>FSR 1735</td><td>TB</td><td>TH</td><td>98</td><td>Terzan 6</td><td>TB</td><td>BL</td><td>138</td><td>NGC 6752</td><td>TB</td><td>TH</td></tr><tr><td>19</td><td>Pal 4</td><td>LR</td><td>HL</td><td>59</td><td>NGC 6235</td><td>TB</td><td>HL</td><td>99</td><td>NGC 6453</td><td>TB</td><td>TH</td><td>139</td><td>NGC 6760</td><td>TB</td><td>HL</td></tr><tr><td>20</td><td>Crater</td><td>IR</td><td>HL</td><td>60</td><td>NGC 6254</td><td>TB</td><td>TH</td><td>100</td><td>UKS 1</td><td>TB</td><td>HL</td><td>140</td><td>NGC 6779</td><td>TB</td><td>HL</td></tr><tr><td>21</td><td>NGC 4147</td><td>LR</td><td>HL</td><td>61</td><td>NGC 6256</td><td>TB</td><td>TN</td><td>101</td><td>VVV CL001</td><td>TB</td><td>TN</td><td>141</td><td>Terzan 7</td><td>PT</td><td>HL</td></tr><tr><td>22</td><td>NGC 4372</td><td>TB</td><td>TH</td><td>62</td><td>Pal 15</td><td>LR</td><td>HL</td><td>102</td><td>NGC 6496</td><td>TB</td><td>TH</td><td>142</td><td>Pal 10</td><td>TB</td><td>HL</td></tr><tr><td>23</td><td>Rup 106</td><td>TB</td><td>HL</td><td>63</td><td>NGC 6266</td><td>TB</td><td>TN</td><td>103</td><td>Terzan 9</td><td>TB</td><td>BL</td><td>143</td><td>Arp 2</td><td>PT</td><td>HL</td></tr><tr><td>24</td><td>NGC 4590</td><td>TB</td><td>HL</td><td>64</td><td>NGC 6273</td><td>TB</td><td>HL</td><td>104</td><td>Djorg 2</td><td>TB</td><td>BL</td><td>144</td><td>NGC 6809</td><td>TB</td><td>HL</td></tr><tr><td>25</td><td>BH 140</td><td>TB</td><td>HL</td><td>65</td><td>NGC 6284</td><td>TB</td><td>HL</td><td>105</td><td>NGC 6517</td><td>TB</td><td>TN</td><td>145</td><td>Terzan 8</td><td>PT</td><td>HL</td></tr><tr><td>26</td><td>NGC 4833</td><td>TB</td><td>HL</td><td>66</td><td>NGC 6287</td><td>TB</td><td>HL</td><td>106</td><td>Terzan 10</td><td>TB</td><td>HL</td><td>146</td><td>Pal 11</td><td>TB</td><td>HL</td></tr><tr><td>27</td><td>NGC 5024</td><td>IR</td><td>HL</td><td>67</td><td>NGC 6293</td><td>TB</td><td>HL</td><td>107</td><td>NGC 6522</td><td>TB</td><td>TH</td><td>147</td><td>Sagittarius II</td><td>PT</td><td>HL</td></tr><tr><td>28</td><td>NGC 5053</td><td>IR</td><td>HL</td><td>68</td><td>NGC 6304</td><td>TB</td><td>TN</td><td>108</td><td>NGC 6535</td><td>TB</td><td>TH</td><td>148</td><td>NGC 6838</td><td>TB</td><td>HL</td></tr><tr><td>29</td><td>NGC 5139</td><td>TB</td><td>HL</td><td>69</td><td>NGC 6316</td><td>TB</td><td>TH</td><td>109</td><td>NGC 6528</td><td>TB</td><td>TN</td><td>149</td><td>NGC 6864</td><td>TB</td><td>HL</td></tr><tr><td>30</td><td>NGC 5272</td><td>IR</td><td>HL</td><td>70</td><td>NGC 6341</td><td>LR</td><td>HL</td><td>110</td><td>NGC 6539</td><td>TB</td><td>TH</td><td>150</td><td>NGC 6934</td><td>LR</td><td>HL</td></tr><tr><td>31</td><td>NGC 5286</td><td>TB</td><td>HL</td><td>71</td><td>NGC 6325</td><td>TB</td><td>TH</td><td>111</td><td>NGC 6540</td><td>TB</td><td>TH</td><td>151</td><td>NGC 6981</td><td>LR</td><td>HL</td></tr><tr><td>32</td><td>AM 4</td><td>PT</td><td>HL</td><td>72</td><td>NGC 6333</td><td>TB</td><td>HL</td><td>112</td><td>NGC 6544</td><td>TB</td><td>TH</td><td>152</td><td>NGC 7006</td><td>LR</td><td>HL</td></tr><tr><td>33</td><td>NGC 5466</td><td>IR</td><td>HL</td><td>73</td><td>NGC 6342</td><td>TB</td><td>TH</td><td>113</td><td>NGC 6541</td><td>TB</td><td>TH</td><td>153</td><td>Laevens 3</td><td>IR</td><td>HL</td></tr><tr><td>34</td><td>NGC 5634</td><td>IR</td><td>HL</td><td>74</td><td>NGC 6356</td><td>TB</td><td>HL</td><td>114</td><td>2MASS-GC01</td><td>TB</td><td>TN</td><td>154</td><td>NGC 7078</td><td>TB</td><td>HL</td></tr><tr><td>35</td><td>NGC 5694</td><td>LR</td><td>HL</td><td>75</td><td>NGC 6355</td><td>TB</td><td>TH</td><td>115</td><td>ESO 280-06</td><td>TB</td><td>HL</td><td>155</td><td>NGC 7089</td><td>LR</td><td>HL</td></tr><tr><td>36</td><td>IC 4499</td><td>IR</td><td>HL</td><td>76</td><td>NGC 6352</td><td>TB</td><td>TH</td><td>116</td><td>NGC 6553</td><td>TB</td><td>TH</td><td>156</td><td>NGC 7099</td><td>TB</td><td>HL</td></tr><tr><td>37</td><td>NGC 5824</td><td>IR</td><td>HL</td><td>77</td><td>IC 1257</td><td>LR</td><td>HL</td><td>117</td><td>2MASS-GC02</td><td>TB</td><td>TH</td><td>157</td><td>Pal 12</td><td>PT</td><td>HL</td></tr><tr><td>38</td><td>Pal 5</td><td>IR</td><td>HL</td><td>78</td><td>Terzan 2</td><td>TB</td><td>BL</td><td>118</td><td>NGC 6558</td><td>TB</td><td>TN</td><td>158</td><td>Pal 13</td><td>IR</td><td>HL</td></tr><tr><td>39</td><td>NGC 5897</td><td>TB</td><td>HL</td><td>79</td><td>NGC 6366</td><td>TB</td><td>TH</td><td>119</td><td>IC 1276<sup>d</sup></td><td>TB</td><td>HL</td><td>159</td><td>NGC 7492</td><td>LR</td><td>HL</td></tr><tr><td>40</td><td>NGC 5904</td><td>LR</td><td>HL</td><td>80</td><td>Terzan 4</td><td>TB</td><td>BL</td><td>120</td><td>Terzan 12</td><td>TB</td><td>TH</td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** Columns (1), (5), (9), and (13) represent the GC index number. Columns (2), (6), (10), and (14) represent the names of all GCs in Vasiliev & Baumgardt (2021). The alternative names of the GCs, which were used in our previous works Ishchenko et al. (2021) and Chemerynska et al. (2022): <sup>a</sup>E 1, <sup>b</sup>BH 184, <sup>c</sup>BH 229, and <sup>d</sup>Pal 7. Columns (3), (7), (11), and (15) show the type of GC orbits (TO), where the tube orbit is abbreviated as TB, the perpendicular tube orbit as PT, the long radial orbit as LR, and the irregular orbit as IR. Examples for the each type are given in Appendix A. Columns (4), (8), (12), and (16) show GCs belonging to the different regions of the Galaxy (GR) at present, according to the classification from Bland-Hawthorn & Gerhard (2016), where the bulge is abbreviated as BL, the thin disk as TN, the thick disk as TH, and the halo as HL.