© The Authors 2026

1 Introduction

At the centres of most galaxies, there lie a supermassive black hole (SMBH) and a massive, dense stellar cluster (Kormendy & Ho 2013; Neumayer et al. 2020). Unlike the SMBHs, the clusters are directly observable, particularly in late-type galaxies (Neumayer et al. 2020). The sizes of these clusters, known as nuclear star clusters (NSCs), are comparable to those of typical globular clusters (GCs), with a half-light radius of approximately 3.5−5.5 pc (Böker et al. 2004; Baumgardt & Vasiliev 2021). However, NSCs significantly surpass GCs in brightness (∼ 10⁷ L_⊙; Matthews & Gallagher 1997; Böker et al. 2002) due to their greater mass (∼ 10⁶−10⁷ M_⊙) and the presence of both young and old stellar populations (Ho et al. 1995; Böker et al. 2002; Lotz et al. 2004) with unique evolutionary features owing to their dense environments, such as stellar collisions and mergers leading to rejuvenated stars, and the depletion of late-type giants due to envelope stripping (Genzel et al. 2003). Certain correlations have been observed between NSCs and their host galaxies (Côté et al. 2006; Wehner & Harris 2006; Ferrarese et al. 2006; Rossa et al. 2006). In particular, a number of studies have noted a correlation between the velocity dispersion of NSCs and that of their host galaxies, as well as between NSC masses and the bulge masses. Some of these relations are also present for SMBHs, raising questions about the mutual role of NSCs and SMBHs in their evolutionary processes.

The morphology of NSCs varies significantly between galaxies and can include both disc-like and spherical components (Seth et al. 2006). Apparently, the morphology of an NSC partially or fully reflects that of its surrounding environment (Capuzzo-Dolcetta & Miocchi 2008a). Understanding the morphology and kinematics of NSCs plays a key role in deciphering the mechanisms of their formation and growth (De Lorenzi et al. 2013). Overall, NSC formation is thought to occur via two main scenarios:

• The accretion of stellar clusters onto galactic centre (often referred to as the dissipationless, cluster-inspiral, or dry-merger scenarios; Tremaine et al. 1975; Lotz et al. 2001; Capuzzo-Dolcetta & Miocchi 2008a,b; Agarwal & Milosavljević 2011; Antonini 2013; Gnedin et al. 2014; Tsatsi et al. 2017; Fahrion et al. 2022; Gao et al. 2024).

• The formation of stars from gas that has previously migrated to the galactic centre (also known as the dissipative, in situ formation, or wet formation scenarios; Loose et al. 1982; Milosavljević 2004; Seth et al. 2006; Nayakshin et al. 2009; Aharon & Perets 2015).

Both scenarios are likely valid (Hartmann et al. 2011; Leigh et al. 2015), with their relative contributions depending on the environment in which the NSC forms. As an example, Fahrion et al. (2022) showed that the dominant NSC formation channel depends on NSC (and host galaxy) mass. The primary question is whether these mechanisms are dominant and (if not) what their role is in the formation and growth of the NSC. Additionally, the mechanisms of NSC formation and subsequent growth may not necessarily be identical. In early-type galaxies, where gas concentration is low, the dominant mechanism for NSC growth is expected to be cluster-inspiral, although in situ formation cannot be excluded during the early stages of such galaxies’ history. The diversity of stellar populations, manifested by the coexistence of different stellar generations and evidence for multiple star-formation episodes in NSCs, also supports the simultaneous operation of both mechanisms. Moreover, as the central region grows in mass, it can more efficiently attract gas, triggering new star formation or enhancing its rate (Schinnerer et al. 2006; Partmann et al. 2025). NSC formation might also result from galaxy interactions and mergers during their dynamic evolution. It is logical to assume that during the accretion of randomly distributed stellar clusters onto the galactic centre would not impart their orbital angular momentum to the forming NSC. However, simulations of randomly distributed stellar cluster accretion have shown that special initial conditions are not required for NSC rotation to develop, although rotation becomes more pronounced if clusters infall with their orientations aligned with the galactic disc’s rotation (Seth et al. 2008).

Unlike NSCs in other galaxies, the NSC of the Milky Way (MW) can be resolved into individual stars due to its relative proximity (∼ 8.1 ± 0.1 kpc; GRAVITY Collaboration 2019; Do et al. 2019). However, interstellar matter along the line of sight to the Galactic centre (GalC) imposes limitations on observing certain NSC properties. The MW’s NSC is thought to have a mass of ∼ 10⁷ M_⊙ (Schödel et al. 2008) and contains the SMBH (Sgr A*) with a mass of ∼ 4.3 × 10⁶ M_⊙ (Gillessen et al. 2009). Observations indicate an increase in star formation rates in the MW’s NSC over the past few hundred million years (Blum et al. 2003; Pfuhl et al. 2011). In the region within ∼ 0.5 pc of the GalC, young stars with ages around 6 Myr dominate (Paumard et al. 2006; Feldmeier-Krause et al. 2015), while the fraction of old stars comparable in age to MW GCs is extremely low (Do et al. 2015). This supports the dominance of the in situ scenario. However, in a ∼ 2.5 pc region around the GalC, 80% of stars are older than 5 Gyr (Blum et al. 2003; Pfuhl et al. 2011). The kinematic features of the MW’s NSC, such as the detection of a dynamically distinct, metal-poor component and the evidence for non-isotropic rotation, favour the cluster-inspiral scenario (Do et al. 2020). The deficit of GCs in the central region compared to the overall stellar distribution in the MW may also be a consequence of GC inspiralling and destruction during NSC formation (Lotz et al. 2001; Capuzzo-Dolcetta & Mastrobuono-Battisti 2009). Additionally, the MW’s NSC exhibits rotation aligned with the Galactic disc’s rotation (Trippe et al. 2008; Schödel et al. 2009). The cluster-inspiral scenario could also partially or fully explain the γ-ray and X-ray excess observed in the GalC as numerous millisecond pulsars, cataclysmic variables, and black holes formed in dense star clusters may have been delivered to the GalC via inspiralling (Arca-Sedda et al. 2018).

