Identifying the possible ex-situ origin of the globular clusters of the Milky Way: A kinematic study

This is the second paper in a series, which studies the likelihood that some globular clusters (GCs) of the Milky Way (MW) could have originated from a dwarf satellite galaxy (DSG). Using a large suite of three-body simulations we determine the present-day orbital properties of 154 GCs that could have escaped from 41 MW DSGs over the past $8\,\mathrm{Gyrs}$. For the MW we considered two sets of static and dynamic models which account for the sustained growth of the MW since its birth. We focus on the Magellanic Clouds and Sagittarius. We compare the apogalactic distance, eccentricity, and orbital inclination of the MW GCs with those of runaway GCs from DSGs, to constrain their possible ex-situ origin. We observe a positive correlation between a DSG mass and the dispersion of its runaway GCs in the orbital parameter space of ($R_\mathrm{ap}$, $e$). We provide tables of the identified MW GCs and their likely associated progenitors. In total, we find 29 (19%) MW GCs which could be kinematically associated with MW DSGs. We report, for the first time, 6 and 10 new associations with the Large Magellanic Cloud and the Sagittarius, respectively. For the Sagittarius we predict a concentration of runaway GCs at large apogalactic distances of $R_\mathrm{ap}\approx275-375\,\mathrm{kpc}$, $e\approx0.8$, and a relative inclination of $\Delta\theta\approx20^{\circ}$. So far, there has not been any observed GCs with such orbital elements. Complemented with photometric and spectroscopic observations, and cosmological simulations, the findings from the present study could conclusively settle the debate over the in-situ vs. ex-situ origin of the MW GCs.


INTRODUCTION
According to the standard Λ cold dark matter cosmology, all galaxies are thought to grow in a hierarchical fashion in which larger structures are formed through the continuous merging of smaller structures (White & Rees 1978).If a dwarf satellite galaxy (DSG) merges with a massive galaxy, its stellar populations are expected to be incorporated into the merged system.Cosmological simulations indicate that star formation can be tagged as occurring within the potential well radius of the host primary galaxy (i.e.in-situ) or within the satellite galaxy that is eventually accreted (i.e.ex-situ).This suggests that today's galaxies should contain contributions from both in-situ and ex-situ formed globular clusters (GCs) depending on the galaxy's assembly history (Oser et al. 2010;Pillepich et al. 2015).
More than 150 GCs have been identified around the Milky Way (MW), whose origin is still an issue.There is observational evidence to imply that the MW has been probably made up of mergers with smaller DSGs so that GCs can be considered as potential tracers of this process (Searle & Zinn 1978).A number of GCs in the MW, ★ E-mail: a.rostami@iasbs.ac.ir especially the so-called 'young halo' population, are hypothesized to have been captured from DSGs (Law & Majewski 2010b;Mackey & Gilmore 2004).The ex-situ GCs were born in DSGs and then captured from their host DSG by the MW.The origin of ex-situ GCs can be broadly divided into two categories, 1) GCs that were born in DSGs and their host galaxy is completely merged into the MW; and 2) GCs that originated from a DSG that is orbiting around MW.
Up until today several studies have attempted to associate some of MW GCs with known merging process of DSGs.The Sagittarius dwarf spheroidal galaxy (Sgr) is the first discovered candidate of a merger with the MW (Ibata et al. 1994).Several studies over the past two decades have associated a number of GCs to the Sgr.A number of studies used photometric data and compared the stellar properties of GCs with those of the Sgr stars, such properties include metallicity, chemical abundances, colour-magnitude diagram and age of stars (Carraro 2009;Carraro & Seleznev 2011;Carretta et al. 2014;Sbordone et al. 2015;Carretta et al. 2017;Bellazzini et al. 2002;Caffau et al. 2005;Sbordone et al. 2005).Other studies utilized numerical models that accurately reproduce the position and radial velocity of stars belonging to the Sgr streams, and identified the GCs that may be associated with the Sgr (Law & Majewski 2010b;Bellazzini et al. 2020).Massari et al. (2019) obtained the energy and angular momentum range for GCs that were likely to be associated with the Sgr.The GCs which are thought to be associated with the Sgr are Whiting 1, NGC 6715 (M54), Ter 7, Arp 2, Ter 8, Pal 12, NGC 6284, NGC 5053, NGC 4147, NGC 5634, NGC 5824, Pal 2, and NGC 2419 (Law & Majewski 2010b;Massari et al. 2019;Bellazzini et al. 2020).Among these GCs, M 54, Arp 2, Ter 7, and Ter 8 are still bound to the Sgr (Bellazzini et al. 2020).Koppelman et al. (2019) identified new members of the Helmi streams (Helmi et al. 1999), using the complete phasespace information combining Gaia Data Release 2 (DR2), and the APOGEE DR2, RAVE DR5, and LAMOST DR4 spectroscopic surveys.Moreover, they performed -body simulations to limit the progenitor's properties and timing of accretion.They then labeled GCs that show overlap in energy and angular momentum space with the stream members and concluded that seven GCs are possible candidates to be associated with the Helmi streams (NGC 4590,NGC 5024,NGC 5053,NGC 5272,NGC 5634,NGC 5904,and NGC 6981).
The Gaia-Enceladus-Sausage is an elongated structure in velocity space discovered by Belokurov et al. (2018).using the kinematics of metal-rich halo stars.They showed that it could be created by a massive dwarf galaxy with a total mass of 5 × 10 10 M ⊙ on a strongly radial orbit that merged with the MW.Myeong et al. (2018) sought evidence for the associated Gaia-Enceladus-Sausage GCs by examining the 91 MW GCs' structural data in action space using the Gaia DR2 catalogue, accompanied with the proper motions obtained from the Hubble Space Telescope proper.They identified eight GCs as belonging to the Gaia-Enceladus-Sausage (NGC 1851, NGC 1904, NGC 2298, NGC 2808, NGC 5286, NGC 6864, NGC 6779, and NGC 7089).Massari et al. (2019) studied the origin of 151 GCs of the MW.They concluded that only 40% of the clusters probably formed insitu, while 35% of the GCs are possibly associated with known merging processes including the Gaia-Enceladus-Sausage (19%), the Sagittarius dwarf galaxy (5%), the progenitor of the Helmi streams (6%), and to the Sequoia galaxy (5%).They did not find any origin for other GCs, of these, about 16 per cent of the GCs were classified in the high energy group And the rest of them low energy group.
Using numerical simulations, Khalaj & Baumgardt (2015, 2016) studied the escape of GC stars from the Fornax DSG, in the context of multiple-stellar populations in GCs.They showed that Fornax could have lost a substantial amount of its stellar mass as a result of the MW tidal field and gas expulsion occurring inside GCs.
Rostami Shirazi et al. (2022), hereafter Paper I, determined the escape fraction of GCs from 13 massive DSGs of the MW.They demonstrated that the escape fractions are not negligible, i.e. for most DSGs this value was at least ∼ 20 per cent and for two of them it was above 80 per cent.They concluded that it is very likely that a number of MW GCs originated from DSGs.Given the number of GCs observed in the Fornax, Sgr, Small Magellanic Cloud (SMC), and the Large Magellanic Cloud (LMC), and their corresponding escape fractions obtained by Paper I, they estimated that at least two GCs to have escaped from the Fornax, two GCs from the SMC, four GCs from the LMC, and about 14 GCs have escaped from the Sgr.Moreover, they found that the average escape time of GCs from DSGs reaches a plateau at  ≈ 8 Gyrs which means that the MW DSGs with an age greater than 8 Gyrs, are not likely to lose any more GCs.
The present paper is a continuation of Paper I. Here, we utilize a method based on the kinematic properties of the MW GCs and the MW DSGs.We classify the orbits of the MW GCs and calculate the probability that these GCs are associated with DSGs.Our method does not take into account the photometric/spectroscopic properties of GCs and DSGs.We focus only on the second category of the origin of ex-situ GCs, i.e.GCs that were born in a DSG that is orbiting around MW, and then escaped from it as a result of the MW tidal field.In particular, we focus on the Sgr, LMC, and the SMC, as the heaviest MW satellite galaxies.We also briefly touch upon merging processes and dissolved DSGs.In the future study, we exclusively investigate the association of the MW GCs with merging processes whose progenitor DSGs are completely merged into the MW.
The present paper is divided into four sections.Section 2 includes the description of our methodology and adopted simulation models.In Section 3 we present our results regarding the association of the MW GCs with some DSGs as well as drawing a comparison with other relevant studies.Finally, in Section 4, we summarize the main outcomes of our work, and propose a direction for future studies.

