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

Context. Recent observational data show that the Milky Way (MW) galaxy contains about 170 globular clusters (GCs). A fraction of them is likely formed in dwarf galaxies accreted onto the MW in the past, while the remaining of clusters are formed in-situ. Therefore, different parameters, including orbits, of the globular clusters is a valuable tool for studying the Milky Way evolution. However, since the evolution of the 3D mass distribution of the MW is poorly constrained, the orbits of the clusters are usually calculated in static potentials. Aims. In this work, we study the evolution of the GCs in several external potentials, where we aim to quantify the effects of the evolving galaxy potential on the orbits of the GCs. Methods. For the orbits calculation we used five MW-like potentials from IllustrisTNG-100 simulation. The orbits of 159 GCs were integrated using a high-order N-body parallel dynamic code phi-GPU, with initial conditions obtained from recent Gaia DR3 catalogues. Results. We provide a classification of the GCs orbits according to their 3D shapes and association with different components of the MW (disk, halo, bulge). We also found that the globular clusters in the external potentials have roughly similar energy-angular momentum distributions at the present time. However, both total energy and total angular momentum of the GCs are not 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 GCs are less affected by the evolution of the TNG potentials as compared to the clusters which are likely formed ex-situ. Therefore, our results suggest that time-varying potentials significantly affect the orbits of the GC, thus making it vital for understanding the formation of the MW.


Introduction
According to the ΛCDM model, the Milky Way (MW) globular clusters (GCs) are the first stellar associations formed in the early Universe as gravitationally bound systems (Gratton et al. 2019).As the result their typical ages are > 10 − 12 Gyr (Marín-Franch et al. 2009;VandenBerg et al. 2013;Valcin et al. 2020) and the present masses 10 5 M (Harris et al. 2013;Kharchenko et al. 2013;Baumgardt et al. 2019;Baumgardt & Vasiliev 2021).Many GCs' survived over long galactic history and, thus, became relics from the earliest era of galaxy formation (Garro et al. 2022).Recent observations reveal the picture where 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 the remnants of either GCs or dwarf galaxies destroyed in the process of "merging history" of our Galaxy (Mateu 2023;Ibata et al. 2021;Ferrone et al. 2023).They undergo complete destruction (or mass loss) under the influence of the complex time-dependent gravitational field of our Galaxy.Study of kinematics and chemical abundances of the GCs captured by the MW is a very good tool for studying the past history of our Galaxy (Xiang & Rix 2022), as well as searching possible progenitors for these GCs (Massari et al. 2019;Malhan et al. 2022).
Thanks to the publication of the Gaia (ESA) catalogues (Gaia Collaboration et al. 2021Collaboration et al. , 2018) full phase-space information is available for almost the entire population of the GCs, which makes possible to investigate their orbits.Several recent works studied the orbits of the GCs in the MW potential assuming its being axisymmetric and non-evolving 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 the z-component of the angular momentum (L z ) of the GCs, which is being used to group some of the GCs and link these kinematically coherent groups with their galaxy-progenitors (either MW or accreted satellite).However, the above-mentioned assumptions are not fully correct because the MW galaxy con-tains the bar and boxy-peanut bulge in the center (Dwek et al. 1995;Wegg & Gerhard 2013;Ness & Lang 2016).In the outer parts, the axisymmetry of the MW potential is believed to be broken by the massive orbiting satellites (Laporte et al. 2018(Laporte et al. , 2020)), including the LMC/SMC 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 its rapid increase at early times (until ≈8 Gyr) followed by a slower evolution (Snaith et al. 2014) which is in agreement with the mass growth of the MW-like galaxies constrained via halo abundance matching (van Dokkum et al. 2013;Fattahi et al. 2019).While some models suggest that the mass-growth of the MWtype galaxies affects the coherency 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. 2022;Pagnini et al. 2022), 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 a time variable external potential.This feature significantly differs our investigation from other similar works (see for example 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), where the authors use the fixed or the analytically time-evolving MW potential.In series of articles about evolution of GC, and in particular this one, we aim to investigate how the evolution of the mass and spatial scales of both stellar and dark matter components of the MW affect the orbital parameters of the GCs.Since the evolution of the MW parameters is poorly constrained, we analyse the dynamics of the GCs in the MW-like evolving potentials directly obtained from IllustrisTNG cosmological simulations (Pillepich et al. 2018;Nelson et al. 2018).
The paper is organised as follows.In Section 2.1 we present the GCs selection with the error distribution from the Gaia data release 3 (DR3) catalogue.In Sections 2.2 and 2.3, we introduce the time-dependent gravitational potentials of the MW-type galaxies that were selected from IllustrisTNG-100 and shortly describe the numerical integrator.In Section 3, we present the individual orbits of the GCs in the evolving TNG potentials.In Section 4 we describe GCs phase space evolution.In Section 5 we summarise our results.