To validate the cluster-inspiral mechanism for the formation and growth of NSCs, detailed numerical simulations of stellar cluster dynamics within the potential field of the host galaxy are typically performed, with special attention paid to cluster-GalC and cluster-cluster interactions (Oh et al. 2000; Bekki et al. 2004; Hartmann et al. 2011). A significant factor in these studies is the consideration of dynamics within a live, time-dependent potential, rather than a static one (Bekki 2010a,b).

In this paper, we investigate the cluster-inspiral scenario for the formation and growth of the NSC via the accretion of the MW’s GCs. Overall, GCs are tightly bound stellar conglomerations associated with all types of galaxies. The survivability of GCs during their passages near the GalC depends on their mass and concentration parameter, defined as the logarithm of the ratio of the tidal and core radii, c=log (r_t/r_c). Larger GC mass and higher concentration allow a cluster to survive more passages near the GalC. For instance, Miocchi et al. (2006) simulated the dynamics of GCs in a triaxial system and found a correlation between the mass-loss timescale, τ, and the concentration parameter, τ ≃ 14 c^6.1 t_dyn. Orbital energy loss and structural changes in GCs (especially for the higher masses ≳ 10⁶ M_⊙) can occur due to the dynamical friction and tidal interactions with the galaxy’s potential (tidal dissipation), resulting in GC deceleration and potential accretion onto the GalC, thereby contributing mass (Pesce et al. 1992; Capuzzo-Dolcetta 1993; Capuzzo-Dolcetta & Vicari 2005). Depending on the compactness of the GCs, the number required for NSC formation varies, ranging from several dozen for compact clusters to several hundred for looser ones. Antonini (2013) used simulations to find that the rapid formation of NSCs from GCs depends heavily on the presence or absence of a SMBH at the centre. Perets & Mastrobuono-Battisti (2014) conducted a series of simulations involving the sequential infall of 12 identical GCs (each −1.1 × 10⁶ M_⊙) onto the inner region of the Milky Way, which contains a SMBH (M_•}=4 × 10⁶ M}_⊙), from circular orbits with a radius of 20 pc and random initial parameters. Their results yielded a total NSC mass comparable to observational data: 1.4 × 10⁷ M_⊙.

Our previous studies (Ishchenko et al. 2023a,b) have demonstrated that during their dynamical evolution in a time-varying MW-like external potential from cosmological database IllustrisTNG-100 (Nelson et al. 2019), some MW GCs might closely approach the GalC on a cosmological timescale, potentially contributing to the growth of the NSC’s mass. The frequency of close GC passages near the central region of the Galaxy was estimated and found to be too low to fully account for the NSC’s mass if we are considering only the periodic capture of stars during such encounters. This led us to conclude that it is essential to consider that some GCs may not have survived to the present day, having been tidally disrupted through interactions with the NSC and thereby contributing significantly more to its mass.

The primary objective of this study is to perform detailed dynamic N-body modelling of theoretical GC models within a time-varying Milky Way-like potential, accounting for stellar evolution. This aims to evaluate the contribution of accreted mass from GCs with varying orbital parameters, including cases of complete cluster disruption.

The paper is organised as follows. Section 2 describes the generation of initial conditions for point-mass models of GCs, integrated backward over 10 Gyr to assess close encounters with the GalC. Section 3 details the full N-body integration of theoretical GCs models, incorporating stellar evolution and physical model parameters. Section 4 presents our analysis of stellar accretion onto the proto-NSC, while Section 5 discusses and summarises the results.

Fig. 1 Evolution of halo and disc masses, along with their characteristic scales for 411321 external potential. Green, red, and blue lines with dots show the parameters recovered from the IllustrisTNG-100 data. Black solid lines correspond to the values after the interpolation and smoothing with a 1 Myr time step that was used in the orbital integration. The dashed lines represent the 5 Gyr duration.

2 Mesh of initial condition distributions and GC orbital integration as a point mass

According to our current understanding, it is clear that the initial condition phase space for Galactic GCs progenitors positions and velocities is quite broad, but obviously limited by the Galaxy assembling history in the early Universe (Kravtsov & Gnedin 2005; Pagnini et al. 2023; Ishchenko et al. 2023a). It is also evident that we cannot numerically model the full set of GCs initial conditions phase space distribution. On this basis, we restricted our study of the potential NSC donor GCs with some limits in angular momentum and energy space. As a starting point, for our initial condition generation, we chose the well-studied 411321 TNG-TVP potential (Ishchenko et al. 2023a), as detailed in Fig. 1.

To obtain the spatial scales of the discs and dark matter haloes, we decomposed the mass distribution using the Miyamoto-Nagai (MN) Φ_d(R, z) potential (Miyamoto & Nagai 1975) and Navarro-Frenk-White (NFW) Φ_h(R, z) potential (Navarro et al. 1997), expressed as $$\begin{array}{rl}{\mathrm{\Phi }}_{\mathrm{t}\mathrm{o}\mathrm{t}}& ={\mathrm{\Phi }}_{\mathrm{d}}(R,z)+{\mathrm{\Phi }}_{\mathrm{h}}(R,z)\\ & =-\frac{G{M}_{\mathrm{d}}}{\sqrt{R^2+{(a_{\mathrm{d}}+\sqrt{z^2+b_{\mathrm{d}}^2})}^2}}-\frac{G{M}_{\mathrm{h}}\cdot \ln (1+\frac{\sqrt{R^2+z^2}}{b_{\mathrm{h}}})}{\sqrt{R^2+z^2}},\end{array}$$(1) where $R=\sqrt{x^2+y^2}$ is the planar galactocentric radius, z is the distance above the plane of the disc, G is the gravitational constant, a_d is the disc scale length, b_d, h are the disc and halo scale heights, respectively, and M_d and $M_{\mathrm{h}}=4\pi {\rho }_0{b}_{\mathrm{h}}^3$ (ρ₀ is the central mass density of the halo) are the masses of the disc and halo, respectively.

In this potential, we step back to the 10 Gyr look-back time and put our theoretical GC’s centres in the Galactocentric Cartesian coordinate frame. In this potential for the selected time, we see that the minimum specific energy level is around – 18.6 × 10⁴ km² s⁻². Thus, in the generation routine, we started our energies from −19 to −14 × 10⁴ km² s⁻². For the specific angular momentum space, we chose the limits: 0−3 × 10² kpc km s⁻¹ for L_tot to cover a full range of orbital ecc almost from 0 to 1 and from −0.5 to +0.5 × 10² kpc km s⁻¹ for L_z to cover a significant range of all prograde and retrograde orbits of the currently observed and modelled Galactic GCs (Ishchenko et al. 2023b, 2024).

