EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 673 (May 2023) A&A, 673 (2023) A152 Full HTML
Open Access
Issue
A&A
Volume 673, May 2023
Article Number A152
Number of page(s) 18
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202245117
Published online 24 May 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 ≳105 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 angular momentum (Lz) 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 − LzGrand 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⊙), 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⊙ 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.

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

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⊙ = 8.178 kpc (Gravity Collaboration 2019; Reid & Brunthaler 2004), a height above the Galactic plane Z⊙ = 20.8 pc (Bennett & Bovy 2019), and the velocity of the local standard of rest (LSR) VLSR = 234.737 km s−1 (Bovy et al. 2012; Drimmel & Poggio 2018). Accordingly, the Sun is located at X⊙ = −8178 pc, Y⊙ = 0 pc and Z⊙ = 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⊙ = 11.1 km s−1, V⊙ = 12.24 km s−1, and W⊙ = 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 Mpc3. In a box of this size, a sufficient number of the MW-mass disk galaxies can be resolved (in our sample, 54 objects1) with a mass resolution of 7.5 × 106 M⊙ for dark matter and 1.4 × 106 M⊙ for the baryonic particles. For our analysis, we identified the MW-like galaxy candidates in the Illustris simulations with at least 105 dark matter particles and at least 103 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).

Table 1.

Parameters of the time-varying potentials selected from the IllustrisTNG-100 simulation at redshift zero.

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−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). 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,

Φ tot = Φ d ( R , z ) + Φ h ( R , z ) = − G M d R 2 + ( a d + z 2 + b d 2 ) 2 − G M h · ln ( 1 + R 2 + z 2 b h ) R 2 + z 2 , $$ \begin{aligned} \begin{aligned} \Phi _{\rm tot}&= \Phi _{\rm d} (R,z) + \Phi _{\rm h} (R,z) \\&= - \frac{GM_{\rm d}}{\sqrt{R^{2}+\Bigl (a_{\rm d}+\sqrt{z^{2}+b^{2}_{\rm d}}\Bigr )^{2}}} - \frac{GM_{\rm h}\cdot \mathrm{ln}\Bigr (1+\frac{\sqrt{R^{2}+z^{2}}}{b_{\rm h}}\Bigl )}{\sqrt{R^{2}+z^{2}}}, \end{aligned} \end{aligned} $$(1)

where R = x 2 + y 2 $ R=\sqrt{x^{2}+\mathit{y}^{2}} $ is the planar galactocentric radius, z is the distance above the plane of the disk, G is the gravitational constant, ad is the disk scale length, bd, h are the disk and halo scale heights, respectively, and Md and M h = 4 π ρ 0 b h 3 $ M_{\mathrm{h}}=4\pi\rho_{0}b^{3}_{\mathrm{h}} $ (ρ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 functions 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 GitHub2.

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 identification, 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.

thumbnail 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 Mh, disk mass Md, NFW halo scale parameter bh, and the MN disk scale parameters ad and bd 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.

2.3. Integration procedure

For the GC orbit integration, we used the high-order parallel dynamical N-body code φ-GPU3 (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 η = 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 (ΔE/Et = 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 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 investigated 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⊙ (green circles) from Baumgardt & Vasiliev (2021)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.

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

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 velocity) 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 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 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.

thumbnail 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 in411321 TNG-TVP.

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 potential, 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 revelations). 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.

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

thumbnail 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 ≤ 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.

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 ≲ 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. 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 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 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%).

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

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 sample, Liller and Terzan 4 are the closest GCs to the Galactic centre.

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

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) onthe 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 website5 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.

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

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 online5. 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 trapezoial 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.

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.

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

thumbnail 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).

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 Δ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 the 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.

thumbnail Fig. 11.

Relative energy Δ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.

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.

  • ΔE/E = 45%…50%: Crater, Pal 3, Sagittarius II (black in Fig. 11) Eridanus, AM 1, and Pal 4;

  • ΔE/E = 42%…44%: NGC 2419, Pyxis, Pal 14, Lae 3 (orange in Fig. 11), and Whiting 1;

  • Δ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 analyse the phase-space distribution of the GCs at present in three phase-space coordinate combinations: total angular momentum (Ltot) versus total energy (E), z-component of the angular momentum (Lz) versus total energy (E), and the perpendicular component of the angular momentum (Lperp) versus total energy (E), where L tot = L x 2 + L y 2 + L z 2 $ L_{\mathrm{tot}}=\sqrt{L_{x}^{2}+L_{\mathit{y}}^{2}+L_{z}^{2}} $ and L perp = L x 2 + L y 2 $ L_{\mathrm{perp}}=\sqrt{L_{x}^{2}+L_{\mathit{y}}^{2}} $. The orientation of the Z-axis in our Cartesian galactocentric coordinates is directed towards 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. They orbit close to the Galactic centre because their total energies are highest (by absolute value). Most of the GCs with TB obits show from zero to some prograde rotation with the smallest perpendicular angular momentum component. This might indicate their common dynamical origin or possible orbital evolution as a consequence of their similar angular momentum loss timescales.