METHODOLOGY
The principal problem we need to solve is obtaining the orbital parameters of GCs that were initially bound to a given DSG and then escaped from it due to the Galactic tidal field.One can model DSGs assuming they are composed of a stellar system of  stars and a number of bound GCs, whose initial conditions are random, albeit in equilibrium.Such an approach requires solving the body problem, i.e. calculating the forces that each individual field star exerts on a GC and then summing over all the forces, for a large number of time steps over a period of ∼ 8 Gyrs.This is computationally expensive, hence not well-fit for our purpose of sweeping a large parameter space.
DSGs are collisionless systems.As a result, one can instead assume a smooth potential field for them.This is an advantage of collisionless systems, allowing us to reduce the -body problem to a three-body problem consisting of the MW, a DSG, and a GC.The usefulness and prominence of such an approach have been previously shown in studies such as Khalaj & Baumgardt (2016).In Section 2.4, we will show that the tail stream distribution of a DSG, which is obtained by our three-body method is in a good agreement with those obtained by more time-consuming direct body method.This is also the same approach we followed in Paper I. We place the MW at the centre of a right-handed Cartesian coordinate system, whose -axis points towards the location of the Sun, the  -axis is in the direction of the Galaxy's rotation, and the -axis is determined by the right-hand rule.The MW is assumed to remain still throughout the simulation owing to its large mass.The trajectory of the DSG is only determined by the MW.We take GCs as point masses whose motion is prescribed by both the MW and the DSG.GCs are spatially distributed according to the density profile of their host DSG.The initial velocities of GCs are drawn from a three-dimensional Maxwell-Boltzmann distribution.We make sure that the GCs are initially bound to their host DSGs.GCs that remain within 2 ×  tidal for about four times their orbital period are considered as initially-bound and their trajectory is followed for a simulation time, using a 10th order Runge-Kutta integrator.At the end of each simulation, we designate those GCs whose final distances exceed 2 ×  tidal from the centre of their host DSG.We consider such GCs as runaways.It should be mentioned that in all simulations we assume that the distribution of GCs follows the distribution of baryonic matter of DSGs.In Section 3.1, we will discuss how the results change if any other initial GC distribution is adopted.Note that only the dark halo component is included in the equation of motion of test particles as it is the most dominant component in dwarf galaxies.There are more (minor) details about the initial conditions of the GCs and the method we use to solve the three-body problem.However, for the sake of brevity, we refrain from repeating them here.One can refer to Paper I for a detailed description of our methodology and the relevant initial conditions.
We perform the simulations for a large ensemble of GCs, i.e. until we end up with 2500 runaway GCs for each DSG and obtain their orbital parameters.Then we compare the orbital parameters of the MW GCs with those of runaway GCs of the given DSG and calculate the probability of their association with the given DSG.Moreover, we examine the effect of the DSG mass on the scattering of GCs in the parameter space, using a semi-analytical method.

The potential field of the MW
For the MW, we adopt the well-known MWPotential2014 model, described in Bovy (2015).In this model, the MW potential field consists of three components, namely bulge, disc, and halo.The bulge density follows a power-law distribution (spherical) with an exponential cutoff given by The disc is a Miyamoto & Nagai (1975) model It also assumes a dark-matter halo with a Navarro-Frenk-White potential (NFW) from Navarro et al. (1997) the potential is time-invariant.In addition to this model, we also consider a dynamic model for the MW potential which is timedependent.This model is motivated by the fact that the mass of the MW had been smaller at birth and has continuously grown since then as a result of e.g.merging processes.To model the sustained growth of the MW over the past 8 Gyrs, we assume that at each point in time, the MW conforms to a potential profile similar to MWPotential2014, where its parameters (masses and scale lengths) are time-dependent.Following Haghi et al. (2015), we can derive the required relations for the masses and scale lengths as a function of time.In particular, the virial mass is given by where  is the cosmological redshift and  c = 0.34.The virial mass at  = 0 (present time) is denoted by  vir (0), and its value is given by MWPotential2014.The values of  vir and  vary as a function of time as given by Haghi et al. (2015), i.e. equations 9 and 6 therein.The relations for the bulge and disc mass as well as their scale lengths are given by five similar equations which can be expressed in a compact mathematical form as follows  {d,b} () =  vir () {d,b} (0) where the {...} notation factors out the variable parts of each equation, while keeping in the identical parts of equations to avoid repetition.Table 1 summarizes the values of the aforementioned parameters (masses and scale lengths) for two epochs of  = 0 Gyr (present time) and  = −8 Gyr.Similar to Fritz et al. (2018a), we consider two values for the halo mass of the MW, i.e  vir = 0.8 × 10 12 M ⊙ (from MWPotential2014), and  vir = 1.6 × 10 12 M ⊙ .All other parameters of the MW potential model remain unchanged.Models with light and heavy halos are denoted by L and H, respectively.

The potential field of the Sgr
We assume a Plummer (1911) model for the present-day potential of the Sgr as follows To determine the values of  Sgr and  sc , we use the results of Vasiliev & Belokurov (2020).They studied the three-dimensional structure of the Sgr using the astrometric and photometric data of Gaia.They found that the total mass of the Sgr enclosed within a radius of 5 kpc is 4 × 10 8 M ⊙ ; and that the peak value of the circular speed of the Sgr is 21 km s −1 which is reached at a radius of ∼ 2.5 − 3 kpc.These findings result in  Sgr = 4.8 × 10 8 M ⊙ and  sc = 1.8 kpc.The half-mass radius ( h ) in a Plummer model relates to  sc via a linear equation, i.e.  h = 1.304 sc .This yields  h = 2.34 kpc for the Sgr.

The time-dependent potential models of the Sgr and the corresponding orbits
We get the present-day equatorial coordinates (, ), proper motions (  cos ,   ), the line-of-sight velocity ( LOS ), and the heliocentric distance ( ⊙ ) of the Sgr from the Gaia DR2 (see Table 5).For each of the potential models adopted for the MW, i.e. dynamic and static models as defined in Section 2.1.1,we trace back the orbit of the Sgr for 8 Gyrs to obtain its initial conditions.In total, we consider eight different models to obtain the orbit of the Sgr.These models are labeled as L1 to L4 (light) and H1 to H4 (heavy).The characteristics of these models are summarized in Table 2.
Models with a static potential for the MW and the Sgr are designated by L1 and H1.The models labeled as L2 and H2 are as same as L1/H1, except that the dynamical friction also enters the equations of motion.We apply the standard form of dynamical friction to the halo (e.g.Binney & Tremaine 2011) where  is the distance of the DSG from the centre of the MW,  is the mass of the DSG, and ln(Λ) = 3 is the Coulomb logarithm.
In the H3 and L3 models, the MW and the Sgr potentials are considered as being dynamic (see Section 2.1.1).In H4 and L4 models, the Sgr and MW potentials are dynamic and the effect of dynamical friction has been considered.Tidal stripping induced by the MW has reduced the total mass (dark + baryonic matter) of the Sgr over the past 8 Gyrs.We have quantified this effect using an -body simulation, where the Sgr is made up of  = 50, 000 equal-mass particles distributed according to the Plummer model.The particles do not undergo stellar evolution.The motion of each particle is determined by the combined potential of all other particles as well as the MW, using an 8th-order Runge-Kutta integrator.We then determine the values of  Sgr ( = −8 Gyr) and  sc ( = −8 Gyr) in such a way that after 8 Gyrs of evolution, they match those of the present-day Sgr.In particular, we obtain an initial dark halo mass of  Sgr ( = −8 Gyr) = 6 × 10 10 M ⊙ for L3, and  Sgr ( = −8 Gyr) = 2.8 × 10 11 M ⊙ for H3.To obtain the initial conditions of the H4 and L4 models, we consider dynamical friction as well.We find  Sgr ( = −8 Gyr) = 4 × 10 10 M ⊙ for L4, and  Sgr ( = −8 Gyr) = 1 × 10 11 M ⊙ for H4, which is in agreement with the results of Gibbons et al. (2017) and Jiang & Binney (2000).Time evolution of the Sgr total mass within the tidal radius in H4 and L4 models is shown in Figure 1.Note that the initial stellar mass of Sgr is estimated to be ∼ 10 9 M ⊙ (e.g.Niederste-Ostholt et al. 2010;Law & Majewski 2010a), which is smaller than the dark halo mass by one to two orders of magnitude.As a result, we ignore the effect of stellar-mass in the equations of motion.
Figure 2 shows the orbit of the Sgr for the aforementioned eight models for the past 8 Gyrs.In H3, H4, L3, and L4 the position of -body particles which constitute the Sgr is known as a function of time.This enables us to calculate the Sgr potential as a function of time and position (ì , ).This time-dependent potential combined with the MW potential is used to determine the trajectory of GCs.