Using the simple rejection method, inside this mesh of energy and angular momentum, we generated limits randomly with a uniform distribution of more than 16 k initial models. For each of these points in a phase space, E, L_tot, L_perp, L_z, we find a corresponding random position and velocity for our theoretical GC centres in the Galactocentric Cartesian coordinate frame. In Fig. 2 (top), we present the initial phase space distribution of GCs in a specific energy and angular momentum space. We are aware of the limitations of our approach for the theoretical model GC centres’ phase space distribution, which is connected with the underlying assumption of the similarities between currently observed and currently already tidally destroyed GCs. However, we use this assumption here as a first step to generate our set of theoretical GCs. In the future, we will also investigate other possibilities.

Using this set of initial positions and velocities, we ran our orbital integration for 5 Gyr, from −10 Gyr look-back time. For this task, we used the high-order parallel dynamical N-body code ϕ-GPU¹ (Berczik et al. 2011, 2013). The code is based on the fourth-order Hermite integration scheme with hierarchical individual block time steps. In Fig. 2, we present mean values for orbital ecc via colour palette, calculated during 5 Gyr of dynamical integration time. The orbital eccentricity we calculated using the simple estimation based on the apo- and peri-centres of the particle orbit: ecc =(r_apo−r_per) /(r_apo+r_per). In Fig. A.1 as an illustration, we present orbits for each representatives of the ecc, from 0.0 to 0.9.

In the next step, we analysed all GC models to see the number of close passages with relative distance less than 100 pc near the GalC during the full 5 Gyr of evolution. We chose this indicator as the first pre-selected option for future GCs list of potential NSC contributors, which will be simulated with full N-body particle systems. From the initial 16 k GCs models, we found ∼ 12 k models that have recorded at least one close pass. In Fig. 2 (bottom panels) via colour palette, we present the number of close passages. As we see from these plots, high numbers of close passages have models, which are distributed in low specific energy limits, from −17 to −18.5 (bright yellow dots). Models with a low probability of close passes are of little interest in our case of NSC formation via GC accretion (dark dots). In the next step, we prepare a set of models, which have more than 100 events of close passage near the GalC per model. We found 6621 such models.

To reduce our number of models for future N-body simulations, we assumed the second pre-selection option, namely, ecc. In this way, we selected five GC models, which are representatives for each ecc bin in steps 0.1, from 0.0 to 0.9 (single blue dots in Fig. 2 in bottom panels). We note that for ecc 0.9, there are no models that have at least one close passage. However, to maintain consistency, we have selected some of these models. In total, based on our two pre-selection options (described above), we selected 50 GC models for full N-body integration. In Fig. 2, we present the initial distribution for these selected models in the phase space E versus L_tot, E L_perp, and E versus L_z as blue dots. In physical space, these GC objects are in a range of distance from the GalC of 100 pc−5 kpc.

Fig. 2 Distribution of the GCs in phase space: E vs Ltot, and E vs Lz in 411321 external potential. Top: initial phase space distribution, where the colour-coding shows orbital mean ecc, based on 5 Gyr of integration. Bottom: interactions (close passages) with the GalC, based on 5 Gyr of integration. Colour-coding represents the number of such events. Black dots represent selected GCs for future N-body simulations.

3 Full N-body integration of GC’s

To carry out the N-body dynamical modelling together with the stellar evolution for selected GCs, we assumed the initial conditions described below. The GCs centres initial coordinates and velocities at the 10 Gyr lookback time were taken from the orbital integration of our theoretical GC from the previous step (see Sect. 2).

For our N-body modelling of the theoretical GCs, we created three sets of models with different half-mass radii: 1, 2, and 4 pc. Other physical parameters, such as initial masses and King concentration parameters, we assumed the same for all models: 60 × 10³ M_⊙ and W₀=8.0. We generated the initial positions and velocities of the stars inside the clusters using the nowadays very popular Agama library (Vasiliev 2019). For the individual initial mass of the stars, we used the Kroupa mass function (Kroupa 2001) with lower-upper mass limits equal to 0.08−100 M_⊙. Each cluster was initially set in a dynamical equilibrium state with the distribution function of the theoretical King model (King 1966). Because our GC models have relatively low mass, for survivability we try to compensate for this with the higher King concentration parameter for all the model systems. In Table 1 we summarise these parameters. In Fig. 3, we show initial density distribution for our models with different initial halfmass radii.

In Fig. 4, we present the cumulative distribution of initial star masses for different r_hm models with the selected Kroupa (2001) initial mass function. As we can see, all the lines are fully consistent. To carry out the full N-body modelling for our three sets (in total 150 GC) of models, we used the same N-body code ϕ-GPU, which was mentioned in Sect. 2, but with the up-to-date stellar evolution prescription (Kamlah et al. 2022).

Based on the initial mass function and total particle numbers described above, for each GC model, we expect to have 174 black holes (BHs), 723 neutron stars (NSs), and ∼ 8000 white dwarfs (WDs); however, the majority (∼ 87 k) will still be made up of low-mass stars on the main sequence after 5 Gyr of evolution. As an initial model for the Galactic NSC, we included a special central proto-NSC object with a mass of 4 × 10⁶ M_⊙ and a radius of 10 pc in our simulations. This object acts as an extra-gravitational term, which we include to our global Galactic potential 411324 TVP. During the N-body simulation, we assumed that a star from the GC has accreted to the NSC when the particular GC star is inside the 10 pc radius of the proto-NSC. When such a star was detected, we added the current mass of the star at the time of the pass to the NSC and considered the star to be accreted. In such a case, we will always have a growing NSC mass. Such a simple star versus proto-NSC merging scenario (based only on distance) might overestimate the number of accretion cases with some high-velocity objects. To overcome this effect, we restricted our study to GC systems only with the orbital velocities less than 150 km/s, which lies in a range of the bounding velocity around proto-NSC with the masses and sizes considered above.

Fig. 3 Initial density distribution for our models with different initial half-mass radii.

