© The Authors 2025

1 Introduction

One of the rapidly developing fields of research in astronomy today is Galactic archaeology (Frebel 2010; Miglio et al. 2013; Spitoni et al. 2020; Koppelman et al. 2020; Xiang & Rix 2022; Deason & Belokurov 2024), which studies the history of the assembly of the Milky Way (MW). Based on observational data (Helmi 2020; Buder et al. 2021; Gaia Collaboration 2023), both astrometric and spectroscopic, as well as numerical simulations (Wang et al. 2024; Merrow et al. 2024). Galactic archaeology investigates the evolutionary processes that shaped our Galaxy from the earliest stages of the Universe’s existence, including a series of mergers between the MW and its neighbours.

Today we are able to identify at least five major merger events: Sagittarius, Gaia Enceladus (often referred to as Sausage, GE/S), Sequoia, Helmi, and Kraken, which are key to understanding the MW’s current structure and dynamics (Belokurov et al. 2006; Helmi et al. 2018; Myeong et al. 2019; Ibata et al. 2021; Kruijssen et al. 2020; Matsuno et al. 2022; Mateu 2023). Each of these events not only brought stars and dark matter but also contributed whole globular cluster (GC) systems. These GCs are today identified as an ex situ GC population in the MW. They are vital for tracing the global dynamical history of the Galaxy (Monkman et al. 2006; Lind et al. 2015; Escudero et al. 2018; Giusti et al. 2024). Their long lifetimes (Vandenberg et al. 1996; VandenBerg et al. 2013; Valcin et al. 2021) and stable internal structure (Ishchenko et al. 2024) make them excellent candidates for studying past interactions and large merger events, providing clues about the types of galaxies that contributed to the growth of MW and the processes involved in these mergers. Moreover, analysis of the chemical compositions and dynamics of these GCs can shed light on the different conditions and environments of their formation, enriching our understanding of galactic evolution. This framework of using GCs as a tracers highlights the intricate and often violent dynamical assembling history of our Galaxy, illustrating how it has evolved from a collection of smaller systems into the large spiral disc galaxy as we see today.

According to cosmological simulations, there is a correlation between massive merger events and possible highly radial orbits of satellite galaxies (Fattahi et al. 2019; Mackereth et al. 2019). The accretion process of large stellar systems, such as GE/S, into the host galaxy’s potential is a non-trivial task. This process leads to a complex distribution of energy and angular momentum in the debris of the accreted galaxy (Quinn 1984; Quinn & Goodman 1986; Tormen et al. 1998; van den Bosch et al. 1999), as well as to a complex chemical gradient (Kirby et al. 2011; Ho et al. 2015) in the Galaxy.

For instance, it is believed that the local stellar disc of the MW was primarily formed as a result of a merger with the massive GE/S satellite dwarf galaxy 8–11 Gyr ago or z ∼1–2 (Helmi et al. 2018; Belokurov et al. 2018; Pfeffer et al. 2020). This dwarf galaxy had a mass of ∼10^8−9.5 M_⊙, accounting for ∼10% of the MW’s mass at that time. This event was one of the most significant mergers in the MW and the last major one (Ruchti et al. 2015; Fragkoudi et al. 2020). Additionally, it has been suggested that the GE/S merger also might have played a significant role in the formation of the MW’s bar; see Kruijssen & Portegies Zwart (2009).

In the study of Koppelman et al. (2020) from the set of extensive N-body simulations we can conclude, that the most comparable results for the Gaia DR2 velocity distribution of the MW solar vicinity (<1 kpc) halo stars we obtain for the disc merger on the retrograde orbit with the orbital inclination angle of ∼30 degree. In addition, more than 20 GCs were added to the MW as a result of this merger (Massari et al. 2019; Kruijssen et al. 2020; Malhan et al. 2022; Chen & Gnedin 2024; Aguado-Agelet et al. 2025; Ceccarelli et al. 2025). These clusters are crucial for studying the Galaxy’s past interactions and dynamics. The merger also triggered intense bursts of star formation, contributing to the thickening of the MW’s disc.

As a first step in our current investigation, we carried out an extensive literature search for possible GE/S related GC systems. Here we mainly focused on the studies of Ceccarelli et al. (2025), Chen & Gnedin (2024), Malhan et al. (2022), Kruijssen et al. (2020), and Massari et al. (2019), which cover all the main aspects of observational and theoretical backgrounds of the MW GCs phase-space distribution analyses.

Based on the Gaia DR3 data for this set of GCs, we ran a 9 Gyr lookback time orbital simulation in the time-variable, cosmologically motivated, IllustrisTNG selected halo potential (Ishchenko et al. 2023a). Using this kinematical model for GCs orbital evolution and applying the energy and angular momentum phase-space limits for the GE/S objects from the above literature, we created our own list of the most probable GCs with a strong GE/S progenitor connection.

As a final step, we modelled, using the high-order dynamical N-body code, our set of probable GE/S GCs forward dynamical evolution (from −9 Gyr) up until today. Based on the detail orbital integration and stellar density distribution study of the escaped stars from these GCs, we finally redefined the current GE/S stream phase-space boundaries using our external time-variable potential (TVP). We also compared our GC’s escaped stars distribution with the catalogues of the very metal-poor (VMP) stars, observed today in the GE/S debris. Based on our calculations of the external galaxy model, we redefined the boundaries of the GE/S stream itself.

The paper is organised as follows. In Sect. 2 we analyse GCs, collected from the literature which may belong to the GE/S dwarf galaxy. In Sect. 3 we reconstruct the orbits of selected GCs up to a lookback time of 9 Gyr in a time-variable external galactic potential. In Sect. 4, we prepare a sample of GCs associated with the GE/S as a progenitor, based on the evolution of the orbital parameters. In Sects. 5 and 6, we carry out an N-body dynamical simulation of the selected GCs and compare the distribution of their escaped stars with the observed stars in the GE/S debris in energy and angular momentum phase-space.

2 Data sample

In this section, we give an overview of the literature and analyse the different authors’ classification prescription of GCs associated with the GE/S dwarf galaxy. The Gaia space telescopes have contributed significantly to this process. Thanks to the last catalogues of DR2 and DR3 Gaia Collaboration (2018); Gaia Collaboration (2020); Gaia Collaboration (2023) of Gaia, it has become possible to combine two methods of classification: kinematics and chemical. The present study employs a comparative approach that integrates the kinematic method with the analysis of chemical patterns. The results offer a comprehensive exploration of the interplay between GC orbital parameters and chemical patterns. In the following, we summarise the classification of the GC membership to the GE/S based on recent publications.