thumbnail Fig. 12.

Distribution of orbital types at redshift zero for 411321 TNG-TVP. From left to right: total angular momentum vs. energy (Ltot, E), zth component of the angular momentum vs. total energy (Lz, E), and perpendicular component of the angular momentum vs. total energy (Lperp, 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 (Ltot, E) and (Lperp, E) phase space, and Sagittarius II is not shown in the (Lz, E) phase space either because their values are extremely high.

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 Lperp with a weak trend for 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 high values of Lperp 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.

thumbnail 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, 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 (Ltot, E) and (Lperp, E) phase spaces and Sagittarius II is not shown in (Lz, E) phase space because their values are extremely high.

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.

thumbnail Fig. 14.

Evolution of the total angular momentum vs energy phase space (Ltot, E) (left panels, top to bottom), zth component of the angular momentum vs energy (Lz, E) (middle, top to bottom), and the perpendicular component of the angular momentum vs energy (Lperp, 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 (Ltot, E) and (Lperp, E) phase spaces, and Sagittarius II is not shown in (Lz, E) phase space because their values are extremely high.

thumbnail Fig. 15.

Same as in Fig. 14, but for 462077 and 474170 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 we 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 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 411321, 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 ≤ 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 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 Δ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 ex situ.

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.

1

MW-like potentials are presented on the web page of the project https://bit.ly/3b0lafw

2

The ORIENT: https://github.com/Mohammad-Mardini/The-ORIENT

3

N-body code φ-GPU: https://github.com/berczik/phi-GPU-mole

4

Error values for PMRA, PMDEC, RV, and D⊙ from https://people.smp.uq.edu.au/HolgerBaumgardt/globular/orbitstable.txt

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.

Acknowledgments

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

  1. Agertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64 [NASA ADS] [CrossRef] [Google Scholar]
  2. Armstrong, B. M., Bekki, K., & Ludlow, A. D. 2021, MNRAS, 500, 2937 [Google Scholar]
  3. Bajkova, A. T., & Bobylev, V. V. 2021, Res. Astron. Astrophys., 21, 173 [CrossRef] [Google Scholar]
  4. Bajkova, A. T., Carraro, G., Korchagin, V. I., Budanova, N. O., & Bobylev, V. V. 2020, ApJ, 895, 69 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bajkova, A. T., Smirnov, A. A., & Bobylev, V. V. 2021, Astron. Lett., 47, 454 [CrossRef] [Google Scholar]
  6. Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957 [NASA ADS] [CrossRef] [Google Scholar]
  7. Baumgardt, H., Hilker, M., Sollima, A., & Bellini, A. 2019, MNRAS, 482, 5138 [Google Scholar]
  8. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [Google Scholar]
  9. Bennett, M., Bovy, J., & Hunt, J. A. S. 2022, ApJ, 927, 131 [NASA ADS] [CrossRef] [Google Scholar]
  10. Berczik, P., Nitadori, K., Zhong, S., et al. 2011, in International Conference on High Performance Computing HPC-UA 2011, 8 [Google Scholar]
  11. Berczik, P., Spurzem, R., & Wang, L. 2013, in Third International Conference on High Performance Computing, HPC-UA 2013, 52 [Google Scholar]
  12. Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349 [Google Scholar]
  13. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  14. Boldrini, P., & Bovy, J. 2022, MNRAS, 516, 4560 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bovy, J. 2011, Ph.D. Thesis, New York University, USA [Google Scholar]
  16. Bovy, J., Allende Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131 [NASA ADS] [CrossRef] [Google Scholar]
  17. Carpintero, D. D., & Aguilar, L. A. 1998, MNRAS, 298, 1 [NASA ADS] [CrossRef] [Google Scholar]
  18. Chemerynska, I. V., Ishchenko, M. V., Sobolenko, M. O., Khoperskov, S. A., & Berczik, P. P. 2022, Adv. Astron. Space Phys., 12, 18 [CrossRef] [Google Scholar]
  19. Conroy, C., Naidu, R. P., Garavito-Camargo, N., et al. 2021, Nature, 592, 534 [NASA ADS] [CrossRef] [Google Scholar]
  20. Davies, J. J., Crain, R. A., Oppenheimer, B. D., & Schaye, J. 2020, MNRAS, 491, 4462 [NASA ADS] [CrossRef] [Google Scholar]
  21. Donnari, M., Pillepich, A., Joshi, G. D., et al. 2021, MNRAS, 500, 4004 [Google Scholar]
  22. Drimmel, R., & Poggio, E. 2018, Res. Notes Am. Astron. Soc., 2, 210 [Google Scholar]
  23. Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716 [CrossRef] [Google Scholar]
  24. Fattahi, A., Belokurov, V., Deason, A. J., et al. 2019, MNRAS, 484, 4471 [NASA ADS] [CrossRef] [Google Scholar]
  25. Fernández-Trincado, J. G., Minniti, D., Souza, S. O., et al. 2021, ApJ, 908, L42 [Google Scholar]
  26. Ferrone, S., Di Matteo, P., Mastrobuono-Battisti, A., et al. 2023, A&A, 673, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Garro, E. R., Minniti, D., Gómez, M., et al. 2022, A&A, 658, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Garrow, T., Webb, J. J., & Bovy, J. 2020, MNRAS, 499, 804 [NASA ADS] [CrossRef] [Google Scholar]
  31. Genel, S., Nelson, D., Pillepich, A., et al. 2018, MNRAS, 474, 3976 [Google Scholar]
  32. Gómez, F. A., Besla, G., Carpintero, D. D., et al. 2015, ApJ, 802, 128 [CrossRef] [Google Scholar]
  33. Grand, R. J. J., Deason, A. J., White, S. D. M., et al. 2019, MNRAS, 487, L72 [Google Scholar]
  34. Gratton, R., Bragaglia, A., Carretta, E., et al. 2019, A&A Rv., 27, 8 [NASA ADS] [Google Scholar]
  35. Gravity Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Habouzit, M., Genel, S., Somerville, R. S., et al. 2019, MNRAS, 484, 4413 [CrossRef] [Google Scholar]
  37. Haghi, H., Zonoozi, A. H., & Taghavi, S. 2015, MNRAS, 450, 2812 [NASA ADS] [CrossRef] [Google Scholar]
  38. Harris, W. E. 2010, ArXiv e-prints [arXiv:1012.3224] [Google Scholar]
  39. Harris, W. E., Harris, G. L. H., & Alessi, M. 2013, ApJ, 772, 82 [Google Scholar]
  40. Ibata, R., Malhan, K., Martin, N., et al. 2021, ApJ, 914, 123 [NASA ADS] [CrossRef] [Google Scholar]
  41. Ishchenko, M. V., Sobolenko, M. O., Kalambay, M. T., Shukirgaliyev, B. T., & Berczik, P. P. 2021, Rep. Nat. Academy Sci. Republic Kazakhstan, 6, 94 [Google Scholar]
  42. Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864 [Google Scholar]
  43. Katz, N., & Gunn, J. E. 1991, ApJ, 377, 365 [NASA ADS] [CrossRef] [Google Scholar]
  44. Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 [Google Scholar]
  45. Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R. D. 2013, A&A, 558, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202344233 [Google Scholar]
  47. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286 [Google Scholar]
  48. Laporte, C. F. P., Belokurov, V., Koposov, S. E., Smith, M. C., & Hill, V. 2020, MNRAS, 492, L61 [Google Scholar]
  49. Łokas, E. L. 2020, A&A, 638, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Lovell, M. R., Pillepich, A., Genel, S., et al. 2018, MNRAS, 481, 1950 [Google Scholar]
  51. Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 [NASA ADS] [Google Scholar]
  52. Malhan, K., Ibata, R. A., Sharma, S., et al. 2022, ApJ, 926, 107 [NASA ADS] [CrossRef] [Google Scholar]
  53. Mardini, M. K., Placco, V. M., Meiron, Y., et al. 2020, ApJ, 903, 88 [NASA ADS] [CrossRef] [Google Scholar]
  54. Marín-Franch, A., Aparicio, A., Piotto, G., et al. 2009, ApJ, 694, 1498 [Google Scholar]
  55. Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
  56. Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Mateu, C. 2023, MNRAS, 520, 5225 [Google Scholar]
  58. McConnachie, A. W. 2012, AJ, 144, 4 [Google Scholar]
  59. McConnachie, A. W., Higgs, C. R., Thomas, G. F., et al. 2021, MNRAS, 501, 2363 [NASA ADS] [CrossRef] [Google Scholar]
  60. Merritt, D. 1999, PASP, 111, 129 [NASA ADS] [CrossRef] [Google Scholar]
  61. Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
  62. Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235 [Google Scholar]
  63. Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
  64. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
  65. Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 [Google Scholar]
  66. Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2 [Google Scholar]
  67. Ness, M., & Lang, D. 2016, AJ, 152, 14 [NASA ADS] [CrossRef] [Google Scholar]
  68. Pagnini, G., Di Matteo, P., Khoperskov, S., et al. 2023, A&A, 673, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Panithanpaisal, N., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 920, 10 [NASA ADS] [CrossRef] [Google Scholar]
  70. Petersen, M. S., & Peñarrubia, J. 2021, Nat. Astron., 5, 251 [NASA ADS] [CrossRef] [Google Scholar]
  71. Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]
  72. Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196 [Google Scholar]
  73. Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872 [Google Scholar]
  74. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [Google Scholar]
  75. Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJ, 781, L31 [CrossRef] [Google Scholar]
  76. Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676 [Google Scholar]
  77. Sun, G., Wang, Y., Liu, C., et al. 2023, Res. Astron. Astrophys., 23, 015013 [CrossRef] [Google Scholar]
  78. Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484, 5587 [NASA ADS] [Google Scholar]
  79. Valcin, D., Bernal, J. L., Jimenez, R., Verde, L., & Wandelt, B. D. 2020, J. Cosmol. Astropart. Phys., 2020, 002 [CrossRef] [Google Scholar]
  80. VandenBerg, D. A., Brogaard, K., Leaman, R., & Casagrande, L. 2013, ApJ, 775, 134 [Google Scholar]
  81. van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35 [NASA ADS] [CrossRef] [Google Scholar]
  82. Vasiliev, E. 2019, MNRAS, 484, 2832 [Google Scholar]
  83. Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978 [NASA ADS] [CrossRef] [Google Scholar]
  84. Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874 [Google Scholar]
  85. Xiang, M., & Rix, H.-W. 2022, Nature, 603, 599 [NASA ADS] [CrossRef] [Google Scholar]
  86. Yun, K., Pillepich, A., Zinger, E., et al. 2019, MNRAS, 483, 1042 [Google Scholar]

Appendix A: Evolution of the circular velocity at the distance of the Sun for TNG-TVPs.

thumbnail 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−1; 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.

Appendix B: Orbit types

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

Appendix C: Globular clusters from the Gaia DR3 catalogue.

Table C.1.

Globular clusters from the Gaia DR3 catalogue.

All Tables

Table 1.

Parameters of the time-varying potentials selected from the IllustrisTNG-100 simulation at redshift zero.

In the text
Table C.1.

Globular clusters from the Gaia DR3 catalogue.

In the text

All Figures

thumbnail 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.
In the text
thumbnail 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 Mh, disk mass Md, NFW halo scale parameter bh, and the MN disk scale parameters ad and bd 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.
In the text
thumbnail 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.
In the text
thumbnail 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 in411321 TNG-TVP.
In the text
thumbnail 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.
In the text
thumbnail 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 ≤ 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.
In the text
thumbnail 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.
In the text
thumbnail 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.
In the text
thumbnail 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.
In the text
thumbnail 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).
In the text
thumbnail Fig. 11.

Relative energy Δ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.
In the text
thumbnail Fig. 12.

Distribution of orbital types at redshift zero for 411321 TNG-TVP. From left to right: total angular momentum vs. energy (Ltot, E), zth component of the angular momentum vs. total energy (Lz, E), and perpendicular component of the angular momentum vs. total energy (Lperp, 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 (Ltot, E) and (Lperp, E) phase space, and Sagittarius II is not shown in the (Lz, E) phase space either because their values are extremely high.
In the text
thumbnail 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, 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 (Ltot, E) and (Lperp, E) phase spaces and Sagittarius II is not shown in (Lz, E) phase space because their values are extremely high.
In the text
thumbnail Fig. 14.

Evolution of the total angular momentum vs energy phase space (Ltot, E) (left panels, top to bottom), zth component of the angular momentum vs energy (Lz, E) (middle, top to bottom), and the perpendicular component of the angular momentum vs energy (Lperp, 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 (Ltot, E) and (Lperp, E) phase spaces, and Sagittarius II is not shown in (Lz, E) phase space because their values are extremely high.
In the text
thumbnail Fig. 15.

Same as in Fig. 14, but for 462077 and 474170 TNG-TVPs.
In the text
thumbnail 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−1; 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.
In the text
thumbnail 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.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.