Fig. 4 Normalised cumulative distribution for initial star masses.

4 Analysis of the GCs stellar accretion processes on to proto-NSC

In this section, we discuss in general terms the stellar accretion rate onto the NSC during the full duration of simulations (5 Gyr) for our three sets of models. We illustrate the accretion rate in the context of the (i) orbital ecc, (ii) initial phase space conditions, and (iii) contribution from the GC’s high-mass stellar remnants.

4.1 Global accretion rate onto proto-NSC

In Fig. 5, we show the dynamics of stellar accretion rates onto NSC for all 50 GCs models with three different values of r_hm =1, 2, and 4 pc. The data are presented as a 2D histogram of simulation time versus cumulative accreted mass onto NSC. To enhance the analysis and reveal trends and regions with the highest contributions, the histogram bins were smoothed using a Gaussian kernel. The plot was constructed using a 2D histogram with 100 × 100 bins, interpolated over five points, and smoothed using a Gaussian kernel with σ=0.2. The used Gaussian kernel affects only the sharpness of the contours, not the overall distributions on the plots.

Figure 5 clearly illustrates the difference in the timing of major accretion episodes onto the proto-NSC as a function of the model cluster’s density (parametrized by half-mass). As shown, for clusters with r_hm=1 pc, accretion events occur throughout the integration time interval. Models with r_hm=2 pc have a maximum accretion rate during the first 1 Gyr and later around 2.5−3 Gyr. The last set of models with r_hm=4 pc also includes an early accretion phase that almost ended after 4 Gyr. For all three sets of models, the proto-NSC accretion process has a clear maximum set around 0.5 Gyr.

To illustrate the accretion process onto the central NSC, in Fig. 6, we present the polar and azimuthal angles of particles at the moment of accretion in the Cartesian Galactocentric coordinate system, assuming a proto-NSC accretion radius of 10 pc. For the presentation, we chose three models with high accretion rates, but with different cluster lifetimes. Two models (No 2065 and 155 in Table B.1) have a short dissolution time of around 0.6−0.9 Gyr. The third model (No 1371) is one of our long-lifetime clusters, but they expel almost all the stars during the 5 Gyr evolution time. The short-lifetime clusters show a strong asymmetric and patchy distribution of the accretion angle. They both have a strong accretion event during the first few hundred megayears (Myr; large blue regions). In contrast, the long-lifetime cluster has a more regular and uniform accretion angle distribution. The accretion events (each dots in these figures are individual accretion moments) are distributed equally in time (from blue to red). Besides the patchy structures in some of these plots, we can generally conclude that the accretion processes on our central 10 pc sphere are quite random and stochastic.

In Table 2, we present the total number of accreted stars and their total mass for each of all 50 GC models for our three sets of GCs. As can be seen, the model with r_hm=1 pc contributed more accreted masses compared to the other models. Based on the values presented in the table, we can conclude that our extended set of 150 numerical simulations provides an average lower limit of the mass contribution to the NSC from all investigated GCs of ≈ 6%.

In Table B.1, we present detailed statistical stellar accretion information for all 50 GCs for each of our three sets. As can be seen, at least 5 GCs are completely dissolutes during 5 Gyr of simulations. We define the dissolution of a cluster as occurring when the remaining tidal mass decreases to less than 5% of the cluster initial mass. With initial r_hm =1 pc we have five dissolved GCs during this 5 Gyr integration period. With an initial r_hm = 2 pc, we have 9 GCs, and with r_hm=4 pc−42, respectively. In most cases, the dissolution of the GCs occurs in the models with the orbital ecc in a range of 0.3–0.6.

Figure 7 illustrates the contribution of GCs stars to the mass of the NSC over a period of 5 Gyr for the case of r_hm=2 pc. The different colours represent the mean orbital eccentricity, as indicated in the legend. The mass is expressed in solar masses M_⊙ and is shown on a logarithmic scale. To better resolve the early accretion phase, we have intentionally divided the first Gyr into two intervals. This is because the highest accretion rates occur within this period for all models. This might indicate that NSC formation, at least in the cluster-inspiral scenario, occurs predominantly in the early stages of Galaxy assembling. The figure also shows that models with mean ecc of 0.4 and 0.5 make the largest contribution in mass of the NSC. This may be attributed to the fact that GC models with these mean eccentricities are also characterized by pronounced orbital eccentricity oscillations over our 5 Gyr timescale, spanning a broad range of ecc values (also see Fig. 9).

In Fig. 8, we present the mass accretion rate to the NSC as a function of specific energy and angular momentums: L_tot, L_perp, and L_z at the initial moment of time (i.e. 10 Gyr ago in 411321 external TVP potential). As we see from the plots, the highest mass accretion rate for GC models is a system, in which the initial phase space distribution has a nearly perpendicular orbital plane to the galactic disc, namely, of ∼ 0.0 in L_z (yellow-red range in palette colours).

Fig. 5 Stellar accretion rates onto the NSC for 50 GCs models for three different values of rhm = 1, 2, and 4 pc, respectively.

Fig. 6 Polar and azimuthal angles of particles at the moment of accretion in the Cartesian Galactocentric coordinate system. The accretion radius sphere of the proto-NSC was set to 10 pc. The points individual colour coding shows the moment of accretion. Each panel corresponds to the accretion from different GC.

Fig. 7 Contribution of GCs stars to the mass of the NSC over 5 Gyr for different orbital eccentricities for rhm=2 pc case. The colours represent the orbital eccentricity, as shown in the legend. The mass is given in solar masses M⊙ and is displayed on a logarithmic scale.

Fig. 8 Initial distribution for 50 theoretical GCs in E- Ltot, E−Lz at −10 look-back time. Different symbols show the rhm values. Colour represent accreted mass from each GC to the NSC over 5 Gyr of integration time.

Fig. 9 Evolution of the mass accretion to the NSC for models with rhm depends on the ecc and inclination angle α. Short black lines: mean values, colour box − ± 1 σ for ecc, thin colour lines: minimum and maximum ecc, and α values during 5 Gyr of integration. Colour: accreted mass to the NSC.

4.2 Mass accretion to the NSC as function of GC orbital parameters distribution