Globular Clusters sample
In our study, we used the recent catalogues of the 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 phasespace information required for the initial conditions of our simulations: right ascension (RA), declination (DEC), heliocentric distance (D ), 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 use also the self-gravity of the clusters together with the Galaxy external potential.We exclude the Mercer 5 object since there is no mass information for this GC from the catalogues above.
Before the orbital integration, we analysed the GCs Gaia measurement errors influence which may change the orbital shape during the orbital integration (see Bajkova & Bobylev 2021).We analysed the errors in PMRA, PMDEC, RV and D for 159 GCs from the catalogue above and in Fig. 1 presented the distribution of the relative errors of listed parameters.As we see, the vast majority (up to 95%) of the errors are well below 5% but some GCs (Crater and Lae 3) have a quite large value in proper motion.Analysing the errors in the heliocentric distance we found only 18 objects that have 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 subsection 3.1.
For the transformation of the positions and velocities to the Cartesian Galactocentric rest-frame (see Johnson & Soderblom 1987;Bovy 2011, for the coordinate transformation equations), we assumed the Galactocentric distance of the Sun R = 8.178 kpc (Gravity Collaboration et al. 2019;Reid & Brunthaler 2004), a height above the Galactic plane Z = 20.8pc (Bennett & Bovy 2019) and the velocity of the Local Standard of Rest (LSR) V LSR = 234.737km s −1 (Bovy et al. 2012;Drimmel & Poggio 2018).Accordingly, the Sun is located at X = −8178 pc, Y = 0 pc and Z = 20.8pc in our Cartesian Galactocentric coordinate system.For the peculiar velocity of the Sun we used the following values with respect to the LSR: U = 11.1 km s −1 , V = 12.24 km s −1 , W = 7.25 km s −1 (Schönrich et al. 2010).

Time-variable gravitational potentials from Illustris simulation
Since the time evolution of the MW potential is poorly constrained from observations, for the GCs orbits integration, we have selected the MW-like galaxies from the most recent Illustris cosmological numerical simulations.For our purpose, we used the publicly available snapshot data from the IllustrisTNG-100 (Pillepich et al. 2018;Nelson et al. 2018;Springel et al. 2018;Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2019).These simulations successfully reproduce the main scaling relations and evolution of galaxy sizes (Genel et al. 2018), including the formation of realistic disc galaxies (Pillepich et al. 2019), the gas-phase mass-metallicity relation (Torrey et al. 2019), and the range of MW-type galaxies observed kinematic properties (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), central black hole activity (Habouzit et al. 2019;Donnari et al. 2021), star formation quenching (Genel et al. 2018;Davies et al. 2020), and many others.These results ensure that the investigation of properties of different types of galaxies from the Illustris cosmological simulations allows us to study the GCs dynamics in an environment similar to the MW.The IllustrisTNG-100 is characterised by a simulation box approximately 100 Mpc 3 .In a box of such size, one can resolve a sufficient number of the MW-mass disk galaxies (in our sample 54 object 1 ) with the mass resolution of 7.5 × 10 6 M for dark matter and 1.4 × 10 6 M for the baryonic particles, respectively.For our analysis, in the Illustris simulations, we identified the MW-like galaxy candidates with at least 10 5 dark matter particles and at least 10 3 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 at redshift zero well reproduce MW parameters (disk and halo masses with the spatial scales) at present the best (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, simply we do not have enough particles to determine such a 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 extra parameter to select the best TNG galaxies which represent the MW-type systems.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 −1 (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).Due to the galaxy mass decrease as a function of lookback time, we clearly see the circular velocity decrease at all galactocentric radii.
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) po-tentials: where R = x2 + y 2 is the planar Galactocentric radius, z is the distance above the plane of the disc, G is the gravitational constant, a d is the disk scale length, b d,h are the disk and halo scale heights, respectively, M d and M h = 4πρ 0 b 3 h (ρ 0 is the central mass density of the halo) are masses of the disk and halo, respectively.
In our case all these mass distribution parameters are functions of lookback time.We performed such a fit for each available snapshot in each of our selected galaxies.The code, which we use to find these parameters, is also publicly available at GitHub 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 the timestep of 1 Myr.Finally, to avoid some short-time scale perturbations, likely linked to the mergers or errors in the galaxy identification, we additionally smoothed the evolution of the main galaxy parameters.The above procedure results are shown for our five time-variable potentials (TNG-TVPs) in Fig. 2. Table 1.Parameters of the time-varying potentials selected from the IllustrisTNG-100 simulation at redshift zero.The last column shows the parameters of the corresponding MW components according to Bennett et al. (2022)