In earlier pioneering investigations, authors Mackey & Gilmore (2004) and Recio-Blanco (2018) suggested that ∼50 GCs may be formed ex situ and fall onto to our Galaxy’s gravitational potential during approximately seven to eight continuous mergers together with their dwarf galaxies. In Kruijssen et al. (2019), the authors make the suggestion that our Galaxy might have undergone at last three major mergers with different dwarf galaxies. The authors also make a sample of GCs that are associated with the remnants of the Canis Major, Kraken, and Sagittarius dwarf galaxies. This sample is based on summarising information about metal abundances from earlier works, and comparing it in the age-metallicity space. We need to note that most of the clusters from their sample of Canis Major GCs are identified in later works as systems that belong to the GE/S dwarf galaxy. In Myeong et al. (2018), the authors used integrals of motion to separate GCs into different groups. In a later publication (Myeong et al. 2019), the authors find 21 potential GE/S GCs based on their own dynamical analysis.

Recently Aguado-Agelet et al. (2025); Ceccarelli et al. (2025) determined the age and metallicity of 13 GCs that are associated with the GE/S dwarf galaxy, applying the data from the HST telescope and theoretical isochron fitting, based on the new

BaSTI

stellar evolution model mesh (Hidalgo et al. 2018). For the association of the selected GCs with the GE/S stream, the authors applied only their chemical analyses. The obtained age-metallicity relationship separates the clusters into two formation epochs, with a time interval of about 2 Gyr. As a final result, the authors present the list of GCs potentially associated with the GE/S stream: NGC 288, NGC 362, NGC 1261, NGC 1851, NGC 2298, NGC 2808, NGC 5286, NGC 5897, NGC 6205, NGC 6341, NGC 6779, NGC 7089, and NGC 7099.

In the study of Chen & Gnedin (2024), the authors proposed a method of classifying progenitors of GCs based on the analysis of kinematic characteristics (i.e. integrals of motion and action angles) together with their chemical composition and age. The authors applied their own developed clustering method using up to ten different parameters, including orbital parameters (actions J_r, J_φ, and J_θ, apocenter r_apo and pericenter distances r_peri, eccentricities,

ecc

, Galactocentric distance, D, and total energy, E) with their chemical properties (age and metallicity). It needs to be noted that for the kinematic analysis the authors applied the ‘static’ galactic potential based on the standard three-component model of Bovy (2015). The authors associated the next GCs with the GE/S dwarf: ESO-SC06, IC 1257, NGC 1261, NGC 1851, NGC 1904, NGC 2298, NGC 2808, NGC 288, NGC 362, NGC 4147, NGC 5286, NGC 6205, NGC 6229, NGC 6341, NGC 6779, NGC 6864, NGC 7089, NGC 7099, NGC 7492, Palomar 2, NGC 5634, NGC 5904, NGC 6981, NGC 7078, NGC 6426, and NGC 6584.

In Malhan et al. (2022), the authors proposed the classification of ∼170 GCs that combined their chemical and kinematic parameters. The authors used extensive sky surveys (LAMOST DR7 and APOGEE DR17) together with astrometric data collected from the Gaia EDR3 catalogue. The classification was carried out by calculating action angles together for GCs and progenitors (streams). They applied the static and axisymmetric model of the Galactic potential (McMillan 2017). The final list of GCs that Malhan et al. (2022) attribute to the GE/S dwarf galaxy is: NGC 7492, NGC 6229, NGC 6584, NGC 5634, IC 1257, NGC 1851, NGC 2298, NGC 4147, NGC 1261, NGC 6981, NGC 1904, NGC 7089, NGC 5904, and NGC 6864.

In earlier papers by Massari et al. (2019) and Kruijssen et al. (2020) on the determination of the GCs membership with the GE/S, the authors mainly used astrometric data that was taken from the earlier Gaia DR2. During their investigations, they applied the same static axisymmetric Galactic potential (McMillan 2017) to compare the orbital parameters of the GCs and the GE/S stream. The list of GCs that Massari et al. (2019) attribute to the GE/S dwarf galaxy is: Djorg 1, NGC 4833, NGC 5897, NGC 6235, NGC 6284, Terzan 10, NGC 5139, ESO-SC06, IC 1257, NGC 1261, NGC 1851, NGC 1904, NGC 2298, NGC 2808, NGC 288, NGC 362, NGC 4147, NGC 5286, NGC 6205, NGC 6229, NGC 6341, NGC 6779, NGC 6864, NGC 7089, NGC 7099, NGC 7492, Palomar 2, NGC 5634, NGC 5904, Palomar 15, NGC 6101, and NGC 3201.

Contrary to the investigation above, in the work of Kruijssen et al. (2020) only the

ecc

and r_apo parameters as we see them today were studied. Here the authors also applied the simple age-metallicity estimation for GCs. The resulting list of GCs is attributed to the GE/S dwarf galaxy: NGC 288, NGC 362, NGC 1261, NGC 1851, NGC 1904, NGC 2298, NGC 2808, NGC 4147, NGC 4833, NGC 5286, NGC 5897, NGC 6205, NGC 6235, NGC 6284, NGC 6341, NGC 6779, NGC 6864, NGC 7089, NGC 7099, NGC 7492, and NGC 5139.

In Table 1 we present the combined sample of GCs associated with the GE/S, based on the described above publications. Our table consists of several pre-sets of GCs. For the first pre-set, we assumed GCs that we see in all five publications (4 GCs). For the second pre-set, we combined the GCs that appeared in at least four papers (12 GCs). The other two pre-sets have objects from the lists of only three, two, or even one publication. Based on this, we assume that these sets possibly contain the less probable (or less confident) GCs from the GE/S stream members.

3 Collected GCs’ orbital reconstruction as a point mass in time-variating potential

In this section, we carry out the orbital reconstruction of 36 GCs from the Table 1 inside the Galactic time-variable external potential. This gives us the opportunity to investigate the orbital evolution inside the MW of our selected GCs on the cosmological timescale.