Orbital parameters of GCs
Among all the properties which describe the orbit of a GC we choose the 3-tuple of ( ap , , ì ), where  ap is the apogalactic distance,  is the eccentricity, and ì  is normal to the orbital plane.Instead of ì , one can equivalently use the angular momentum vector ì , since ì  =  n.Within the context of this study, these are sufficient to determine the momentum and energy state of each GC.However, one should note that this 3-tuple does not uniquely determine the orbit of each GC.Moreover, for GCs whose orbital plane as well as the orbital inclination changes, we use the average of ì  vectors over the last 1 Gyr of the simulation.

Probability of GCs association with the DSGs
Using the data from Gaia DR2 (Baumgardt et al. 2019), we can calculate the position and velocity vectors of MW GCs at the present time in our coordinate system.As a result, the orbital parameters of GCs ( ap , , ì ) can be fully obtained.The uncertainties in proper motions (  cos ,   ) and line-of-sight velocity ( LOS ) of GCs and DSGs are usually not negligible and yield large uncertainties in their orbits.In comparison, errors in , , and  ⊙ are negligible.
To examine the correlation of orbital parameters of MW GCs with the runaway GCs of a DSG, we proceed as follows.First, for each of the observed parameters (proper motions and the line-ofsight velocity) of MW GCs, we consider a set ( ) of three data points, which includes the mean value of the parameter along with its corresponding lower and upper bounds, i.e.

𝑆( 𝑝)
where  ∈ {  cos ,   ,  LOS } and   is the corresponding error of each parameter.For each MW GC, this gives us a set  of 3 3 = 27 elements, each element being a different 3-tuple of observed data points.In other words,  =  ( ), where denotes the three-fold Cartesian product of sets since we have three observed parameters.This further translates into a set of 27 3-tuples of orbital parameters ( ap , , ì ) for each GC.It is clear that if the error of data is small, these 27 different orbits of a GC are very close.
Next, depending on whether the DSG is the Sgr, the LMC, the SMC, or else we follow different procedures described in one of the following subsections.

Association with the Sgr
The orbital parameters of the Sgr are well-measured, hence the corresponding errors are relatively negligible.As a result, in all simulations, we only use the mean values for its orbital parameters without considering the errors.
We construct a three-dimensional parameter space whose components are  ap ( kpc), , and Δ • , where Δ is the orbital inclination of a GC measured with respect to the Sgr orbital plane.Each of the Sgr runaway GCs and the MW GCs, occupy a single point in this parameter space.We utilize the multivariate (threedimensional) probability density function (PDF) of runaway GCs in the aforementioned parameter space as the basis to quantify associations.We employed (Gaussian) kernel density estimation to derive the PDF. Figure 3 illustrates the scaled PDF for the Sgr L4 model, showing the isodensity contours.There seems to exist two peaks in the PDF, where one has a spherical-like form (right peak), while the other is more elongated (left peak).Upon close inspection, it becomes evident that the elongated peak is a blend of two sub-peaks.These three peaks correspond to regions of the highest concentration of runaway GCs in the parameter space and coincide with the apogalactic distance of the Sgr as it completes ∼ 3 orbits around the MW and spirals inwards due to dynamical friction over the past 8 Gyrs (ref. Figure 2).As we move away from the peaks, the PDF asymptotically reaches zero, indicating regions where the probability of association is extremely low.To consider the errors into account, we first calculate the PDF values for all elements in the associated set  that we have for each MW GC.We then average over these PDF values to obtain a single PDF value for each MW GC.Based on the PDF values we then classify MW GCs into three categories, namely Flag 1, Flag 2, and nonassociations.Flag 1 corresponds to MW GCs with a very high association probability with the Sgr.In comparison, Flag 2 GCs have a lower association with the Sgr.These respectively refer to all MW GCs which lie within a boundary (contour), encompassing 65% of all runaway GCs for Flag 1, and between 65% and 95% for Flag 2. Non-associated GCs are highly unlikely to have originated from the Sgr as they lie outside the 95% boundary, i.e.where the probability of association is less than 5%.

Association with the LMC and the SMC
The LMC is the most massive DSG of the MW.So far, 19 GCs have been observed in the LMC (Law & Majewski 2010b).According to Paper I, it is estimated that at least 16 per cent of LMC GCs should have been stripped off by the MW, with an average escape time of ≈ 8 Gyrs.This implies that the escape process of the LMC GCs has been completed and we don't expect any more runaway GCs in the future.These results translate into four runaway GCs for the LMC.Assuming that the runaway GCs have not yet been dissolved, they should be still orbiting the MW.
To investigate the association between the MW GCs and the LMC, we repeat the same procedure as we did for the Sgr.We only consider H1 (without the dynamical friction) and H2 (with the dynamical friction) models.We assume a Plummer density profile for the LMC and omitting the (negligible) errors on its orbital parameters.According to the Bekki & Chiba (2009), the total dynamical mass and scale length of the LMC are 2 × 10 10 M ⊙ and 3 kpc, respectively.Two more simulations are performed for SMC in H2 model.One with the presence of the LMC and the other without it.Also for the SMC, we assume a Plummer model with a total dynamical mass of 3 × 10 9 M ⊙ and a scale length of 2 kpc (Bekki & Chiba 2009).However, it should be mentioned that since the LMC and SMC are irregular galaxies, considering Plummer density profile does not necessarily match the observed profiles.We adopted the Plummer model for the density profile of DSGs as it is a good approximation for spherical objects such as globular clusters and galactic dark matter halo.In Section 2.3, we will show that the total mass of DSG is the most effective parameter in the distribution of runaway GCs rather than the exact density profile of the host DSGs.