Integration procedure
For the GCs' orbits integration we used a high-order parallel dynamical N-body code ϕ−GPU3 (Berczik et al. 2011(Berczik et al. , 2013)).The code is based on the fourth-order Hermite integration scheme with hierarchical individual block time steps.For the force calculation, acting on the particles during the integration, we combine the particles self gravity with the time variable external potential described in the previous subsection.For the self gravity calculation (as a first order simplification) we apply the current mass of the individual GCs.We applied the integration timestep parameter η=0.01 (Makino & Aarseth 1992), which gives for our investigation the required integration accuracy.As a test case, we run a typical simulation with the fixed for today (that is, for t = 0) values of the external potential #411321 with up-to-date halo and disk masses and scale lengths.During this test case the total relative energy drift (∆E/E t=0 ) over a 10 Gyr backward integration was below ≈ 2.5 × 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 for position and velocity the central values from Baumgardt & Vasiliev (2021) and Vasiliev & Baumgardt (2021).

Influence of the measurement errors
Prior to the analysis of the orbital types of the GCs, we investigate the impact of the initial conditions uncertainties on the initial velocities values and structure of the orbits.For each GC in our sample, we generated extra 10 initial conditions where the initial velocities take into account the normal distribution of the errors for PMRA, PMDEC, RV (blue circles), and D (green circles) from Baumgardt & Vasiliev (2021) 4 , see Fig. 3.As we see, both randomisations have a similar effect on the GCs' initial velocities.In general, the errors by velocities are even slightly dominated over the errors due to the distance determination.
The effect from the velocity determination errors is even slightly larger 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 proper motions and radial velocity) as for the NGC 6121 (with the one of the smallest values: 0.1%, 0.12%, and 2% respectively).These two GCs are marked by crosses in Fig. 3.The velocity uncertainties significantly change the orbital characteristics 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 vertical motions are less affected.In contrast, for the NGC 6121, which parameters are well measured, the orbits are less affected by the velocity uncertainties, as we can see on the bottom row of Fig. 4. Thus, since the velocity uncertainties for most of our GCs are similar to the NGC 6121, we assume that the orbits of the majority of the GCs are robust to the variations of the initial condition.
As we can see in Fig. 4 the influence of the relative errors of the heliocentric distances on the orbital shape of our GC sample is indeed less significant compared to 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).As we can see, the effect from the relative distance error is significantly smaller compared to the effect from the proper motions and radial velocity errors, presented on the same plot.Next, we focus our attention only on the GC-point orbits obtained using the central values for positions and velocities from Baumgardt & Vasiliev (2021) and Vasiliev & Baumgardt (2021).