For orbital reconstruction and numerical integration, we used the high-order parallel N-body code φ-GPU, which is based on fourth-order Hermite integration together with the hierarchical individual block timestep scheme Berczik et al. (2011, 2013). This code was already successfully applied to the MW GCs global dynamical evolution in our previous works Ishchenko et al. (2023a,c,b). To reconstruct orbits and calculate the main orbital parameters for our pre-sets of the 36 GCs during 9 Gyr of lookback evolution, we integrated each of them as a single point mass. The current positions, proper motions, radial velocities, and heliocentric distances were selected from the catalogue of Baumgardt & Vasiliev (2021)¹ and presented in Table A.1.

To calculate the Cartesian Galactocentric co-ordinates, we used the Sun position X_⊙ = −8.178 kpc, Y_⊙ = 0 kpc, and Z_⊙ = 0.0208 kpc Gravity Collaboration (2019); Bennett & Bovy (2019). For the local standard of rest (LSR) velocity, we used 234.737 km s⁻¹ Reid & Brunthaler (2004), and for the peculiar velocity of the Sun with respect to the LSR, we used U_⊙ = 11.1 km s⁻¹, V_⊙ = 12.24 km s⁻¹, and W_⊙ = 7.25 km s⁻¹ (Schönrich et al. 2010). In detail the co-ordinate and velocity transformations are described in Johnson & Soderblom (1987). For the equatorial position of the north galactic pole (NGP), we used the updated values from Karim & Mamajek (2017): RA_NGP = 192°.7278, Dec_NGP = 26°.8630, and θ_NGP = 122°.9280.

For our purpose we applied the time-variable external MW potential Ishchenko et al. (2023a), based on the IllustrisTNG-100 cosmological simulation database Nelson et al. (2019b,a). From our pre-selected list of the MW-like potentials, we used the

411321

TVP presented in Fig. 1. This potential at present has global parameters very similar to our Galaxy today (for more details, see Sect. 2.2 in Ishchenko et al. 2023a). The potential consists of two components: a dark matter halo, described by the potential, Φ_halo, from Eq. (1) of Navarro et al. (1997), and a disc, described by the potential, Φ_disc, from Eq. (1) of Miyamoto & Nagai (1975): $$\begin{array}{r}\Phi ={\Phi }_{halo}+{\Phi }_{disc}\\ =-\frac{G{M}_{halo}\ln \left(1+\frac{\sqrt{R^2+z^2}}{b_{halo}}\right)}{\sqrt{R^2+z^2}}-\frac{G{M}_{disc}}{\sqrt{R^2+{\left(a_{disc}+\sqrt{z^2+b_{disc}^2}\right)}^2}}.\end{array}$$(1)

In the above equation, $R=\sqrt{x^2+y^2}$ is the radial distance from the Galactic center, M_halo and M_disc are the mass of the MW halo and the mass of the MW disc, respectively, and b_halo, a_disc, and b_disc are the spatial scales of the halo and disc potentials. Figure 1 shows the time evolution of the main TVP potential

411321

parameters M_halo, M_disc, b_halo, a_disc, and b_disc during the whole 12 Gyr time evolution of the MW system. Recently, our TNG TVP potentials (particularly the

411321

model) were also included in the very popular N-body Galactic dynamical code

PeTar

(Wang et al. 2020). One can find the recent implementation of the updated

PeTar

code in a project GitHub link².

As a illustration of our GCs typical orbits, we present in Fig. 2, the 9 Gyr backward orbital evolution of the selected clusters from our list (see more detail explanation in the Sect. 4.4). As we see for example orbits in our TNG TVP potential, they have a significant long term orbital evolution including the strong precession of the orbits along the Galactic Z axis.

Also, in Table A.1 we show the relative errors for each GCs in proper motions, radial velocities, and heliocentric distances based on Gaia DR3 measurements (Baumgardt & Vasiliev 2021). Most of the GCs have high-accuracy measurements in the 6D kinematic parameters, with relative errors below ∼2%. In some specific cases, due to, the low proper motion of NGC 6584, the error in µ_α above ∼10%. For the NGC 7089 cluster radial velocity error, we go up to ∼5%. The detailed analyses about errors are represented in Appendix A. In Appendix A we also discuss the error values of 6D components for each GC and their possible influence on the long-term orbital evolution.

4 Procedure of identification GCs as GE/S membership in TVP

Here, we describe the (i) concept of identifying GCs as possible members of the GE/S stream based on (ii) energy–angular momentum evolution and orbital parameters in a TVP. In subsequent sections, we present (iii) our final sample of GCs connected with the GE/S, classified in this way together with the discussion (iv) about uncertainty surrounding the origin of some of these GCs.

4.1 Concept of identification

Table 2 shows the possible limits in the phase-space parameters of GCs with GE/S connection, which we collected from recent publications (Sun et al. 2023; Malhan et al. 2022; Massari et al. 2019; Limberg et al. 2022). In works Sun et al. (2023), Massari et al. (2019), and Malhan et al. (2022), the authors mainly use the current distribution in phase-space of kinematic parameters and integrals of motion for the selected GCs. Based on these criteria, they define the possible ranges for the association of concrete GCs with the GE/S stream. In article Limberg et al. (2022) the authors set the limits by distribution of possible individual GE/S stellar sample that have been selected by their chemical abundances. The authors made a cleaner GE/S sample by discarding stars with chemical abundances compatible with the in situ regions of the [Al/Fe] versus [Mg/Mn] space. In contrast to our current investigation, all of these previous studies of course only defined the current dynamical parameters, i.e. they did not study the possible changes in GCs and GE/S stream orbital motion in the past.

As we see from Table 2, the limits (or boundaries) for L_z, total energy, and

ecc

are defined in most of the studies. The boundaries for L_perp and r_apo are defined in only a few of these studies. In Figs. 3 and 4, we highlighted these parameter boundaries as different colour areas. One can note that the total energy range is quite different between these publications. The Malhan et al. (2022) study especially has a quite narrow energy range, as does the GCs from the GE/S stream members.

In the previous section, we describe our N-body modelling with the aim of orbital motion reconstruction of 36 GCs inside

411321

TVP up to 9 Gyr lookback time. We show the reconstructed orbital evolution for each GC (colour thin lines), grouped according to the pre-sets from Table 1.

We assume that GCs can have a GE/S dwarf galaxy origin if more than 80% of their trajectories during the full 9 Gyr evolution fall within the selected coloured areas; see Figs. 3 and 4. These areas represent the parameters for L_z and L_perp, together with the total energy as the three main identification criteria. The second criterion is the orbital parameters a,