After running the full set of models, we analysed the orbital parameters of each individual run and checked the correlation between the model’s orbital parameters and the amount of mass, accreted onto proto-NSC from each model. As main orbital parameters, we chose the orbits ecc and inclination angle (defined as cos (α)=L_z/L_tot). We chose these parameters because of the strong dependence of the accreted stellar mass on these values. Due to the time-variable nature of our external Galactic potential in each orbital revelation, we derived the ecc and inclination angle, α. We also prepared the statistics among each orbital revelation and derive the average and the standard deviation for these orbital parameters. As we can already see from Fig. A.1, the time evolution of the orbital shapes is quite complex and, as a consequence, we have quite significant standard deviations, especially with respect to ecc.

In Fig. 9, we show the models’ mass accretion onto NSC (colour coded) as a function of the model’s initial ecc and inclination angle, α. Each model has been numbered (X-axis) according to its average ecc. As we can see, the models with the highest mass accretion (model numbers from 10 to 37) also have the highest ecc changes during the simulation. The most prominent GC mass donors have a ecc change from 0.1 to 0.9, represented in red colour.

One particularly notable feature of the accretion process in our set of models is the highest accretion for the runs with close to 90 degree of orbital inclination. These models (from 23 to 35) have almost exactly a 90 degree with small deviation (thick colour box). It is worth mentioning that these models have also a quite large of minimum-maximum inclination angles, during the time evolution (thin colour lines). It is also interesting that angle range have a quite strong asymmetry, namely, the indication to the significant counter rotation of these runs, from 90–180 degrees.

For more deep illustration of the dynamical evolution for orbital parameters, such as ecc, α, and semi-major axis, we present several plots for the models with different accretion rates onto the proto-NSC (numbers in right panels) in Fig. C.1. Here, in the middle panels, we see a certain compression of the semimajor axis of the orbit towards the GalC around 1.5−4 Gyr. The dynamics of the orbits is directly related to the time evolution of the Galactic potential, namely, to the halo mass and scale length for the same time period. For GC models that have evolution in the Galactic Pale, the Myamoto-Nagai scales a_d and b_d also have an influence, squeezing the scales for the disc (see Fig. 9). The left panels show the orbital ecc evolution of the individual models. As we can see, the most mass-donor GCs have a wide range of ecc changes, ranging from 0.1 to 0.9.

Based on Figs. 9 and C.1, we can conclude that same most mass donor models have close to constant orbital inclination angle of almost 90 degrees. A combination of ecc and inclination angle α gives us a best candidates for the GC models, which can be a source of the stellar mass accretion onto proto-NSC.

Fig. 10 Total accretion of the mass of stellar remnants (WD, NS, and BH) from 50 GCs for all the three cases: rhm=1, 2, and 4 pc. The mass in given in logarithmic scale.

4.3 Accretion rate for the stellar remnants from GCs

We analysed the accretion of stellar remnants (WDs, NSs, and BHs) from 50 GCs for three different values of r_hm=1, 2, and 4 pc. Table 3 shows the total number and masses of stellar remnants accreted onto the NSC for these three sets of models.

Since the r_hm=1 pc case has the highest amount of total accreted mass between three sets of models, it is not surprising that it also contributes the most mass to stellar remnants. Figure 10 shows the total mass of accreted stellar remnants. For NSs and BHs, the mass contributions are comparable among all three models, reflecting their low absolute numbers in the stellar population. In contrast, the significantly higher mass of WDs accreted in the r_hm=1 pc model likely results from their higher intrinsic abundance, as WDs are the dominant stellar remnants in clusters with the standard type of initial mass functions. The enhanced central density of the stellar remnants in such a compact cluster further amplifies their contribution.

In Appendix D, we present a more detailed analysis of the accretion of stellar remnants onto the NSC. As expected, for all three cases of r_hm=1, 2, and 4 pc, the largest contribution to the NSC mass over 5 Gyr comes from stellar remnants accreted from GCs with mean orbital eccentricities around 0.4–0.5. In terms of the temporal evolution, the highest accretion rate for all types of remnants occurs in the interval between 0.5 and 1 Gyr. The analysis of individual accretion events shows that the bulk of the accreted objects are WDs and low-mass stars.

Fig. 11 Global tidal mass loss and accretion mass maps over all our theoretical GC models. The X and Y axes are the orbital parameters of our GC’s. The Z axis is a tidal mass of the clusters after 5 Gyr of dynamical integration. In the bottom row, the colour-coding corresponds to the accreted mass onto proto-NSC from the individual clusters.

5 Discussions and conclusions

In this section, we discuss the possible effects that can lead to the stellar accretion of GCs onto the proto-NSC during the long-term cosmological time scale dynamical evolution of the theoretical GCs system. In more detail, we discuss (i) the role of the GC mass loss itself and the possible accretion onto the GalC, (ii) the role of stellar remnants, and (iii) the physical effects of the time-variable potential on the accretion process.

5.1 Role of the GC’s mass loss and accretion rate for low- and high-mass cases

As a basic GC model in our set of simulations, we choose systems with an initial mass of 60 × 10³ M_⊙. The general picture of the cluster’s global mass loss for three different initial half-mass radii is presented in Fig. 11. The top row presents the global mass loss of clusters with different orbital parameters, while the bottom row (colour-coded) shows the proto-NSC-accreted mass from the GCs. A clear effect of the half-mass radius is visible: clusters with larger half-mass radii experience significantly greater total tidal mass loss across all orbital parameters (see the ‘blue plateau’ in the top-right panel). In general, the low-eccentricity models for all the different half-mass radii lose significantly less mass compared to the models with the higher ecc.

The models with r_hm=4 pc are almost completely dissolved; however, these dissolved clusters do not contribute significantly to the mass accreted onto the proto-NSC. In other words, the global tidal mass loss of some GC models does not necessarily translate into substantial stellar accretion onto the central region of the Galaxy. As mentioned in Sect. 3, we did not include any dynamical friction to the prescription of the accretion process. The possible dynamical friction effects can change the actual values of the accreted stellar masses (in the direction to the higher masses), but our general conclusions regarding the dependency of the NSC accreted mass from models half-mass radii, will be qualitatively valid even in this case.