Association with other MW DSGs
The errors on the orbital parameters and masses of other MW DSGs are not necessarily negligible (see Table 5).In addition, performing simulations for the other 38 MW DSGs is computationally expensive.As a result, the method we used for the Sgr, the LMC, and the SMC cannot be readily applied to other DSGs.For them, we proceed as follows.First, for each DSG and the set of its runaway GCs ( DSG ), we define the filter of an orbital parameter F DSG ( ) as the maximum difference between the value of that parameter for the DSG ( DSG ) and that of its runaway GCs ( GC ), i.e.
where  ∈ { ap , , Δ}, and Δ is the orbital inclination of the runaway GC with respect to the orbital plane of its host DSG.The value of F DSG ( ) gives us a measure of the parameter dispersion in the parameter space for the given DSG.We will demonstrate in Section 2.3.2 that F DSG ( ) is positively correlated with the mass of the DSG.This indicates that for lighter DSGs, the orbit of runaway GCs is very similar to that of the DSG.After the LMC and the SMC, The scaled PDF of the Sgr runaway GCs, colour coded with respect to isodensity contours.The black and white contours correspond to boundaries enclosing 65%, and between 65% and 95% of all data, respectively.These provide the basis for categorizing the MW GCs into Flag 1 (black circles), Flag 2 (black triangles), and non-associations (white circles).In the bottom panel, the red crosses mark the location of the Sgr in the parameter space as the Sgr spirals inwards as a result of dynamical friction.Likewise, the red arrowheads mark the apogalactic distances of the Sgr in the top panel, where we have Δ = 0 • for the Sgr.One can interpret non-associations as GCs with an association probability of less than 5%.The apparent placement of some points (e.g.non-associations) within the white or black contours is due to the fact that panels display projections.As a result, one needs to consider both panels together in order to infer the association categories.The elongated peak on the left consists of two sub-peaks, i.e. there exists three peaks in unison with the ∼ 3 apogalactic passages of the Sgr around the MW (see Figure 2).
the Sgr is the heaviest DSG of the MW.The values of F DSG ( ) obtained for these three DSGs can be considered as an upper limit on F DSG ( ) for all other (lighter) DSGs.Among them, the Sgr provides the most stringent bounds.Therefore, we pick F Sgr ( ) as a measure to compare the orbit of MW GCs and runaway GCs of a DSG.This alleviates the need for performing more simulations similar to those of Section 2.1.3.
Having found an estimation on F DSG ( ) for all MW DSGs, we check if both of the following conditions are satisfied for all pairs of a MW GC and a DSG where we would like to reemphasize that GCs refer to MW GCs and not the runaway GCs of DSGs.We will show in Section 2.3.3, if the mass and the orbital eccentricity of a DSG is fixed and it is placed in an orbit with larger values for  ap , the corresponding F DSG () and F DSG (Δ) will not change.However, F DSG ( ap ) grows as  ap (DSG) increases, hence condition  in Equation ( 11) is defined.
As mentioned earlier, unlike the Sgr and the LMC, the uncertainties on the orbital parameters of other DSGs are not negligible and cannot be simply omitted.Similar to the set  of 27 initial conditions that we made for each GC, we generate 27 initial conditions for each DSG.This translates into 27 × 27 = 729 sets of orbital parameters for each pair of a MW GC and a DSG, for each of which we check if the conditions given by Equation ( 11) hold.
Finally, for each pair of a MW GC and a DSG, we define the association probability as follows where  true is the number of elements in the set of all orbital parameters for which Equation ( 11) holds and 729 is the total number of elements.Based on this probability we assign association flags to GCs and DSGs, as we did in Section 2.2.In particular, a probability of P ≥ 0.6 is designated as Flag 1 indicating a high association probability, and Flag 2 with 0.2 ≤ P < 0.6 corresponds to a lower association probability.

dichotomy in the semi-major axes of runaway GCs
We consider the MW, the DSG, and the GC as point masses whose masses are ,  1 , and  2 , respectively.For simplicity, we restrict their motion to two dimensions only, i.e. the  −  plane.Initially, the GC is bound to the DSG and forms a binary with a semi-major axis of  12 and an eccentricity of  12 .In addition, the common centre of mass (CM) of the DSG-GC system orbits the MW, with an eccentricity and semi-major axis of  CM and  CM , respectively.
Once the GC becomes unbound, the DSG and the GC have their own orbits around the MW with orbital parameters of ( 1 ,  1 ) and ( 2 ,  2 ), respectively.The conservation of energy yields where Since the motion of objects is on the  −  plane only, the angular momentum is in the  direction ( ì  =  ẑ).The angular momentum for a two-body system is where  and  are the reduced mass and the total mass, respectively.The conservation of angular momentum yields where Given the values for {,  1 ,  2 ,  CM ,  12 ,  CM ,  12 ,  1 ,  2 } one can solve the system of equations 13 to 16 for unknowns  1 and  2 .As an example, suppose  = 1.6 × 10 12 M ⊙ (H1 model),  1 = 10 8 M ⊙ , and  2 = 10 5 M ⊙ .The GC initially orbits the DSG with ( 12 = 1 kpc,  12 = 0).Moreover, the DSG-GC system initially orbits the MW on an orbit with ( CM = 50 kpc,  CM = 0).Owing to the negligible mass of the GC compared to the mass of DSG, it has a marginal effect on the DSG orbit upon becoming unbound, hence  1 ≈  CM .We assume that the orbit of the GC after the escape is circular ( 2 = 0).Solving the system of equations, yields two solutions for unknowns  1 and  2 .One solution is { 1 = 49.9958, 2 = 54.3845},describing an orbit on which the cluster reaches higher distances.The other solution is { 1 = 50.0042, 2 = 46.0204}corresponds to an orbit with a smaller semi-major axis for the GC.This is a pattern that we observe in all of our simulations, i.e. there exist two populations of runaway GCs, scattered either closer to or further from the MW, with respect to their host DSG.This dichotomous pattern is due to GCs escaping in the vicinity of the two Lagrange points L 1 and L 2 .
It is worthwhile to mention that for some given values, these equations do not have any solutions.For example, consider the case of  CM = 0.5 and  2 = 0 for which the equations are unsolvable.In Section 2.3.2,we will examine the effect of  1 and  2 on the distribution of GCs.

The effect of DSG mass on the distribution of orbital parameters of runaway GCs
To illustrate the effect of the DSG mass on the distribution of orbital parameters of runaway GCs, we assume that the values of ,  2 ,  CM ,  12 , and  12 are given and equal to those of Section 2.3.1.
If we consider two values for  CM ∈ {0.0, 0.5}, then for different masses of a DSG with different eccentricities of GCs, we obtain two values for  2 .Figure 4 shows the values of  2 for each  2 in the range of 10 8 to 10 10 M ⊙ for  1 .As  1 increases, the dispersion in the semi-major axis ( 2 ) of GCs around the MW increases.
As an example, for  CM = 0.0 and  2 = 0.3, if  1 = 10 8 M ⊙ , then  2 ∈ {36, 75} kpc, whereas for  1 = 10 10 M ⊙ we have  2 ∈ {25, 114} kpc.As another example, if  CM = 0.5 (the top panel of Figure 4) and 1 < 10 10 M ⊙ , the equations do not have any solutions for  2 and  1 if  2 < 0.5.In all these calculations, we assume  2 = 10 5 M ⊙ .If we set  2 = M ⊙ , the values of  2 do not change much.This means that the distribution of GCs that escaped from DSG is similar to the field stars that possibly escaped from it (cf.Khalaj & Baumgardt 2016).
Moreover, we consider two DSGs with masses of 10 7 and 10 9 M ⊙ .The potential profile of these DSGs is the Plummer model, and we determine their scale lengths () such that the mass of the DSG divided by  3 is equal to 0.5 M ⊙ pc −3 .This is a reasonable density for the MW DSGs, because according to Table 3 of Paper I, the average total mass density of DSGs is close to 0.5 M ⊙ pc −3 .We placed these DSGs at a distance of 30 kpc on the MW disc and selected their velocities such that their eccentricities are  = 0.5.We use the H1 model for the simulation and obtain the orbits of 1000 runaway GCs from each of the DSGs after 8 Gyrs.The orbital parameters of these GCs can be seen in Figure 5.It is evident that the dispersion of orbital parameters of runaway GCs increases as the mass of the DSG increases.In other words, for heavier DSGs, the orbits of escaped GCs have a lower similarity to their host DSG orbit.The dynamical masses of most MW DSGs are less than 10 8 M ⊙ .As a result, the orbits of the escaped GCs are expected to be very similar to the orbit of their host DSG.
The Sgr is one of the DSGs that is close to the centre of the MW.Most of the other DSGs are located at larger radii, and consequently, the MW has stripped a smaller fraction of their initial mass in comparison.As a result, their present-day mass is not significantly different from their initial mass.This means that the choice of a static potential is a good approximation for simulating these DSGs.Moreover, one can argue that dynamical friction has a marginal effect on these DSGs due to their low mass, compared to e.g. the Sgr or the LMC.We showed that there is a positive correlation between the dispersion of the orbital parameters of the runaway GCs, and the mass of their host DSGs.Combined with the fact that the Sgr is heavier than other DSGs (except for the LMC and the SMC), we can conclude that the Sgr filters in the static models can be considered as a high limit for the filters of other DSGs.This justifies our methodology in Section 2.2.3.It should be noted that repeating this test for the lower-mass DSGs leads to much narrower filters.Future investigations based on more accurate data of density profiles and orbital parameters of other 38 MW DSGs can improve the results.