ecc

, r_apo, and r_peri. In cases of disputed identification, whether or not a GC belongs to GE/S, chemical patterns will play a decisive role.

4.2 GC orbital parameters matching GE/S limits

In Fig. 3 we show the evolution of the total specific energy and angular momentum components L_z and L_perp for individual GCs, obtained in TVP. As can be seen in upper row for 1–3 pre-sets, all of the cclusters mainly fall within the limits, especially in L_z (Limberg et al. 2022; Malhan et al. 2022 and Sun et al. 2023) boundaries. All of the clusters also mainly fall within the limits of the regions by energy from −17 to −8×10⁴ km² s⁻². In pre-set 4 the next five clusters do not correspond to the limits: NGC 6235, NGC 6426, and NGC 7078, and especially NGC 3201 and NGC 6101, which are completely out of the L_z limits ≈ −3 × 10³ kpc km s⁻¹.

Analysing the bottom row in Fig. 3, it can be seen that all of the GCs in the first pre-set fall within the limits. In contrast to the evolution of the GCs in L_z, we can see clusters that do not correspond to the specific energy limits: NGC 7099 (black line) from the second pre-set, and NGC 5139, NGC 6205, and NGC 5897 (green, light blue, and beige) from the third pre-set. The GCs from the fourth pre-set – NGC 6284, Terzan 10, NGC 4833, and NGC 3201, which falls within the boundaries by L_perp – are also not compatible with the specific energy limits. Most of them have, for almost all the evolutionary time, a larger bounding specific energy from −22 to −18.5×10⁴ km² s⁻².

Evolution of the orbital parameters during a 9 Gyr lookback time integration that we present in Fig. 4. The coloured areas represent the distribution limits of the GC’s GE/S for r_peri, r_apo (upper row), together with a and

ecc

illustrated in the lower row. Note that the limits of the orbital elements are only defined only in two from our main four papers: see Malhan et al. (2022) and Sun et al. (2023).

For the pre-set 1, one can see a good correspondence of all the orbital parameters of individual GCs with the predefined ranges from the publications from Table 2. As can be seen for pre-sets 2–4, the clusters with

ecc

< 0.7 completely out of the ranges. These are mainly: NGC 7099, NGC 6235, NGC 7078, NGC 6101, NGC 5139, NGC 5897, NGC 6205, NGC 288, and Terzan 10. The GC’s distribution in the limits of r_peri and r_apo also show a large mismatch for almost all of the GCs in the second to fourth pre-sets. The only exceptions are a short list of GCs from the second pre-set that are almost entirely inside the r_peri and r_apo ranges: NGC 4147, 5634, 5904, 6864, and 7492.

4.3 Resulting GE/S GCs sample

Based on our numerical integration in the TVP external potential of the GCs, which we present in Figs. 3 and 4 with a corresponding analysis, we can summarised our pre-sets in the tree category. Namely, the first category represents the ‘most probable’. Here we represent the GCs that have a strong association with six parameters that determine the limits of GE/S: L_z, L_perp, the total energy, r_apo, r_peri, and

ecc

. We get 15 GCs that fall into that category: NGC 1261, NGC 1851, NGC 7089, NGC 2808, NGC 2298, NGC 6864, NGC 6779, NGC 4147, NGC 5634, NGC 5904, NGC 7492, NGC 1904, IC 1257, NGC 6229, and NGC 6981. All of these clusters have a high conformity in orbital parameters and energy versus angular momentums during all 9 Gyr of integration of their orbital evolution. NGC 2808, NGC 1904, and IC 1251 do not satisfy the r_peri limit according to Malhan et al. (2022), but in our view it is not a critical parameter. This is first of all due to the quite large uncertainties in their definition and also the fact that this parameter was derived for our GC samples in only one publication Malhan et al. (2022). The clusters we selected are also well suited in terms of chemical parameters. They are quite compatible with the GE/S general chemical pattern. Here only NGC 6229 and NGC 6981 can be under some discussion; see more details about the membership of these GCs to our ‘most probable’ group in Sect. 4.4.

We assume the second category to be ‘tentative’. Here, we include clusters that do not demonstrate full consistency in all six parameters, but that could correspond to the GE/S progenitor based on the total energy and angular momentum components, L_z, L_perp, or that have a significantly compatible chemical pattern with the GE/S system. In the end, we assume the next nine clusters to be ‘tentative’: NGC 362, NGC 5286, NGC 6341, NGC 5897, NGC 7099, NGC 6426, ESO SC06, Djorg 1, and NGC 6426. The main reason for classifying the cluster in this group is its low energy value: < −18×10⁴ km² s⁻² (see Fig. 3, bottom row).

Finally, in the ‘excluded’ category we collected the clusters that in most of the parameters do not correspond to the specified limits or for which at least one of the parameters has a large discrepancy with the limits. In Table 3 we present the resulting sample of GE/S clusters. In total, of the 36 GCs that were selected from the recent literature, we find 24 GCs that are the most probable and tentative candidates and that we are highly confident have a GE/S dwarf galaxy origin. Also in Table B.1 we present detailed statistics on the orbital components and energy versus angular momentum values in phase-space in the context of the number of publications in which a particular cluster was mentioned.

In our current investigation, we use one of the five available MW-like TVP potentials selected from the IlustrisTNG-100 database; see Fig. 2 and Table 1 in Ishchenko et al. (2023a). To make our results in Table 3 more robust, we selected an additional external potential –

441327

– and examined the impacts of these changes on the evolution of GCs’ orbital parameters and energy versus momentum values. We carried out the same procedure as is described in Sect. 3 for our pre-selected GCs samples, collected in Table 1. In Figs. D.1 and D.2, we present the evolution of the GCs’ phase-space parameters during the 9 Gyr lookback time integration in

441327

TVP external potential for our created pre-sets of GCs.

We also provide a comparison of stellar mass loss in two simulations using both the time-dependent (TVP) and time-independent (TVP-FIX) potentials. We found slightly different mass loss, but during the whole simulation the difference between the TVP and TVP-FIX results is not more than a few percentage points; see the grey curves in Fig. C.1. Also for more details see Appendix C.