In addition to the main set of GC models, each with a a total initial mass of 60 × 10³ M_⊙, we additionally simulated the dynamical evolution of four models from the total set with r_hm = 2 pc and ecc values from 0.3 to 0.6. For these models, higher initial masses of 120 × 10³ M_⊙ and 180 × 10³ M_⊙ were adopted (see Table 1). The general mass loss of these additional models is presented in Fig. 12. As expected, the lower-mass models (among the set with the same initial half-mass systems) have a slightly larger general mass loss ≳ 90%, compared to the higher-mass models. The models with 120 × 10³ and 180 × 10³ M_⊙ exhibit very similar global mass loss behaviour, with values around ∼ 60−80%.

The main results of the stellar accretion simulations for the more massive GC models are presented in Table 4, where the relative fractions of the accreted mass and the number of particles are given as percentages of their initial values (also see Fig. 13). As the initial mass of the GC models increases, their total mass contribution to the proto-NSC also grows. However, for the model with ecc of 0.4, the relative fraction of accreted stars decreases when the mass increases from 60 × 10³ to 120 × 10³ M_⊙, although at 180 × 10³ M_⊙ the model no longer deviates from the established trend. It is worth noting that although the model with ecc = 0.6 yields only a small number of stars, their average mass is approximately 4.4 M_⊙, compared to ≤ 3 M ⊙ for the other models.

As a side remark, we also ought to mention the generally low average mass of the accreted stars. For models with ecc 0.4 and 0.5, where the accretion statistics are more robust, the average masses of the accreted stars are ≈ 0.33 M_⊙. Compared with the average IMF mass of our clusters (≈ 0.57 M_⊙), we see that during our 5 Gyr dynamical simulations, the accreted population is primarily dominated by low-mass stars. One possible reason for this mass selection is a significant mass segregation among the clusters, resulting in a higher fraction of high-mass stars and remnants residing in the more tightly bound central region. As an example, we checked the average cumulative masses in our typical theoretical GCs (2858 with r_hm=1 pc) over the 5 Gyr of simulations and we find that after ∼ 1 Gyr the average cumulative mass within 0.1 pc is almost around ∼ 20 M_⊙ and it quickly drops below the level of ∼ 0.3 M_⊙ after 10 pc.

Fig. 12 Evolution of the mass loss in the GC models due to mass loss and orbital type after 5 Gyr of integration. Values Mini are in M⊙.

Fig. 13 Relative accretion from low- and high-mass GC models with different ecc as percentages.

5.2 Role of the stellar remnants

Our analysis of stellar remnants accretion onto the NSC from GC models reveals a strong dependence of this process on the initial half-mass radius, r_hm, and orbital eccentricities, ecc (see Fig. 10). The time-resolved accretion history shows a pronounced early peak within the first Gyr (see Fig. D.2).

As seen from Fig. 10, the largest total contribution of remnants to the accreted mass comes from more compact GC models (r_hm=1 pc, violet colour). It is worth noting that orbital eccentricity also plays a key role in the process of accretion efficiency. The highest rates remnant accretion are observed at intermediate values of ecc = 0.4−0.6 (see Fig. D.1). The accretion models (Figs. 10 and D.2) also show a clear distinction among the remnant types in terms of both timing and mass contribution, with WDs being the dominant component across all scenarios.

We also note the significantly higher accreted mass of BHs compared to that NS remnants in the global accreted stellar remnant mix (by almost a factor of 10; see Fig. 10). However, according to the adopted cluster IMF, the number of forming BHs is about four times smaller than the number of NS remnants. This apparent contradiction can be explained by our selective supernovae kick mechanisms (Kamlah et al. 2022), where NS remnants always receive a significant velocity kick; meanwhile a substantial fraction of BHs, due to the fallback prescription, have a zero velocity kick at the moment of formation. These results highlight the relationship between NSC-GC interactions and the presence of stellar remnants in the NSC during its formation and growth.

5.3 Role of the external potential

As we see from our study, the actual accretion onto the central galactic-scale proto-NSC is a complex dynamical process. In our picture of accretion onto the proto-NSC, we modelled the accretion (dynamical capture) of individual stars onto a forming the NSC. The detailed statistical picture of the dynamical accretion process is presented in Fig. 5. In Sect. 4.1, we describe the general accretion trends for the different sets of models. Here, we would like to discuss the direct influence of the modelled time-variable Galactic potential on the growth of the NSC. The mass and size evolution of our selected potential 411321 is shown in Fig. 1.

The dark matter halo mass and characteristic scale exhibit a fairly monotonic growth over time. Since the halo characteristic scale is quite large (15−25 kpc in the selected integration interval from −10 to −5 Gyr), the direct influence of the extended halo potential is not significant for our set of the theoretical GC models (the maximum apocentre is 5 kpc). In contrast, the baryonic disc potential shows a markedly different behaviour. The ratio of a_d and b_d changes significantly over the entire integration interval, from ∼ 1 to 3. Such behaviour of the inner Galactic potential has a strong influence on our theoretical set of GCs.

In Fig. C.1, we see the direct effect of these changes on the evolution of the semi-major axis. Since these changes occur around 2–3 Gyr of integration time and we observe significant additional mass accretion in that same timeframe (see Fig. 5), we can directly connect these changes.

5.4 General remarks

Based on the total accreted mass from all three sets of our theoretical 150 GCs during the first 5 Gyr of Galactic evolution, the forming NSC acquires ∼ 4 × 10⁵ M_⊙, see Table 2. This mass corresponds to roughly 2% of the current NSC mass of ∼ 2 × 10⁷ M_⊙ (Bland-Hawthorn & Gerhard 2016). If we also include the set of models with higher initial masses, 120 and 180 × 10³ M_⊙, in our global accretion estimate and assume the same frequency distribution of such clusters as for the lower-mass GC models (60 × 10³. M_⊙), then we can expect at least several times more stellar mass to be accreted onto the forming NSC. So, in the end, according to our set of three relatively low-mass GC models, we can gain around ≈ 6% of the current NSC mass via pure stellar accretion from the theoretical GC models (very conservative, lower limit mass estimation).