The effect of DSG orbital distance on the distribution of orbital parameters of runaway GCs
We assume a DSG with the mass of 10 9 M ⊙ , an orbital eccentricity of  = 0.5, and two different apogalactic distances of  ap = 30 and 100 kpc.Similar to Section 2.3.2,we use the H1 model and obtain the orbits of 1000 runaway GCs after 8 Gyrs. Figure 6 shows the orbital parameters of runaway GCs.A comparison with the right panel of Figure 5, shows that as  ap (DSG) increases the values of F DSG () and F DSG (Δ) do not change.In the case of  ap = 30 kpc, F DSG ( ap ) is ∼ 20 kpc, while for  ap = 100 F DSG ( ap ) increases to approximately 60 kpc.This indicates that the values of  ap (DSG) and F DSG ( ap ) are positively correlated.

Comparison of direct N-body and our large ensemble of three-body simulations: distribution of tail stream of DSGs
First, we use NBODY6 (Aarseth 2003) direct -body code, to obtain orbital parameters of runaway stars of DSGs due to the Galactic tidal field, assuming three DSG masses of 10 7 , 10 8 , and 10 9 M ⊙ .All modelled DSGs embody  = 10 5 equal-mass particles.The initial positions and velocities of the particles in the DSGs are set such that the mass-density obeys a Plummer profile in virial equilibrium.We determine the corresponding scale lengths () such that the mass of the DSG divided by  3 is equal to 0.5 M ⊙ pc −3 .The particles do not undergo stellar evolution.The DSGs move on an orbit with  = 0.5 through a host galaxy that consists of three components, a central point-mass with  = 1.5 × 10 10 M ⊙ , a Miyamoto & Nagai (1975) disc potential with numerical constants  d = 5 × 10 10 M ⊙ ,  = 4 kpc and  = 0.5 kpc (see Equation ( 2)), and a logarithmic potential for the halo as follows The constant  c is chosen such that the combined potential of the three components yields a circular velocity of  ∞ = 220 km s −1 in the disc plane at a distance of 8.5 kpc from the galactic centre.
We placed these DSGs at a distance of 30 kpc on the MW disc.At the end of simulations, we obtain the orbital parameters of runaway particles, i.e. particles of which final distances exceed 2 ×  tidal from the centre of their host DSG.Next, using our three-body method we perform a simulation for the H1 model, in which the MW potential is adjusted to match that of NBODY6 and obtain the orbital parameters of 2500 runaway particles from the DSGs.
Figure 7 compares  ap and  of runaway particles from modelled DSGs as determined by NBODY6 and our threebody methods.Evidently, the resultant distributions from different methods perfectly match.In particular, the region that runaway particles occupy in the orbital parameter space expands as the DSG mass increases.The observed consistency between the two approaches extends beyond  ap and , and includes the inclination of runaway particles as well.As a result, our large ensemble of three-body methods are an acceptable alternative for direct -body methods which are computationally expensive, and therefore can be readily utilized to obtain the distribution of tail streams of DSGs.

Runaway GCs of the Sgr
Figure 8 shows the orbital parameters of 2500 runaway GCs in each of the adopted simulation models.It evidently depicts two populations of GCs, one with  ap ≤  ap (Sgr) and another one with  ap >  ap (Sgr).This has been already shown semi-analytically in Section 2.3.This dichotomy resembles a butterfly pattern in which each wing corresponds to one of the aforementioned populations.Interestingly, one can see that the runaway GCs occupy a wider range in the dynamic models (the bottom two rows), compared to static ones (the top two rows).This is due to the larger initial mass of the Sgr in dynamic models, in which some of the runaway GCs can travel as far as 400 kpc away from the MW, while some others can reach the MW bulge.Moreover, the orbital inclination can be as large as 170 • .This implies a flip in the angular momentum vector, meaning that some runaway GCs have retrograde orbits around the MW.This can be compared with the static models (L1, L2, H1, H2), where the maximum value of the orbital inclination is about 17 • .
In the H4 and L4 models, where dynamical friction is also at play, the Sgr starts with a larger value of  ap and .It then spirals inward and drags (runaway) GCs along with itself.This can be seen in the bottom row of Figure 8, where the distribution of runaway GCs follows the Sgr trajectory.As expected, the presentday distribution of orbital parameters of runaway GCs obtained from the static model (L2) over 8 Gyr yields similar results to that of the dynamic model (L4) over the last 1 Gyr.We also investigated the possible interference of the LMC on the orbital parameters of the Sgr runaway GCs.We considered masses of 10 10 to 10 11 M ⊙ (Erkal et al. 2019) with and without dynamical friction for the LMC.The distribution of the Sgr runaway GCs remained almost unchanged in all these cases due to the large distance of the Sgr from the LMC.
To calculate the equations of motion for GCs, we need to consider the total mass (stellar mass of the galaxy + dark-halo mass) of DSGs.The GCs follow the baryonic component of the DSGs, i.e. the GC system size is close to the DSG size (e.g.Caso et al. 2019).Estimating the initial size of the stellar component of a DSG is not trivial.As a result, we consider three different radii for the distribution of GCs to investigate the sensitivity of the results to this choice.Figure 9 shows the present-day distribution of orbital parameters of runaway GCs for L4 model assuming three different values for the ratio of barynoic to dark-halo scale lengths, for the initial Plummer distributions of GCs, i.e.  sc,bar =  sc,halo , where  ∈ {0.25, 0.5, 1.0}.Owing to the fact that GCs located at large radii from the centre of their host DSG are only loosely bound, one might expect a larger dispersion of orbital parameters in the parameter space for larger  values.Interestingly, this is not what we observe in our simulations (Figure 9).Despite the fact that the escape rate of GCs are significantly different in the three  cases, their distribution in  ap and  is similar (Figure 9).This is due to the initial bounding condition of GCs (Section 2).In particular, for  = 1.0 a large fraction of GCs are initially located in the outer part of the DSG.Such GCs do not satisfy the initial bounding criterion, i.e. they do not remain within 2 ×  tidal for about four times their orbital period.Therefore, they are not counted as runaway GCs and excluded.In other words, GCs that ultimately contribute to the distribution of runaway GCs, form a sub-system with an effective  value which is less than 1.0.On the other hand, for  = 0.25, the GCs mainly reside in the inner region of the DSG, which has a high escape speed.Therefore, runaway GCs emanating from these areas must have a high speed.This leads to a more scattered distribution of orbital parameter for such DSGs, i.e. such DSGs have a larger effective  value.Therefore, cases of  = 1.0 and  = 0.25 yield distributions which are only marginally different from  = 0.5, if any at all.As a result, the distribution of orbital parameters of runaway GCs does not depend on the initial distribution of GCs and the most important factor in the dispersion of runaway GCs is the total mass of the DSG.In addition, in Section 3.3 we demonstrate that it does not even depend on the adopted density profile.Hereafter, we assume  = 0.5 for all simulations.