Our additional tests incorporating a simplified dynamical friction model reveal that its influence on orbital evolution is modest over 9 Gyr, even for clusters with low pericentric distances. As is shown in Fig. D.4. The inclusion of dynamical friction causes only minor deviations in the galactocentric distance and orbital shape, without significantly altering energy or angular momentum; for more details see Appendix D.3.

By analysing the behaviour of orbital parameters, it is possible to identify clusters that demonstrate an even more striking discrepancy with the boundaries of the GE/S; namely, NGC 5139, NGC 6205, NGC 5897 from third pre-set and NGC 4833, NGC 6284, Terzan 10, NGC 6584, and NGC 7078, which have an even stronger bound connection with the Galaxy core than in the

411321

TVP. However, aside from the stronger bounding energies, the orbital and angular momentum parameters do not change significantly in the

441327

potential relative to the

411321

potential. So, the distribution of the orbital parameters and angular-momentum variables of the ‘most probable’ and ‘tentative’ clusters satisfies all the limits of Table 2.

We also present summaries of the 3D plots in Fig. 5, with the GCs distribution in the context of phase-space parameters: L_z, L_perp, and the total energy and orbital parameters, r_apo, r_peri, and

ecc

. For each GC we represent the time evolution during the integration up to 9 Gyr. We clearly see the separation between our tree categories, especially in orbital parameters. The GCs from the ‘most probable’ and ‘tentative’ categories have, as was expected, high values in the

ecc

: 0.7–0.99, while the ‘excluded’ category demonstrate the large variation in

ecc

from 0.3 to 0.9. Also GCs from this category appear to have a low total energy: from −22 to −18×10⁴ km² s⁻². The other two categories range from −18 to −12×10⁴ km² s⁻² during their entire evolution, which is closer to the GE/S stream distribution in the Galaxy.

4.4 GCs under discussion

Here, we would like to discuss in more detail our resulting lists in Table 3. The clusters NGC 6229 and NGC 6981, which we assume to be the ‘most probable’, stand out from the other GCs in the group due to the lack of their chemical abundance observations; see Table B.1. Note that the authors of references Massari et al. (2019) and Malhan et al. (2022) consider NGC 6229 and NGC 6981 to correspond to the GE/S progenitor, whereas a later publication (Callingham et al. 2022) associates these GCs with the group of the Helmi stream. On the other hand, they demonstrate a strong association with GE/S in all kinematic parameters also based on our own numerical simulations in TVP.

In the tentative category we want to note a group of clusters – NGC 362, NGC 5286, NGC 5897, NGC 6341, and NGC 7099 – that have a strong chemistry association with GE/S, but kinematic parameters that have lower energy values (from −20 to −16×10⁴ km² s⁻²) compared with others GCs here. Almost all of the objects from this group have quite similarly shaped orbits. For example, NGC 5286 and NGC 5897 illustrate well the orbits from this group (see, Fig. 2). For these GCs we have a typically strong and multiple passages in the Galactic thick disc. During the GE/S accretion event, it is possible that some clusters, before they become a part of our Galaxy, were stripped by some dwarf galaxies (which themselves earlier were accreted by the MW) and as a result they can form a subgroup of GCs inside the current GE/S stream. These objects, as a result of such a subsequent mergers, can have stronger bounding energies, smaller apocentre distances, and also lower eccentricities.

Note that the authors in Malhan et al. (2022) classified this group (the five systems mentioned above) as members of the Pontus dwarf galaxy. However, the subsequent analysis in Ceccarelli et al. (2025) shows that this group is indeed strongly associated with GE/S by age and also metal abundances.

GC Djorgovsky 1 has much in common with GE/S clusters in kinematic parameters. However, we find only one publication that has studied the chemical abundances in this object; see Table B.1. So we consider to assume this cluster to be in the tentative group of our sample.

In the excluded group, we also have a few objects with some specific characteristics. With regard to the clusters NGC 288 and NGC 6205 shown in Aguado-Agelet et al. (2025) [Fig. 3], the authors found that these two clusters were much older than the GE/S dwarf galaxy. At the same time, the kinematic coincidence of the parameters of the clusters with GE/S can be explained by the fact that during the falling of the GE/S galaxy, the orbits of these two clusters were significantly perturbed (Ceccarelli et al. 2025). Based on these speculations, we excluded them from being probable GE/S stream members.

We also need to explain, in more detail, why in the excluded group we have Palomar 2 and Palomar 15. These clusters have angular momentum and energy values that are very typical of GE/S clusters. Also we have a quite limited knowledge of the chemical composition of Palomar 2 and 15 (similar to the case of Djorgovsky 1). We exclude these objects from being probable GE/S members, mainly based on their quite large semi-major distances (∼20–25 kpc) and high eccentricity values (∼0.95); see the bottom row of Fig. 4.

NGC 5139 (also known as Omega Centauri) is the most discussed cluster. Some of the authors (Massari et al. 2019; Kruijssen et al. 2020) consider this object to be a possible stellar core of the former GE/S dwarf galaxy or Sequoia dwarf based on the age-metallicity relation. In Fig. 2 we also present the orbit of this cluster. In a sense of binding energy, this object is one of the most bound systems in our sample (Fig. 3). At the same time, as can be seen from the Table B.1, this cluster has typical angular momentums compatible with GE/S only in L_Z. In the absence of other compelling factors, we consider leaving this cluster in the excluded group.

We also note that in the most probable sample we have ten GCs that have clockwise rotation and five that have anticlockwise rotation (in the Cartesian Galactocentric reference frame), while in the tentative sample six GCs demonstrate anticlockwise rotation. At the same time, the author states that the GE/S dwarf galaxy was falling towards our Galaxy in an orbit with anti-clockwise rotation Helmi et al. (2018); Koppelman et al. (2020).

5 N-body modelling of the GCs GE/S members

In this section we carry out numerical modelling of 24 GCs in order to estimate the distribution of stars in the phase-space that no longer belong to the GE/S clusters today.

5.1 GCs GE/S initial conditions and integration procedure