At first glance, our findings seem to be in contradiction with the observationally motivated NSC accretion estimates (around 50%) reported by Fahrion et al. (2022) for the MW NSC with a mass ∼ 10⁷ M_⊙. However, we need to take into account the fact that the initial GC masses adopted in our models of theoretical clusters are quite small. In this case, we can estimate indeed only a lower boundary of the NSC accretion mass. In the future, we will employ significantly higher initial GC mass models or a mixed population of low- and high-mass clusters and, in addition, we probably will be able to once again achieve higher larger accretion rates by a factor of few. In addition, in our next set of runs, we will include the simplified Chandrasekhar type of dynamical friction (Chandrasekhar 1943a,b) with some physically motivated theoretical parameters (Just et al. 2011).

In addition to the current set of relatively low-mass GC N-body models, which do not include possible dynamical friction, our estimate of the accreted mass represents a lower limit due to another simplification: we omitted the bulge component in the external time-variable potential model, which is currently based on the Illustris-TNG100 cosmological database.

Also, we need to take into account that probably none of the two accretion scenarios mentioned in the introduction (i.e., stellar or/and gas accretion followed by in situ star formation) are likely to operate in isolation. Therefore, in reality, we are dealing with a same combination of both processes (Fahrion et al. 2022). In our current scenario (low-mass GC destruction and star accretion), the remaining 80−90% of the NSC mass is expected to form during a phase of significant gas accretion onto the NSC.

This view is consistent with the high gas fraction in the GalC during the first few Gyr of Galactic evolution.

Acknowledgements

The authors thank the anonymous referee for a very constructive report and suggestions that helped significantly improve the quality of the manuscript. This work was funded by the Science Committee of the Ministry of Education and Science of the Republic of Kazakhstan (Grant no. AP22787256). PB, MI, OV and OS thank 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.

References

Appendix A GCs orbital evolution for different ecc

In these plots we show GCs orbital evolution during 5 Gyr of forward integration in 411321 external time-variable potential. Each GC was integrated as one physical particle. We present only 11 cases from 50 GCs to demonstrate the evolution of orbital ecc from 0.0 to 0.9. We discuss grey colour models (with dynamical friction) in App. E.

Fig. A.1 Orbital evolution during 5 Gyr of integration, from −10 Gyr, represents orbits with different ecc, from 0.0 to 0.9. Orbits with dynamical friction (grey colour) we discuss in App. E.

Appendix B Detailed statistic for GC models

Table B.1 presents the mean orbital parameters – semi-major axis (a), eccentricity (ecc), and pericentre (r_per) – of all our 50 GC models with M_ini=60 × 10³ M_⊙, obtained from their preliminary backward integration (to −10 Gyr) as point objects, as well as the results of subsequent integration (to −5 Gyr) of theoretical full N-body GC models with different half-mass radii r_hm=1, 2, and 4 pc. As can be seen, the largest contribution to the proto-NSC comes from models with ecc values of 0.4 and 0.5. It can also be noted that the more diffuse the cluster, the higher the probability of its complete dissolution (ds). For example, at r_hm=4 pc, more than 80% of the GC models are fully dissoluted (we consider a GC dissoluted if less than 5% of its initial mass remains).

Appendix C Evolution of the ecc, a, and inclination angle

Fig. C.1 shows the evolution of orbital parameters (eccentricity, semi-major axis, inclination angle) over 5 Gyr for representative GC models that contributed less than 10 M_⊙, between 10 and 100 M_⊙, and more than 100 M_⊙. As shown, models whose eccentricities and semi-major axes oscillated over a broad range of values made a larger contribution to the accreted mass. However, this trend does not extend to the inclination angle.

Fig. C.1 Evolution of the orbital parameters ecc, a, and α over 5 Gyr of GCs simulation. First two plots represent models with accretion less then 10 M⊙, two middle panels – 10−100 M⊙ and two bottom – more than 100 M⊙.

Appendix D Additional analysis of the accretion rate for the stellar remnants from GCs

In Fig. D.1, we present the distribution of the NSC accreted masses as a function of the GCs initial orbital ecc. At all evolutionary stages, the largest contribution to accretion comes from GC models with mean orbital eccentricities of 0.4 and 0.5. These trends are clear for all three cases of r_hm=1, 2, and 4 pc. Interestingly, for white dwarfs in all three cases, the accreted mass drops sharply beyond the peak at eccentricities 0.6 and 0.7, then it increases sharply at eccentricity 0.8. In contrast, on the lower-eccentricity side, the mass decreases more gradually, without a pronounced dip.

Fig. D.1 Masses of accreted stellar remnants (WDs, NSs, and BHs, indicated by colours) from GCs models as a function of the orbital eccentricity of their host clusters, shown for rhm=1, 2, and 4 pc (from left to right).

In Fig. D.2 we illustrate the contribution of individual types of stellar remnant to the NSC mass over 5 Gyr (in different time bins), based on the 50 GC models with initial r_hm=2 pc. The results are shown for different mean orbital eccentricities. The highest accretion rate for all stellar remnants occurs in the time interval between 0.5 and 1 Gyr. Then it gradually declines, with a smoother decrease in the case of WDs and a more pronounced drop for NSs and BHs.

Fig. D.2 Contribution of individual types of stellar remnants to the NSC mass (from left to right: WDs, NSs, and BHs) from all 50 GC models over 5 Gyr for different orbital eccentricities in the case of rhm=2 pc.

Fig. D.3 presents the dynamics of individual stellar accretion events onto the NSC for 50 GCs in the case with r_hm=2 pc. For models with initial r_hm=1 and 4 pc, we have similar accretion plots. The data are shown as a 2D histogram of time versus the mass of the accreted star, with a bin size of 100 for both axes. The left panel highlights the moments of accretion only of the stellar remnants, while the right panel includes all stars. Based on the stellar mass, we can clearly distinguish the relative frequencies and accretion timings of black holes, neutron stars, and white dwarfs. The left panel in Fig. D.3 clearly illustrates the boundary separating white dwarfs at the lower end from neutron stars and black holes at higher masses. The right panel demonstrates that the dominant contribution arises from low-mass stars with masses up to approximately 1.4 solar masses, corresponding to the Chandrasekhar limit.