GCs associated with the Sgr
Using the runaway GCs of the Sgr given in Section 3.1 and the method described in Section 2.2, we determine the probability of association between MW GCs and the Sgr. Figure 10 shows the orbital parameters of the Sgr runaway GCs in H1 and L1 models, as well as those of 13 MW GCs, studied in Massari et al. (2019), Bellazzini et al. (2020), andLaw &Majewski (2010b), considering their uncertainties.These MW GCs are believed to have escaped  from the Sgr.The Sgr orbital uncertainty is very small, hence it occupies a small area in the parameter space.As evident from the figure, many of the MW GCs that are thought to have originated from the Sgr, are far away from the Sgr runaway GCs, in terms of the orbital parameters.In particular, NGC 6284, NGC 5053, NGC 4147, NGC 5824, and Pal 2 are highly unlikely to have originated from the Sgr.For the H1 model, only six GCs have similar orbits to the Sgr runaway GCs.These are Whiting 1, NGC 6715 (M54), Ter 7, Arp 2, Ter 8, and Pal 12. NGC 5824 lies within the range of Sgr filters in terms of  and  ap .However, its orbital inclination is ∼ 100 • , whereas the maximum orbital inclination of the Sgr runaway GCs in the H1 is only 17 • .In the L1 model, in addition to the GCs mentioned for the H1, NGC 2419 and NGC 5466 are candidates for being associated with the Sgr.Table 3 lists the possible association of MW GCs with the Sgr in different models.In this list, we are reporting several Sgr-GCs associations for the first time.The dynamic models yield more associations as a result of their wider filter distributions.
In agreement with other studies, Whiting 1, NGC 6715 (M54), Ter 7, Arp 2, Ter 8, Pal 12, and NGC 2419 are likely to be associated with the Sgr.Moreover, our current study indicates that NGC 5466 has the possibility of an association in all L models as well as H3 and H4 models.In addition, Rup 106, NGC 4590, NGC 5634, and Pal 4 are likely to be associated with the Sgr in all dynamic models.NGC 6101, NGC 5897, NGC 6235, NGC 6934, NGC 6426, and IC 4499 are other GCs likely to be associated with the Sgr, albeit with a lower probability.Pal 4 is one of the GCs that has large orbital energy with an eccentricity of  ≈ 0.9, and our results indicate that it belongs to the Flag 2 category.This is in agreement with Figure 8, where one expects GCs with  ap ≈ 300 kpc and an  ≈ 0.9 to have originated from the Sgr.
Assuming that M 54, Arp 2, Ter 7, and Ter 8 are still bound to the Sgr (Bellazzini et al. 2020), Paper I predicted that about 14 GCs could have escaped from the Sgr.Here we have introduced 18 GCs that are likely to be associated with the Sgr.
It should be noted that to calculate the orbit of the MW GCs, we have considered only the gravitational potential of the Galaxy and ignored the potential of other DSGs.As also shown in Figure 10, this explains the reason why the orbital parameters of M 54, Arp 2, Ter 7, and Ter 8 lie within the range of orbital parameters of runaway GCs but they are not inside the Sgr.If we also consider the potential of Sgr for a cluster such as M 54, its location in Figure 10 will remain close to the Sgr after the evolution.This is in agreement with the fact that M 54 is thought to be the nuclear cluster of Sgr (Ibata et al. 1994;Bellazzini et al. 2008).
As a final remark, our simulations predict a concentration of runaway GCs for the Sgr at large distances of  ap ≈ 275 kpc for the L4 model and  ap ≈ 375 kpc for the H4 model.Both of these concentrations have high eccentricities of  ≈ 0.8 and relative inclinations of Δ ≈ 20 • .These are not associated with any observed MW GCs which is an interesting discovery.We further elaborate upon this in Section 4.

GCs associated with the LMC and the SMC
Figure 11 shows the distribution of the LMC runaway GCs in H1 and H2 models.Similar to the Sgr, the dichotomous pattern can be seen in the figure.In particular, a population of GCs can be seen in the range of [30,80] kpc.These GCs have a larger orbital inclination compared to the GCs at larger distances.Moreover, compared to H1 model, the distribution of the orbital parameters of GCs expands more in H2 model in which the dynamical friction is considered in the orbital history of the LMC, causing the LMC orbit to spiral inwards.
We should point out that the main caveat of our approach is that the LMC and SMC are assumed to be in an equilibrium state by adopting a simplified Plummer density profile, which does not necessarily match the observed profiles.In order to show the sensitivity of the results to the adopted density profile, new calculations were performed based on the Hernquist model for the density profile while both mass and half-mass radius remained consistent with Bekki & Chiba (2009).As illustrated in Figure 12, there is not a significant difference in the orbital parameters of runaway GCs.Therefore, we conclude that the total mass of a DSG is the most effective parameter in the distribution of its runaway GCs rather than the density profile of the host DSGs.
Figure 13 shows the LMC orbit and two runaway GCs in H1 model.This figure implies that there could be GCs that have escaped from the LMC and are likely to be found at smaller distances from the MW.In the same manner, there could be runaway GCs located at very large distances from the MW.It is clear from Figure 11 that despite the wide distribution of orbital parameters of runaway GCs, their maximum orbital inclination is less than 38 • .This condition is a strong constraint, leading to only a few MW GCs meeting the criteria for association with the LMC.LMC, the association of MW GCs with the SMC alone, cannot be studied in the context of this study.

GCs associated with other MW DSGs
As explained in Section 2.2.3 we use the Sgr filters for other DSGs as well, except for the LMC and the SMC.In particular, we consider the Sgr filters for H1, H3, L1, and L3 models.For the orbital eccentricity  4. The table shows that there exists more than one ex-situ origin for a number of GCs.As a result, there are GCs for which we cannot reliably pinpoint only one DSG for their origin.We refer to these cases as multiple identifications.However, the category of association likelihoods, i.e. flag numbers, can rectify this issue to some extent.
Paper I estimated that about two GCs should have escaped To identify the possible progenitors of MW GCs, Boldrini & Bovy (2022) performed a set of comprehensive orbit integrations to track 170 GCs and 11 MW DSGs backwards in time in a combined potential of the Milky-Way and satellites.To evaluate possible past associations, they proposed a globular-cluster-satellite binding criterion based on the satellite's tidal radius and escape velocity.They found that only 6 GCs were effectively associated with the Sgr.NGC 6715, Arp 2, Terzan 7 and Terzan 8 still belong to Sgr, whereas Whiting 1 and Palomar 12 were accreted by the MW less than 0.3 Gyrs ago.According to their results, the other 164 GCs, have not been associated with MW DSGs.The difference between our results and their conclusion could be due to the binding criteria that is used in Boldrini & Bovy (2022).They assumed that the orbital radius of a GC centered on a galaxy at a specific time  c GC ( i ) to be smaller than the satellite tidal radii at  i , and the GC relative velocity with respect to its putative parent galaxy to be lower than the escape velocity of the satellite at  i .It should be noted that in the binding criterion method of Boldrini & Bovy (2022),  c GC () parameter is highly sensitive to the assumed model of MW-mass growth as well as DSGs mass-loss histories and dynamical friction, while the orbital parameter space for runaway GCs is not sensitive to the initial conditions of the host DSGs and MW models.It is expected that the binding criterion method (Boldrini & Bovy 2022), will give correct results in short periods of time (0.5 Gyr).It will be important to continue testing both methods with Gaia DR4 and further fullscale -body modeling of GCs evolving within a dynamic triaxial galactic potential of DSGs and the MW to reach a consensus.