The aim of this N-body simulation is to study the distribution of stars in the Galaxy that have been lost by the GCs GE/Ss due to their own relaxation processes during the full 9 Gyr orbital evolution. Since our primary focus in the modelling was the dynamical behaviour of escaped, sub-solar-mass, long-lived stars within the cosmological timescale and a varying external Galactic gravitational potential (IllustrisTNG-100 motivated –

411321

TVP). In the current simulation, we turned off the individual particles’ stellar evolution during the N-body simulations. While a detailed investigation into GC’s internal stellar mass loss evolution would indeed require the incorporation of the stellar evolution module, our current study specifically targets stars with lifetimes of the order of 9 Gyr and more. These objects can travel, after escape from GCs, independently within the whole Galaxy. For this specific subset of stars, the time-dependent stellar evolution mass loss can be reasonably neglected. A similar approach was already successfully applied in our previous paper Berczik et al. (2024) in which we study the dynamical evolution of the escaped stars from the Reticulum II dwarf galaxy.

As a basic initial physical model for all the 24 GCs at 9 Gyr ago, we chose the equilibrium King model (King 1962) with the concentration parameter W₀ = 8.0 and with the fixed total mass of 10⁵ M_⊙. We chose these parameters to initially have a maximally compact King object in terms of the half-mass radius to tidal radius ratio. This allows us to have systems with the maximum lifetime and long survivability in the TVP external tidal field.

We varied the clusters’ remaining two physical parameters; namely, the number of particles and the initial half-mass radius, depending on the types of cluster orbits inside the Galaxy. Thus, for clusters whose orbits are relatively close to the centre of the Galaxy (r_apo less then 10 kpc) we used a half-mass radius of 2 pc and the initial number of particles N = 45–50k (depending on the clusters current masses).

For clusters that have long, elongated radial orbits and a small number of passages near the Galactic centre, we applied models with a larger half-mass radius of 4 pc and with the initial number of particles 40–45k. Using these two physical parameters (half-mass radius and initial number of particles), we tried to provide the more or a less comparable stellar escape rate of low-mass stars from the clusters with significantly different initial masses.

Besides low particle numbers in our N-body models (compared to the real number of stars in the Galactic GCs), we managed to provide the long-term (up to 9 Gyr) survivability of our dynamical models, which allow us to have a long term ‘source’ of escaped stars from these clusters. We used these stars in our study as ‘indicators’ (tracers) of the long-term orbital motion of our selected GCs. This approach – using a relatively low N for each GCs – enables us to run a significantly larger set of N-body simulations within a realistic computer running time frame.

To evaluate the possible influence of initial conditions on the dynamical evolution of GCs, we performed a number of additional simulations with different numbers of particles for several GCs from our list. We discuss the evolution of the star loss due to the orbital motion and relaxation time for these GCs in Appendix C.

5.2 Distribution of stars in phase-space at the present time from the GE/S GCs sample

In Fig. 6 we show the distribution of stars in the Galaxy at the present time, based on our sample of 24 GCs from the GE/S groups (most probable and tentative) during the 9 Gyr of their evolution. As we see, the global stellar distribution (escaped stars and stars that still belong to the cluster) in the L_tot space have quite a wide distribution range: from 0 to 1.8×10³ kpc km¹ s⁻¹ and from −18 to −12.2×10⁴ km² s⁻² in energy values. In the L_z component, as was expected, we see a quite narrow distribution of all stars for each GC. Most of these stars today are already not bound gravitationally to the progenitor GC itself. Precisely these unbound stars form our GE/S stellar debris in the phase-space region; see Fig. 6 (thick red box) and ‘in this paper’ in Table 4.

As we see from Fig. 6, the densest stellar regions overlap quite well with the possible limits of the GE/S stream in the energy–angular momentum phase-space, presented in the publications of Ye et al. (2024) and Malhan et al. (2022). As a side remark, we need to note that the limits used in the paper Liu et al. (2024) are obviously not compatible with our GCs GE/S stellar debris distribution.

6 Phase-space comparison of the GE/S GCs’ escaped stars with the metallicity-selected individual GE/S stars

Our comparison is based on the GE/S VMP star sample, which are investigated in detail in recent publications: Ye et al. (2024), Ernandes et al. (2024), and Molaro et al. (2020). All of these samples focus on their selection criteria of stars that are connected with the GE/S galaxy based on chemical abundances. We make a comparison between these GE/S stellar samples (selected by chemistry) and our resulting energy and angular momentum distribution of the stars from GCs of GE/S for today in our dynamical external TVP.

GE/S stars from Ye et al. (2024). This study (besides the chemical selection criteria) is based on the VMP stellar samples using kinematic and action-angle analysis of the individual stars’ orbits. The authors analysed the VMP stars sample based on the Gaia DR3 catalogue, together with the application of the radial velocities from the LAMOST DR9 catalogue. In summary, the authors found 1070 stars that are probably associated with the original GE/S progenitor. From this sample we obtained 343 stars with full 6D astrometry data, available in the Gaia DR3 catalogue. Another 727 GE/S VMP stars have low-accuracy (compared to Gaia DR3) radial velocities based on the LAMOST catalogue.

GE/S stars from Ernandes et al. (2024). Here, the GE/S VMP stellar sample based on the data from the SAGA base. The authors make a stellar sample based on the metallicity abundance range that is expected for GE/S stars: as −2.2 < [Fe/H] < −0.5. From this sample, we obtained 73 stars, which have Gaia DR3 astrometry data.

GE/S stars from Molaro et al. (2020). The authors investigated the stellar abundances

Li

and

Be

from the possible members of GE/S. In this work, to identify stars that are connected with GE/S progenitor, the authors additionally defined the limits in the angular momentum L_z component. Note that the sample is based on the data from Gaia DR2 catalogue together with the APOGEE sky survey. From this sample we collected 171 stars, with Gaia DR3 astrometry data.

The total distribution of the combined 587 GE/S stars presented in Fig. 7 in energy and angular momentum co-ordinates. To calculate these values for each star, we used their 6D kinematic data from Gaia DR3, and our external Galactic potential

411321

TVP. In Fig. 7 we plot these 587 GE/S stars together with the ejected stars distribution of GE/S GCs for today. The Ernandes et al. (2024) sample shows the greatest overlap with the phase-space distribution of the ejected stars from our two GE/S GCs categories. In the top row we put the 15 most probable GE/S GCs. In the middle row we put the nine tentative GE/S GCs. In the bottom row, we show the combination of these two GCs samples in one figure. The red line on the bottom panels represents the limits in the energy and angular momentum space, according to our modelled ejected stars from the GCs of the GE/S. The corresponding values of the possible GE/S debris boundaries we also present in the Table 4. As one can see from Fig. 7, all the GE/S stellar samples overlap quite well with the groups of escaped stars from the GE/S GCs.