Appendix E Influence of the dynamical friction on the accretion rate to the proto-NSC

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 expressed as (Binney & Tremaine 2008): $$\frac{\mathrm{d}{\mathbf{V}}_{\mathrm{GC}}}{\mathrm{dt}}=-\frac{4\pi G^2\rho M_{\mathrm{GC}}}{V_{\mathrm{GC}}^3}\chi ·ln\Lambda ·{\mathbf{V}}_{\mathrm{GC}}\text{with}\chi =\frac{\rho (<V_{\mathrm{GC}})}{\rho }.$$(E.1)

Fig. D.3 Stellar remnants accretion rate onto the NSC for 50 GCs models for rhm=2 pc. Left: the mass-time distribution for stellar remnants (WDs, NSs, and BHs); right: for all stars. From the left graph, the boundary between white dwarfs at the bottom and neutron stars and black holes at the top is clearly visible (compact object mass gap). It is also clearly seen in the right figure that the main accretion number contribution comes from low-mass stars, approximately up to ∼ 1.4 solar masses (according to the Chandrasekhar limit).

In general, the functions, χ, and the Coulomb logarithm, Λ, depend on the velocity of the massive object and on 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 calculations we use a fixed Coulomb logarithm value ln Λ=5 and a fixed value for χ=0.5.

In Fig. A.1 we present the comparison of the evolution of the GC orbits over 5 Gyr of integration time for models including dynamical friction (grey colour) and without it (coloured lines). As can be seen from the plot, dynamical friction only slightly affects the orbits. In our opinion, applying dynamical friction in its present form does not significantly change the GC orbits, and so the position of lost stars from these GCs on the phase-space diagrams.

Also we would like to mention that one of the key points in our investigation is a quite strong dependence of mass accretion from clusters onto the proto-NSC for orbits with the high inclination angles, see Sect. 4.2 and Fig. C.1. For such high-Z orbits, the effect of dynamical friction is also much weaker compared to models where the GC orbits lie inside the Galactic disk.

All Tables

All Figures

Fig. 1 Evolution of halo and disc masses, along with their characteristic scales for 411321 external potential. Green, red, and blue lines with dots show the parameters recovered from the IllustrisTNG-100 data. Black solid lines correspond to the values after the interpolation and smoothing with a 1 Myr time step that was used in the orbital integration. The dashed lines represent the 5 Gyr duration.

In the text

Fig. 2 Distribution of the GCs in phase space: E vs Ltot, and E vs Lz in 411321 external potential. Top: initial phase space distribution, where the colour-coding shows orbital mean ecc, based on 5 Gyr of integration. Bottom: interactions (close passages) with the GalC, based on 5 Gyr of integration. Colour-coding represents the number of such events. Black dots represent selected GCs for future N-body simulations.

In the text

	Fig. 3 Initial density distribution for our models with different initial half-mass radii.
In the text

	Fig. 4 Normalised cumulative distribution for initial star masses.
In the text

	Fig. 5 Stellar accretion rates onto the NSC for 50 GCs models for three different values of rhm = 1, 2, and 4 pc, respectively.
In the text

	Fig. 6 Polar and azimuthal angles of particles at the moment of accretion in the Cartesian Galactocentric coordinate system. The accretion radius sphere of the proto-NSC was set to 10 pc. The points individual colour coding shows the moment of accretion. Each panel corresponds to the accretion from different GC.
In the text

	Fig. 7 Contribution of GCs stars to the mass of the NSC over 5 Gyr for different orbital eccentricities for rhm=2 pc case. The colours represent the orbital eccentricity, as shown in the legend. The mass is given in solar masses M⊙ and is displayed on a logarithmic scale.
In the text

	Fig. 8 Initial distribution for 50 theoretical GCs in E- Ltot, E−Lz at −10 look-back time. Different symbols show the rhm values. Colour represent accreted mass from each GC to the NSC over 5 Gyr of integration time.
In the text

	Fig. 9 Evolution of the mass accretion to the NSC for models with rhm depends on the ecc and inclination angle α. Short black lines: mean values, colour box − ± 1 σ for ecc, thin colour lines: minimum and maximum ecc, and α values during 5 Gyr of integration. Colour: accreted mass to the NSC.
In the text

	Fig. 10 Total accretion of the mass of stellar remnants (WD, NS, and BH) from 50 GCs for all the three cases: rhm=1, 2, and 4 pc. The mass in given in logarithmic scale.
In the text

	Fig. 11 Global tidal mass loss and accretion mass maps over all our theoretical GC models. The X and Y axes are the orbital parameters of our GC’s. The Z axis is a tidal mass of the clusters after 5 Gyr of dynamical integration. In the bottom row, the colour-coding corresponds to the accreted mass onto proto-NSC from the individual clusters.
In the text

	Fig. 12 Evolution of the mass loss in the GC models due to mass loss and orbital type after 5 Gyr of integration. Values Mini are in M⊙.
In the text

	Fig. 13 Relative accretion from low- and high-mass GC models with different ecc as percentages.
In the text

	Fig. A.1 Orbital evolution during 5 Gyr of integration, from −10 Gyr, represents orbits with different ecc, from 0.0 to 0.9. Orbits with dynamical friction (grey colour) we discuss in App. E.
In the text

	Fig. C.1 Evolution of the orbital parameters ecc, a, and α over 5 Gyr of GCs simulation. First two plots represent models with accretion less then 10 M⊙, two middle panels – 10−100 M⊙ and two bottom – more than 100 M⊙.
In the text

	Fig. D.1 Masses of accreted stellar remnants (WDs, NSs, and BHs, indicated by colours) from GCs models as a function of the orbital eccentricity of their host clusters, shown for rhm=1, 2, and 4 pc (from left to right).
In the text

	Fig. D.2 Contribution of individual types of stellar remnants to the NSC mass (from left to right: WDs, NSs, and BHs) from all 50 GC models over 5 Gyr for different orbital eccentricities in the case of rhm=2 pc.
In the text

Fig. D.3 Stellar remnants accretion rate onto the NSC for 50 GCs models for rhm=2 pc. Left: the mass-time distribution for stellar remnants (WDs, NSs, and BHs); right: for all stars. From the left graph, the boundary between white dwarfs at the bottom and neutron stars and black holes at the top is clearly visible (compact object mass gap). It is also clearly seen in the right figure that the main accretion number contribution comes from low-mass stars, approximately up to ∼ 1.4 solar masses (according to the Chandrasekhar limit).

In the text