Runaway GCs from dissolved DSGs
As mentioned in Section 1, we expect that a number of MW GCs come from DSGs that have been completely merged with the MW.Such DSGs must have two important characteristics.First, they must be massive enough to form GCs, and second, they must be close enough to the MW to merge with it completely despite their large masses.
In this section we study the distribution of orbital parameters of runaway GCs from a dissolved DSG.For this purpose, as we did for the Sgr in the H3 models, we must first solve the -body problem to obtain the potential field of the DSG as a function of position and time, i.e. (ì , ).We perform the simulation from  0 = −12 Gyrs until the present time.As described in Section 2.1.3for the Sgr, we consider 50,000 particles of the same mass distributed according to the Plummer model with a total mass of 5×10 10 M ⊙ and a half-mass radius of 4 kpc.We place the DSG on the MW disc with an orbit of  ap = 24 kpc and  = 0.2.As the MW potential grows stronger, the eccentricity of the DSG orbit does not change significantly, whereas  ap reaches 17 kpc at the time of its dissolution, which is  ≈ 6 Gyrs.Now that (ì , ) is known, we proceed to determine the orbital parameters of runaway GCs at the present time.Figure 15 shows the distribution of orbital parameters of these GCs.
As shown in the figure, the orbits of escaped GCs cover a very wide range and there is almost no correlation between the orbital parameters of these GCs and the DSG.Due to the strong gravitational force of the DSG and its proximity to the MW, the orbits of GCs get severely disturbed, leading them to be thrown in different directions.Another interesting point is that there are runaway GCs with very large apogalactic distances.This indicates that even GCs that are very far from the MW, could come from a dissolved DSG.In Section 2.3.2 we showed that for light DSGs, the distribution of orbital parameters of GCs covers a smaller range.However, it should be noted that the initial mass of these DSGs can not be considered too low, since they can not form GCs. The wide distribution of GCs orbits makes it very difficult to find their origin (host DSG).As an example, the origin of a GC with  ap = 200 kpc and  = 70 • , can be the same as that of a GC located on the MW disc at a distance of  ap = 8 kpc. Figure 16 shows the presentday spatial distribution of 500 GCs escaped from the dissolved DSG.As evident from the figure, the runaway GCs are widely distributed in space.For such GCs, complementary photometric and spectroscopic observations can help in identifying their possible ex-situ origin.
Using the data from Gaia DR2 and the proper motions from the Hubble Space Telescope, Myeong et al. (2018) studied 91 MW GCs in energy-action space, to investigate their possible association with the Gaia-Enceladus-Sausage structure.They showed that eight of the high-energy, old halo GCs are strongly clustered in both vertical and azimuthal axes, but are well spread out radially.In addition, the GCs exhibit a large radial anisotropy with highly eccentric orbits, i.e.  > 0.8.Their apogalactic distances and orbital inclination range from 12 to 20 kpc, and 5 to 170 • , respectively.They concluded that the characteristics of these GCs are consistent with a DSG merged with the MW, which could be the origin of the Gaia-Enceladus-Sausage.They estimated an initial mass of ∼ 5 × 10 10 M ⊙ for the DSG.We argue that their compiled list of candidates could be incomplete.We demonstrated in Figure 16 that for a dissolved DSG, the runaway GCs occupy a wide area in the parameter space.For example in Figure 16, the eccentricity of the runaway GCs can be any value from 0.0 to 1.0.This is in contrast to the results of Myeong et al. (2018), where all their candidates have  > 0.8.This implies that they could be, in principle, excluding some GCs.However, due to different methods and selection criteria, we cannot readily compare their results with ours and further elaborate on this discrepancy.This could be the subject of a future study, where we exclusively investigate the association of the MW GCs with merging processes whose progenitor DSGs are completely dissolved.

CONCLUSION
As a follow-up to Paper I, we carried out a large number of threebody simulations to obtain the distribution of possible runaway GCs from the MW DSGs.To obtain the distribution we assumed GCs as point masses and DSGs as a Plummer model (Section 2.1.2).For the MW we considered two sets of static and dynamic models which account for the sustained growth of the MW since its birth (Section 2.1.3).The MW models are constituted by three proved the existence of such a dichotomy using a semi-analytical approach.Moreover, the values of  ap and  for runaway GCs are positively correlated with each other.In other words, GCs with larger (smaller) values of  ap are expected to have larger (smaller) values of  on average or vice versa.
• There exists a positive correlation between the mass of a DSG and the dispersion of its runaway GCs in the parameter space.This implies that for massive DSGs such as the LMC, an associated runaway GC could be as far as  ap = 400 kpc or as close as  ap = 30 kpc.In contrast, the values of  ap and  for runaway GCs of a light DSG such as Carina I are expected to lie within a very small range from  ap and  of Carina I.This correlation renders the classical methods of identifying associations between MW GCs and DSGs error-prone.In other words, one cannot simply look for clusters in the parameter space and associate them with a neighboring DSG.For massive DSGs there could be a large tail of GCs in the parameter space which are not clustered around the DSG and their association will be probably ruled out.
• Related to the previous point, the dispersion in the parameter space for the runaway GCs of dissolved DSGs is even larger.The runaway GCs occupy a large fraction of the parameter space for dissolved DSGs.For a dissolved DSG with ( ap = 24 kpc,  = 0.2), the apogalactic distances of runaway GCs can be any value from  ap = 5 to  ap = 300 kpc and their orbital eccentricities can range from  = 0.0 to  > 0.9.The only statement which can be made about such DSGs is that their runaway GCs still exhibit the aforementioned positive correlation between  ap and .Apart from this, there are not any other patterns or correlations between the distribution of runaway GCs and their dissolved host DSG.This makes it extremely difficult, if not entirely impossible, to identify associations for dissolved DSGs, within the context of purely kinematic studies like ours.
• For DSGs that are close to massive DSGs such as the Sgr or the LMC, the orbital parameters of runaway GCs can be significantly perturbed by the massive DSG, to the extent that, some of the GCs might be even recaptured by the massive DSG.This is why for the SMC, we needed to take the presence of the LMC into account.
• We identified 18 GCs which are associated with the Sgr.Combined with the fact that the Sgr still has four bound GCs, this is in agreement with Paper I where we anticipated about 14 GCs should have escaped from the Sgr so far.According to a number of studies the following five GCs have been suggested to have originated from the Sgr, NGC 6284, NGC 4147, NGC 5053, NGC 5824, and Pal 2. However, we could not confirm their kinematic association.We report a number of new associations with the Sgr for the first time.They are NGC 5466, NGC 4590, Rup 106, and Pal 4, all which have a high probability of association.Moreover, we identified NGC 6101, NGC 5897, NGC 6235, NGC 6934, NGC 6426, and IC 4499, albeit with a lower probability.Other associations we identified are Whiting 1, NGC 6715 (M 54), Ter 7, Arp 2, Ter 8, Pal 12, NGC 2419, and NGC 5634, all of which are in agreement with other studies.
• In Paper I we predicted four possible associations with the LMC.Here, for the H1 model, the high probability associations are NGC 5024, Pyxis, and Pal 3, and lower probability ones are Pal 4, NGC 7492, and NGC 5053.For the H2 model, the results are the same as H1, except that Pal 4 is added to Flag 1 associations as well.Intriguingly, Massari et al. (2019) has categorized Pal 3 and Pyxis in the high-energy group and their origin as uncertain.
• According to Paper I, we expect two runaway GCs from the SMC.In the present paper we demonstrated that the LMC has a non-negligible influence on the dispersion of the SMC runaway GCs.In particular, we observed that the SMC runaway GCs could be recaptured by the LMC.As a result, we cannot conclusively constrain the SMC associations in the context of our current study.
• For the rest of the MW DSGs, which have lower masses, we designated 19 associations, of which eight are of Flag 1, and 11 of Flag 2. For the Fornax, we had predicted two runaway GCs (Paper I).We found Crater as a possible association with the Fornax.Furthermore, another Fornax runaway GC might now reside in the LMC (Mucciarelli et al. 2021).This is consistent with our findings for the Fornax.
• In total, we identified 29 MW GCs which could have originated from DSGs.This indicates that a maximum of 19 per cent of all MW GCs could have an ex-situ origin.
• For a number of GCs, we found several ex-situ origins.We refer to these cases as multiple identifications.Our categories of association likelihood, aka flags, can rectify this issue to some extent but not completely.As a result, there are GCs for which we cannot reliably pinpoint only one DSG and their origin still remains an open question.
• Finally, we find a concentration of simulated runaway GCs from the Sgr which are clustered around  ap ≈ 275 kpc for the L4 model, and  ap ≈ 375 kpc for the H4 model, both of which have  ≈ 0.8 and Δ ≈ 20 • .These values correspond to the orbital conditions of the Sgr at the beginning of the simulations, i.e. at  = −8 Gyr, where the Sgr was located at larger distances due to the effect of dynamical friction.So far, there are not any MW GCs observed with such orbital parameters.This could be due to an observational bias since such GCs might not be sufficiently bright to be easily detected.Another possibility is that such GCs have entirely dissolved owing to their low initial mass.This in principle can lead to a dispersed stream with a low surface brightness which also poses an observational challenge.A third option could be that such GCs