7 Conclusions

High-precision astrometry data from Gaia in the last few years has allowed us to identify the most massive stellar debris that was supposed to be the remnant of the major merger with the dwarf galaxy named GE/S (Helmi et al. 2018). Further investigations allowed us to combine stellar streams from such events with GCs. Based on the recent literature review, we collected and combined catalogue of the GCs, which potentially can have the GE/S dwarf galaxy origin. We performed N-body backward integrations for 36 collected from literature GCs, up to the probable time of the merging event of the GE/S with our MW, i.e. up to 9 Gyr. We carried out the GC orbital reconstruction in the time-variable external MW-like galactic potentials, selected from the IlustrisTNG-100 database. Based on Figs. 3 and 4, we selected GE/S GCs using the limits of their distribution in orbital parameters and energy with angular momentum for GE/S GCs that we find in the literature; see Table 2.

Note that applying a time-varying potential to define the boundaries of GCs offers a more flexible and realistic framework compared to static potentials. Unlike static models, which constrain GCs to fixed energy and angular momentum values and thus impose stricter membership criteria, for the case with time-dependent potential we have more natural variations in orbital parameters and as a result more variability in the membership criteria too.

Therefore, we analysed the evolution of orbital parameters and angular momentum and created a sample of 24 GCs that most probably have an association with the GE/S progenitor. We distributed the collected GCs into two categories: most probable and tentative. As a result we have a total list of 24 individual GCs: NGC 1261, NGC 1851, NGC 2298, NGC 7089, NGC 1904, NGC 2808, NGC 4147, NGC 5634, NGC 5904, NGC 7492, NGC 6229, NGC 6864, IC 1257, NGC 6981, NGC 6779, NGC 362, NGC 5286, NGC 6341, NGC 6584, NGC 6426, NGC 7099, ESO SC06, Djorg 1, and NGC 5897; see Table 3. In the next step, we ran a set of dynamical N-body simulations to study the stellar ejection processes due to their dynamical evolution and orbital motion in the TVP external potential.

We identified a possible dynamical connection between individual stars with extremely low metallicities (VMP stellar samples) and stars, ejected from GCs that exhibit a significant connection with GE/S; see Fig. 7. This connection is explored using the current phase-space information for these stars, which provides detailed insights into their positions and velocities within the MW, applying the time-variation external Galactic potential. Furthermore, we applied the advanced dynamical star loss model to simulate the behaviour of GCs together with the ejected stars that were associated with the GE/S over a period of 9 Gyr of orbital motion. This model allows us to understand the evolution of the our 24 selected GCs, including tidal stripping and other gravitational influences, which affect the global stellar distribution over time.

We defined the GE/S debris limits based on the escaped stars from our selected 24 GE/S GCs that were dynamically modelled in external time-variable Galactic potential. Based on this modelling, we obtain the next phase-space GE/S limits: specific energy −18 to −12.2 ×10⁴ km² s⁻², in L_z from −0.98 to 0.72 × 10³ kpc km s⁻¹; for L_perp: from 0 to 1.8 ×10³ kpc km s⁻¹; see Fig. 7. Also we redefined the possible limits of the GE/S GCs in Galactic area with the following orbital parameters: an apocen-tre distance of 10–28 kpc, a pericentre distance of 1–4 kpc, and limits of <18 kpc in the Galactocentric radius and <|15| kpc in the Z direction; see Table 4.

Our current study shows us a possible future extension of our models with the inclusion of the more recent IllustrisTNG-50 cosmological database for the detailed simulation of proto-MW and GE/S mergers. To conduct a detailed study on the origin of GCs, it is also necessary for us to create a comprehensive model of the two galaxies (proto-MW and GE/S), including their GCs subsystems too. This complex model of the GE/S merger involves the highest possible particle resolution, with GCs acting as reference points to track the flow of the merger. In particular, we are interested in the distribution of GCs in phase space 9–11 Gyr ago. This will enable us to create a more robust sample of GE/S GCs. At the same time, it will address the question of whether GCs in phase-space can provide solid evidence of the host galaxy’s accretion history.

Acknowledgements

The authors thank the anonymous referee for a very constructive report and suggestions that helped significantly improve the quality of the manuscript. M.I. and P.B. thanks the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023. The work was carried out with the support of the Ministry of Science and Higher Education of the Republic of Kazakhstan within the framework of Projects No AP26100669 “Galactic Archaeology of Gaia-Enceladus Globular Clusters as Indicators of the Formation History of the Milky Way Disk”.

References

---

Appendix A Estimation of the GC measurement errors

Here we in detail investigated the error measurements and their influence on the orbital evolution of GCs in the past. In Table A.1 we present the 6D kinematic parameters for selected GCs and their computed relative errors that we take from the catalogue Baumgardt & Vasiliev (2021) and on the web page³. All astrometry data is based on Gaia DR3. It should be noted that the observational data (R_⊙, R_V, µ_α, and µ_δ) have corresponding measurement errors, this is especially true for radial velocity and parallax of the GC’s. Such errors can lead to large uncertainties during reconstruction of the GC orbits on a cosmological timescale. For some of the GCs, the orbit may have large discrepancies in its shape. As can be seen, most GCs have a high precision in their measured kinematic parameters. In general, the relative errors are below 2%. However, we also see some examples for larger errors. For example, NGC 6584 has in µ_α a relative error of 10%. The NGC 7089 has an error in the radial velocity R_V that is 7. 9%.

To evaluate the influence of measurement errors on orbit shape during reconstruction, we generated four additional initial conditions for GCs taking into account the errors for parameters R_⊙, R_V, µ_α, and µ_δ within ±σ. For analysis we choose two extremes. One system, NGC 7099, is an example with a poor measurement accuracy and one system, NGC 5286, is a good example with a high measurement accuracy. In Fig. A.1 we show the result of the orbital evolution for both systems. As can be seen, the shape of the orbits is preserved and does not have a large deviations. Also, in the right panels of the Fig. A.1 we demonstrate the influence of measurement errors for the dynamics of total energy and angular momentum for selected GCs. As can be seen from Fig. A.1, even for the GC’s with the high relative errors in a particular component, this does not lead to a strong variation of the orbital turns over long intervals of backward integration in time.

Appendix B Additional selection parameters