The evolution of the Globular Clusters orbital elements
Assuming the cosmologically motivated time variable potential we automatically taking into account the possible influence of the Milky Way satellite galaxies on to the GCs subsystem.If we compare the individual GCs' orbital elements (semi-major axis a and eccentricityecc) variation over the whole integration time (sometimes over the hundreds or thousands of GCs' orbital revelations).In Fig. 5 we present the orbital elements relative changes in five TNG-TVPs.In Fig. 6 we have shown the time evolution of the a and ecc.The strong influence on the GCs' orbits have already discussed in details in the papers Garrow et al. (2020) and Boldrini & Bovy (2022).Our Fig. 5 and Fig. 6 are directly comparable with the Figures 1, 2 and 3 from the Garrow et al. (2020).As we can see the relative changes of a and ecc behave in a similar range.In our case the a evolution is even larger in a factor of two.We also notice a very similar distribution of the relative quantities for the different TNG-TVP potentials.
As a illustration of the time evolution of orbital elements of all GCs we have shown such a data for the #411321 TNG-TVP in Fig. 6.Here we can clearly see the separation between the 'inner' and 'outer' GCs.The 'inner' GCs (a 3 kpc) have a 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.This Fig. 6 can be directly compared with the Figure 3    GCs, even with large a during the whole time of integration -10 Gyr.
For example, in the works of Garrow et al. (2020) and Boldrini & Bovy (2022) during their GCs orbit investigations the MW major satellites were added as an extra gravitational perturbation in the fixed MW potential approximation.In contrast in our case, the IllustrisTNG-100 mass assembling history of our theoretical MW-like galaxies (see Fig. 2) provides 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 GCs' main orbital parameters evolution: apocenters and pericenters (as for a similar plot one can look of the work Bajkova et al. 2021, see Figure 9).On the Xaxis we plotted the apo and peri values of each GCs in the static (fixed) potential, that is, using 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) are quite similar.Up to roughly 5 Gyr lookback time we we did not notice significant changes in these quantities.But af-  ter approximately 8 Gyr the differences (especially in apocenters values) are 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 the #411321 TNG-TVP.For the visualisation, we selected several GCs with different heliocentric distances.In our representative sample, the Liller and Terzan 4 are the closest GCs to the Galactic centre.
As we see in Fig. 8, the Liller and Terzan 4 as closest GCs have a clear decreasement of the apocenter with the minimum values at the 8 Gyr with further increasement.The NGC 6717 has a slow decline with the minimum at the same time.Meanwhile, we don't see such a clear picture of pericenter evolution.
Such a behaviour can be understood taking into account the more strong dependence of the apocenter value (compared to the pericenter) from the mass and size changes of our TNG-TVPs.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.The evolution of the apocenters and pericenters for all GCs we present on the web-site5 , for all five TNG-TVPs and for each GC separately.
Also, we present in Fig. 9 the comparison of apocenters and pericenters of 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 more strong influence of the time variable masses and sizes of MW-like potentials to the orbits.