Figure 1 .
Figure 1.Time evolution of the Sgr total mass within the tidal radius in H4 and L4 models.

Figure 2 .
Figure 2. The orbit of the Sgr DSG in our eight simulation models for a period of 8 Gyrs.Orbits corresponding to light and heavy halo masses are shown in red and black, respectively.The orbits are depicted for the  −  plane of the Galactocentric coordinate.

Figure 4 .
Figure 4.The semi-major axis ( 2 ) of GCs around the MW versus their orbital eccentricities ( 2 ) in the range of 10 8 to 10 10 M ⊙ for  1 .The DSG-GC system initially orbits the MW on an orbit with  CM = 50 kpc and eccentricities of  CM = 0 (bottom panel) and  CM = 0.5 (top panel).

Figure 5 .
Figure 5. Orbital parameters of 1000 GCs escaped from two DSGs located on the MW disc after 8 Gyrs in the H1 model.The eccentricities of DSGs are 0.5.The colour coding represents the orbital inclination of GCs with respect to the DSG.The left and right panels correspond to DSG masses of 10 7 and 10 9 M ⊙ , respectively.

Figure 6 .
Figure 6.Same as the right panel of Figure 5 but for DSG orbits at a higher apogalactic distance  ap = 100 kpc.

Figure 7 .
Figure 7. Distribution of  ap and  of runaway particles from three DSGs located on the MW disc with orbital parameters of  ap = 30 kpc and  = 0.5, obtained from NBODY6 (yellow) and our three-body method (black), for three DSG masses of 10 7 , 10 8 and 10 9 M ⊙ .

Figure 8 .
Figure8.The orbital parameters of 2500 runaway GCs, from the Sgr in eight simulation models.The horizontal axis is the apogalactic distance ( ap ) in kpc and the vertical axis is the orbital eccentricity of these GCs.The colour coding corresponds to the orbital inclination of the GCs in degrees, with respect to the orbital plane of the Sgr.The large squares mark the Sgr.For models with dynamical friction, the position of the Sgr in the parameter space changes as a function of time.This has been shown by the magenta-filled square for the present day, and the white-filled squares for earlier epochs.For H2, H4, L2, and L4 models, the orbital parameters of the Sgr change over time as a result of dynamical friction.

Figure 9 .
Figure 9.The present-day distribution of orbital parameters of runaway GCs for the L4 model assuming three different values for the ratio of baryonic to dark-halo scale lengths, for the initial Plummer distribution of GCs, i.e.  sc,bar =  sc,halo , where  = 1.0 (black),  = 0.5 (green), and  = 0.25 (yellow).The distribution of orbital parameters of runaway GCs is not sensitive to the adopted initial GCs spatial distribution.
Figure11shows the distribution of the LMC runaway GCs in H1 and H2 models.Similar to the Sgr, the dichotomous pattern can be seen in the figure.In particular, a population of GCs can be seen in the range of[30, 80]  kpc.These GCs have a larger orbital inclination compared to the GCs at larger distances.Moreover, compared to H1 model, the distribution of the orbital parameters of GCs expands more in H2 model in which the dynamical friction is considered in the orbital history of the LMC, causing the LMC orbit to spiral inwards.We should point out that the main caveat of our approach is that the LMC and SMC are assumed to be in an equilibrium state by adopting a simplified Plummer density profile, which does not necessarily match the observed profiles.In order to show the sensitivity of the results to the adopted density profile, new calculations were performed based on the Hernquist model for the density profile while both mass and half-mass radius remained consistent withBekki & Chiba (2009).As illustrated in Figure12, there is not a significant difference in the orbital parameters of runaway GCs.Therefore, we conclude that the total mass of a DSG is the most effective parameter in the distribution of its runaway GCs rather than the density profile of the host DSGs.Figure13shows the LMC orbit and two runaway GCs in H1 model.This figure implies that there could be GCs that have escaped from the LMC and are likely to be found at smaller distances from the MW.In the same manner, there could be runaway GCs located at very large distances from the MW.It is clear from Figure11that despite the wide distribution of orbital parameters of runaway GCs, their maximum orbital inclination is less than 38 • .This condition is a strong constraint, leading to only a few MW GCs meeting the criteria for association with the LMC.For the H1 model, NGC 5024, Pyxis, and Pal 3 belong to Flag 1, and Pal 4, NGC 7492, and NGC 5053 belong to Flag 2. For the H2 model, NGC 5024, Pyxis, Pal 3, and Pal 4 belong to Flag 1, and NGC 5053 with NGC 7492 belong to Flag 2. According to Massari et al. (2019) the Pyxis and Pal 3 are considered to belong to high-energy group GCs.The left panel of Figure 14 shows the SMC runaway GCs in the H2 model.GCs are distributed in the range of [50, 180] kpc and  = [0.0,0.4].The maximum orbital inclination of GCs from the SMC orbit is about 10 • .With these conditions, there are not any MW GCs associated with the SMC.This is in contrast to Paper I

Figure 10 .
Figure10.Orbital parameters of the Sgr runaway GCs (dots) and 13 MW GCs in L1 and H1 models, which have been previously suggested to be associated with the Sgr.The area occupied by each MW GC is due to its uncertainty of proper motions and the line-of-sight velocity.The Sgr is also depicted, considering the uncertainty of its observational data.

Figure 11 .
Figure 11.Same as Figure 8 but for the LMC runaway GCs in H1 (without the dynamical friction) and H2 (with the dynamical friction) simulation models.The large squares mark the LMC, and Δ • has been measured with respect to the orbital plane of the LMC.Dynamical friction in H2 model, changes the orbital parameters of the LMC (shown as three squares).As a result, the distribution of the orbital parameters of GCs is more extended compared to H1 model.

Figure 12 .
Figure 12.The orbital parameters of 2500 runaway GCs from the LMC in H1 model, assuming two different density profiles for the LMC, namely Plummer (black) and Hernquist (yellow).

Figure 13 .
Figure 13.The orbit of the LMC (solid black line) as well as two runaway GCs after 8 Gyrs in the H1 model.The orbits are depicted for the  - plane of the Galactocentric coordinate.The solid red line shows an example where the GC is thrown towards the MW with  ap ≈ 30 kpc.The dashed red line is an example where the  ap of the GC reaches 380 kpc.

Figure 14 .
Figure 14.The orbital parameters of 2500 escaped GCs from the SMC in the H2 model without (left panel) and with (right panel) the presence of the LMC.The colour coding indicates the orbital inclination of the GCs with respect to the SMC orbit (in degrees).The orbital parameters of the SMC and the LMC are depicted by large squares.The white-filled squares show the orbital parameters of the SMC at earlier epochs as it changes due to the dynamical friction.

Figure 15 .
Figure 15.The orbital parameters of 2500 GCs escaped from a dissolved DSG in the H3 model after 12 Gyrs.The DSG was placed on the MW disc with an orbital eccentricity of  = 0.2 and an apogalactic distance of  ap = 24 kpc.The DSG was dissolved after 6 Gyrs.The orbital parameters of the DSG at the time of dissolution are shown by the large square.The colour coding indicates the orbital inclination of these GCs with respect to the DSG orbit (in degrees).

Table 2 .
The characteristics of the simulation models.

Table 3 .
List of GCs that could be associated with the Sgr in all simulation models.The numbers indicate the association flags.Flag 1 GCs have the highest probability of being associated with the Sgr (ref.Section 2.2.1).

Table 4 .
The MW GCs that are likely to be associated with a DSG, in H1, H3, L1, and L3 models.Flag 1 and Flag 2 correspond to high and low probabilities of association, respectively.