In Table B.1 we provide a detailed list of our GC sample, organised in the same way as in Table 1. The number of articles in which the GC satisfies the GE/S limits (as assumed by every author) in different kinematic parameters is shown in the columns after the slash. In the last row of the Table labelled ‘Chem.’, we show the number of articles in which GCs were considered as GE/S progenitors according to chemical abundances and the age-metallicity relation Ceccarelli et al. (2025); Callingham et al. (2022); Malhan et al. (2022); Kruijssen et al. (2020); Massari et al. (2019), and Massari et al. (2019).

Appendix C Estimation of relaxation time for the GCs with different particle numbers

To investigate the dependence of star loss on particle number, we performed N-body simulations using a set of runs with larger number of particles (doubling the N). So, we created a new GC models with a particle number of N = 100k with the same half-mass radius. We selected three GCs that are most strongly associated with GE/S: NGC 1851, NGC 2298, and NGC 7089. To calculate the current relaxation time, we use the following equation (Binney & Tremaine 2008): $$t_{rel}=\frac{780[\text{Myr}]}{\ln (\lambda \cdot N)}{\left(\frac{M}{{10}^5{M}_{\odot }}\right)}^{1/2}{\left(\frac{r_{hm}}{[\text{pc}]}\right)}^{3/2}\left(\frac{M_{\odot }}{m}\right),$$(C.1) where m – average mass of the single stellar particle, M – total current mass of the cluster, N – number of particles, r_hm – half-mass radius; λ = 0.1 value we also took from Giersz & Heggie (1994). Figure C.1 shows the evolution of the tidal mass of the clusters as a function of the normalised time (i.e. evolution time divided by relaxation time). The relaxation time is calculated at every 10 Myr over a whole computation period. As can be seen for NGC 1851, NGC 2298, and NGC 7089, the mass evolution does not depend on the number of particles in our models. Models, shown with grey curves, represents the GC’s stellar loss in another realisation of the galactic potential; see more details in Appendix D.

Appendix D Some additional test simulations with different conditions

D.1 Evolution of the orbital parameters in the 441327 TVP potential

As stated in Sect. 4.3, to make our result more robust, we investigate the energy vs. angular momentum evolution together with the orbital elements changes in one more external TVP. We selected

441327

external time-variable MW-like potential to determine the possible dependency of the resulting GE/S GC sample from the external potential. We carried out a point-mass integration of 36 GCs in the same way as we described in Sect. 3 but in the

441327

TVP potential. In Fig. D.1 we show the evolution of the total specific energy and the angular momentum of L_z and L_perp for individual GCs. In Fig. D.2 we also present the evolution of the orbital parameters for these GCs.

D.2 Evolution of the orbital parameters in the time-independent potential based on 411321 TVP parameters

For our additional simulations we analyse mass loss from GCs with another realisation of the time-independent potential,

411321 TVP-FIX

. We take the current masses and scale parameters for 411321 TVP as it is today. For disc we assume the following values: M_d = 7.110×10¹⁰ M_⊙; scale length a_d = 2.073 kpc; scale height b_d = 1.126 kpc. For halo: M_h = 1.190×10¹² M_⊙; scale height b_h = 28.48 kpc; see also Fig. 1 in Sect. 3. We fixed these values of external potential parameters and carried out the N-body orbital reconstruction up to −9 Gyr for the three GCs as a single particle: NGC 1851, NGC 2298, and NGC 7089. For comparison of the potential influence on the orbital reconstruction, we show two types of simulations: red represents the orbital evolution in TVP potential, and gray in TVP-FIX; see Fig. D.3. As we see from these plots, applying these two different potentials, we get almost similar orbital shapes only with some minor differences, which. In our opinion, these changes do not generally have a strong influence on the orbits of these GCs. Of course, the TVP with time-dependent parameters is more physically motivated compared to the TVP-FIX for the long term dynamical simulations of GC dynamics in MW.

Also, to investigate the dependence from selected TVP and TVP-FIX we check the dependence from particle number in the models (NGC 1851, NGC 2298, and NGC 7089). For these clusters, we created completely new models with a particle number N = 50k, the same half-mass radii and integrated them forward for 9 Gyr in TVP-FIX. As we see in Fig. C.1, for clusters NGC 1851 and NGC 2298, the models that evolve in TVP and TVP-FIX, variations in mass loss are almost the same. The reason for such a behaviour is probably related to the compactness of the cluster models with r_hm = 2 pc. The NGC 7089 model with a larger r_hm = 4 pc has a slightly different mass loss, but during the entire simulation the difference between the TVP and TVP-FIX mass loss results is not greater than a few percentage points.

D.3 Evolution of the orbital parameters in the 41321 TVP potential with dynamical friction included

In our current work, we apply additional calculations to determine the influence of dynamical friction on our orbital parameter analysis for the GCs. Dynamical friction in astrophysical context indicates the collective deceleration exerted on a moving massive body by the fluctuating force of field stars. The existence of such an effect was first demonstrated by Chandrasekhar and von Neumann (Chandrasekhar & von Neumann 1942, 1943) in their pioneering works. Later, Chandrasekhar in the studies of (Chandrasekhar 1943a,b) developed a more quantitative theory of dynamical friction. The resulting drag acceleration acting on the GC can be written as (Binney & Tremaine 2008) $$\frac{\text{d}{V}_{\text{GC}}}{\text{dt}}=-\frac{4\pi G^2\rho M_{\text{GC}}}{V_{\text{GC}}^3}\chi \cdot \ln \Lambda \cdot V_{\text{GC}}\text{with}\chi =\frac{\rho \left(<V_{\text{GC}}\right)}{\rho }.$$(D.1)

In general, the functions χ and the Coulomb logarithm Λ depends on the velocity of the massive object and the properties of the background system. For more details on this dynamical friction description, we refer the reader to Just et al. (2011). Based on the results of our previous study Just et al. (2011), in our current calculation we use the fixed Coulomb logarithm value ln Λ = 5 and the fixed value for χ = 0.5.

In Fig. D.4 we present the comparison of the evolution of the GC Galactic orbits during 9 Gyr of integration time for the models including dynamical friction and without it. As can be seen from the plot, dynamical friction only slightly affects on orbits. In our opinion, applying the dynamical friction in a present form does not significantly change the orbit of GC and so the position of lost stars from these GCs on the phase-space diagrams.

All Tables

All Figures