Types of orbits
We calculated the orbits of 159 GCs lookback time for 10 Gyr in each of our TNG-TVP potentials.The orbital evolution for our GC sample is presented online 5 .The visualisation was carried out in three projections of the Cartesian Galactocentric restframe: (X, Y), (X, Z), and (R, Z), where R is the planar Galactocentric radius.As we see, 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, etc., depending on specifics of the potential (e.g., oblate/prolate) and the particular combination of 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 divide them into four main types (categories): -Tube orbit (TB) -110 (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 in the (R, Z) has a trapezoidal shape.-Perpendicular tube orbit (PT) -8 (5%) orbits.The orbital plane is close to the meridional, also contains a hole in the centre and in the (R, Z) plane has a trapezoidal shape.-Long radial orbit (LR) -19 (12%) orbits.These are longperiod, near-hyperbolic orbits characterised by large eccentricity with no prominent hole in the centre.-Irregular orbit (IR) -22 (14%) orbits.The orbit can not be classified as one of the described above.
We listed the orbit types for each GC in Table C.1.In Fig. B.1 we show some typical examples of the orbits used in the classification above.

Association with the regions of the Galaxy
We assumed that Galaxy has several spatial components: bulge, thin disk, thick disk and main halo.Each GCs' orbit can be associated predominantly with one of them (based on the distance criteria).Such 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 at the past We associated nine GCs. -Thin disk (TN).Scale height in Z directions -∼0.3 kpc.Scale length in R plane -∼2.6 kpc.With such a scales we determine the boundary dimensions for the thin disk as ∼7.8 × ∼0.9 kpc.We associated ten GCs.-Thick disc (TH).Scale height in Z directions -∼0.9 kpc.
Scale length in R plane -∼2 kpc.With such a scales we determine the boundary dimensions for the thick disc as ∼6 × ∼2.7 kpc.We associated 47 GCs.-Halo (HL).Over then ∼2.7 kpc in Z direction and more than ∼7.8 kpc in R plane.We associated 94 GCs.
As we see from Fig. 10 most GCs at the present stage (last billion year) belong to the halo and thick disk, at 59% and 30%, respectively.The bulge and thin disk contain just 5% and 6% GCs.But in the distant past there is a high probability 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 listed 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 made, in our point of view, these potentials quite representative.
The majority of the GCs from the bulge, thin and thick have tube orbits, see Table C. Possibly, it underwent extreme encounter(s) with another GC, which changed its orbit.In contrast, halo GCs display all types of orbits, with the most common being tube, irregular and long radial orbits.Also, the halo harbours objects with all GCs that have perpendicular tube orbits.

Energy changes of Globular Clusters during evolution
Since we consider 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 ∆E/E during the evolution of GCs over 10 Gyr.In Fig. 11 we present the relative energy changes of individual GCs during the orbit integration.In general, we see the effect of Galaxy mass changes in our selected TNG-TVPs (Fig. 2).We found that, the maximum rela-tive energy changes are ∼45-50%.As an example, in Fig. 11 we highlighted five GCs in colour: BH 140 and VVV CL001 (red and blue colours) associated with the Galactic disk and likely formed in-situ; AM 4, Sagittarius II, and Lae 3 associated with the halo and likely born outside the Galaxy (see the definition of association Malhan et al. 2022;Fernández-Trincado et al. 2021).All other 154 GCs are marked with a grey colour.As we can see from Fig. 11, the in-situ GCs are less affected by the changes in the gravitational potential compared to the accreted clusters.This trend is especially visible on the last panel of the Fig. 11 for the #474170 TNG-TVP.
All of them belong to the halo (HL) and, in most cases they have either irregular (IR) or perpendicular tube-type (PT) of orbits, see Table C.1.

Present-day properties of the Globular Clusters
In our study we analyse the GCs' 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 GCs' mass.In particular, we analyse the phase-space distribution of the GCs at present in three phase space coordinates combinations: total angular momentum (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 = L 2 x + L 2 y + L 2 z and L perp = L 2 x + L 2 y .The orientation of the Z-axis in our Cartesian Galactocentric coordinates is directed towards to the Galactic South Pole.By this definition 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, orbiting close to the Galactic center since their total energies are largest (by absolute value).Most of the GCs with TB obits show from zero to some prograde rotation with the smallest perpendicular angular momentum component.It can indicate their common dynamical origin or possible orbital evolution as a consequence of their similar angular momentum losses timescales.
Globular clusters with IR orbits show a wide range of the total energy and the angular momentum components and overlap with the GCs of other types of orbits.Their origin is unlikely the same based on the quite wide range of their angular momentum.
The GCs with LR orbits have a wide range of the total energies but tend to have smaller absolute values of the angular momentum.Interestingly these GCs have a relatively narrow range of L perp with a weak trend to have a retrograde rotation.
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 large values of L perp might suggest their common origin.As it was mentioned 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-codded 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 GCs' parameters only in the #411321 TNG-TVP.As we can see from the top panel in Fig. 13, as expected, the GCs that belong to the bulge (BL) and thick disk (TH) have the maximum negative energy values in the entire sample.But the bulge GCs have a relatively narrow energy range compared to the thick disk GCs.The GCs associated with the halo form a system which is less bound with the Galaxy.They have wide energy and angular momentum ranges.
In Fig. 13 we present a comparison of the GCs' positions in the phase space for other four #441327, #451323, #462077, and #474170 TNG-TVPs, bottom panel.The colours (as on 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 be easily explained with different parameters of the galactic potentials (see Table 1) at the present.At the same time, we can note that our sample of external potentials does not give any strong mutual discrepancies.

The Globular Cluster subsystem phase space evolution in five time-varying potentials
Using the time-evolving potentials, we can explore the effect of the Galactic parameters time-variation on the GCs' phase space distributions in energy and angular momentum.In Figs. 14 and 15 we show the GC systems snapshots for different times (0.0, 2.5, 5.0, 7.5, and 10 Gyr lookback time) for each of our TNG-TVPs.
As we can see from Figs. 14 and 15, all our external potentials provide roughly similar phase space energy values for individual GCs.The difference between potentials becomes more defined during 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 (that is, their bounding energies are larger; magenta) and at the end of the simulation (10 Gyr, red) their energies are more positive -this would mean GCs are less bounded which is not the case.Our expectations are fully met by the potential #441327 (Fig. 14, second row).Other TNG-TVP external potentials show slightly different behaviour.For example, the #451323 -(Fig.14, third row), #462077 -(Fig.15, first row) external potentials have maximum bounding energies at 5 Gyr (green).We can explain such behaviour if we check the evolution of the mass and characteristic scales from Fig. 2. For these two MW-like potentials around 5 Gyr, we can observe higher mass and smaller scale sizes (as for the halo and disk).
For the potential #411321 -(Fig.14, first row) the energy evolution shows a more complex behaviour.First the energy values for GCs are going down (approximately 3 Gyr; orange and approximately 5 Gyr; green) and only after ≈ 7 Gyr (cyan) they starts to raise.Such a complex behaviour might be also under-  stood based on the evolution of the halo and disk masses visible in Fig. 2 for this TNG-TVP.
In a case of the #474170 TNG-TVP (Fig. 15, second row).In contrast to the previous external potentials here in the beginning (up to approximately 3 Gyr; magenta-orange) we have almost a constant halo mass and characteristic scale.
In general, we conclude that the evolution of the GCs' phase space distributions weakly varies in the five selected TNG-TVPs.The differences are obviously caused by slightly different accretion history, which originally come from the IllustrisTNG-100 simulation.

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 GCs, we used the values from the catalogues (Vasiliev & Baumgardt 2021), which 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 a 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 Milky Way parameters at present.To obtain the spatial scales of the disks and dark matter halos, we decompose the mass distribution using the MN and NFW potentials.Also, to verify the similarity of the chosen TNG potentials to our Galaxy, we added an extra parameter -the observed circular velocity of rotation at the solar radius.Finally, we selected the five TNG galaxies #411321, #441327, #451323, #462077, and #474170 as time-variable external potentials for our study.
The main results of this study are as follows.
-We integrated the orbits of 159 MW GCs up to 10 Gyr evolution in five external potentials.According to our results, we select a set of GCs which (at present time) belong to four main structural components of the Galaxy.We define nine GCs as currently belonging to the bulge (BL), ten -to the thin disk (TN), 47 -to the thick disc (TH), and 94 -to the halo (HL).-We also proposed the orbit classification according to the GCs' 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), eight 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 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 ≤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.
-We analysed the evolution of the apocenters and pericenters for GCs and as a result, we can conclude the more prominent dependencies of the apocenter from the evolution of the MWlike potentials (masses and sizes).In case of pericenters, we don't see such strong variations.-All our external potentials provide similar energy distributions for individual GCs.Comparing the GCs in different potentials at present, we can also conclude that there are no significant differences in the angular momentum space.Some differences were found in the energy distribution, which can We found that the in-situ formed GCs are less affected by the evolution of the TNG potentials as compared to the clusters, which are likely formed ex-situ.
Our study clearly shows that in order to constrain better the origin of the MW GCs it is vital to use time-variable Galactic potentials for the long-term (> 3 − −4 Gyrs) analysis of the GC orbital evolution.

Fig. 1 .
Fig. 1.Distribution of the relative errors (from left to right): for the radial velocity (eRV), proper motions in right ascension (ePMRA) and declination (ePMDEC), and in the heliocentric distance eD .Blue histograms represent the error distributions.The following GCs are not shown due to the large values of these errors: Crater, Lae 3, NGC 6760, BH 261, NGC 6553, and Pal 13, except the graph on the bottom right panel.

Fig. 2 .
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.Halo mass M h , disk mass M d , NFW halo scale parameter b h , MN disk scale parameters a d and b d are presents 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 a 1 Myr time step that was used in the orbital integration.
from theGarrow et al. (2020).As we can see, our potential keep bound the individual

Fig. 3 .
Fig. 3. Influence of error randomisation to the initial velocity components in the Galactocentric Cartesian reference frame.Distance determination -green circles and velocity determination -blue circles.Red cross marked Crater, black -NGC 6121.

Fig. 4 .
Fig. 4. Changes in the orbits caused by the uncertainties in the radial velocity and proper motions 4 for Crater and NGC 6121.The presented orbits were integrated for 10 Gyr lookback time.The different colours with transparency are represent the orbits for five random realisations of the initial conditions in the #411321 TNG-TVP.

Fig. 6 .
Fig.6.Evolution of the GCs' orbital semi-major and eccentricity during the whole backward integration time in the case of #411321 TNG-TVP.We can clearly see the separation between the 'inner' and 'outer' GCs.The 'inner' GCs (a ≤3 kpc) have a more regular and larger eccentricity changes during the evolution.The 'outer' GCs (a >3 kpc) have a much smaller eccentricity changes during the whole backward integration time.

Fig. 7 .
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.

Fig. 8 .
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 .

Fig. 12 .
Fig. 12. Distribution of orbital types at redshift zero for the #411321 TNG-TVP.From left to right: total angular momentum versus energy (L tot , E), z-th component of the angular momentum versus total energy (L z , E), the perpendicular component of the angular momentum versus total energy (L perp , E).The colour corresponds to the type of the GCs' orbit: blue -tube orbit (TB), orange -perpendicular tube orbit (PT), greenlong radial orbit (LR), and red -irregular orbit (IR).Pal 3, Crater, and Sagittarius II are not shown in (L tot , E) and (L perp , E) phase spaces, also Sagittarius II is not shown in the (L z , E) phase space due to their extremely large values.

Fig. 13 .
Fig. 13.Position of GCs in different regions of the Galaxy at the present time for the #411321 TNG-TVP, top panel.From left to right: total angular momentum, z-th component of the angular momentum and perpendicular component of the angular momentum vs energy.Colour coding represents the GCs that belong to the different regions in the Galaxy: orange -bulge (BL), blue -thin disk (TN), green -thick disk (TH), and redhalo (HL).Bottom panel: the 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 also Sagittarius II is not shown in (L z , E) phase space due to of their extremely high values.

Fig. 14 .
Fig. 14.Evolution of the total angular momentum vs energy phase space (L tot , E) (left panels, top to bottom), z-th 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 is T = 0, orange is T = 2.5, green is T = 5, cyan is T = 7.5 and red is 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 also Sagittarius II is not shown in (L z , E) phase space due to of their extremely high values.

Fig. A. 1 .
Fig. A.1.Evolution of the circular velocity at the distance of the Sun (R ≈ 8.12 kpc; Gravity Collaboration et al. 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 Sun rotation (V ≈ 235 km s −1 ; Mardini et al. 2020).Green dotted lines represent the original data from IllustrisTNG-100.Dotted blue lines represent the interpolation and smoothing with a 1 Myr time step.Second and fourth panels: evolution of the TNG-TVP rotation curves to the initial (at present -red colour), average (5 Gyr -green) and final moments of time (10 Gyr -blue).From left to right from top to bottom: #411321, #441327, #451323, #462077, and #474170.

Fig. B. 1 .
Fig. B.1.Examples of orbits of different 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, where the blue dot is the initial location and the red dot is the end location.
at present.