Globular Clusters Contribute to the Nuclear Star Cluster and Galaxy Center 𝛾 -Ray Excess, Moderated by Galaxy Assembly History

Two unresolved questions at galaxy centers, namely the formation of the nuclear star cluster (NSC) and the origin of the 𝛾 -ray excess in the Milky Way (MW) and Andromeda (M31), are both related to the formation and evolution of globular clusters (GCs). They migrate towards the galaxy center due to dynamical friction, and get tidally disrupted to release the stellar mass content including millisecond pulsars (MSPs), which contribute to the NSC and 𝛾 -ray excess. In this study, we propose a semi-analytical model of GC formation and evolution that utilizes the Illustris cosmological simulation to accurately capture the formation epochs of GCs and simulate their subsequent evolution. Our analysis confirms that our GC properties at 𝑧 = 0 are consistent with observations, and our model naturally explains the formation of a massive NSC in a galaxy similar to the MW and M31. We also find a remarkable similarity in our model prediction with the 𝛾 -ray excess signal in the MW. However, our predictions fall short by approximately an order of magnitude in M31, indicating distinct origins for the two 𝛾 -ray excesses. Meanwhile, we utilize the catalog of Illustris halos to investigate the influence of galaxy assembly history. We find that the earlier a galaxy is assembled, the heavier and spatially more concentrated its GC system behaves at 𝑧 = 0. This results in a larger NSC mass and brighter 𝛾 -ray emission from deposited MSPs.


INTRODUCTION
A compact bright star cluster is commonly observed at the centers of galaxies of all types, known as the nuclear star cluster (NSC) (e.g.Light et al. (1974); Kormendy (1985); Matthews & Gallagher (1997); Hughes et al. (2005)).They take up the innermost a few to tens of parsecs at most, and have masses 10 5 ∼ 10 8  ⊙ (Georgiev et al. 2016;Spengler et al. 2018), making them the densest known star clusters.Observationally, they are distinguished by prominently brighter luminosity on top of the disk or bulge component (e.g.Böker et al. (2002); Kim et al. (2004)).Besides, some NSCs are also observed to co-exist with a supermassive black hole (SMBH) at the galaxy center (e.g.Filippenko & Ho (2003); Nguyen et al. (2019)).
NSCs consist of a mixed stellar population in terms of age, metallicity, etc (e.g.Böker et al. (2001); Kacharov et al. (2018)).The complexity in stellar population has complicated the quest of NSC formation mechanisms.The young and metal rich stars suggests an in-situ formation scenario, where local star formation is triggered by the inflow of gas induced by various mechanisms, such as bar-driven infall (Shlosman et al. 1990) and the action of instabilities (Milosavljević 2004).The fact that this young stellar population is usually flattened, rotating, and concentrated at the center of NSCs (Georgiev ★ E-mail: hliastro@tsinghua.edu.cn& Böker 2014;Carson et al. 2015) also favors the in-situ formation scenario.On the other hand, the old and metal poor population can naturally arise from massive globular clusters (GCs) which migrate into the galaxy center due to dynamical friction.This was actually the first proposed NSC formation scenario (Tremaine et al. 1975) one year after the ground-breaking observation of the NSC in M31 (Light et al. 1974).Observationally, this scenario is also supported by evidences such as the deficit of massive GCs (Capuzzo-Dolcetta & Mastrobuono-Battisti 2009) and the nucleation fraction tracing the fraction of galaxies that have GCs (Sánchez-Janssen et al. 2019).
Despite observational motivation for both NSC formation mechanisms, a comprehensive modeling is lacking.For in-situ formation, direct simulation is challenging.There are limited studies focusing on different aspects of the process, such as gas inflow (Hopkins & Quataert 2010), momentum feedback and self-regulation (McLaughlin et al. 2006), and stellar two-body relaxation (Aharon & Perets 2015).Ex-situ formation was more extensively studied (e.g.Tremaine et al. (1975); Capuzzo-Dolcetta (1993); Lotz et al. (2001); Gnedin et al. (2014); Leveque et al. (2022)).In particular, Hartmann et al. (2011); Antonini et al. (2012); Tsatsi et al. (2017) used direct n-body simulation to study the final parsec-scale evolution of spiraled-in GCs and resulting NSC morphological and kinematic properties.The general picture has been well-established, but previous studies adopted crude treatments of the initial conditions of GCs due to a lack of knowledge on their formation.Besides, the evolution of the old GC systems is closely correlated with the evolution of the host galaxy, which hasn't been taken good care of in previous studies.Up to today, it is still uncertain the roles of in-situ and ex-situ formation mechanisms (e.g.Guillard et al. (2016); Fahrion et al. (2022)), with ongoing debates on topics such as a potential transition between the two mechanisms indicated by the NSC mass or galaxy stellar masses (Lyubenova & Tsatsi 2019).
While the spatial distribution of DM is already expected to peak at galaxy centers (e.g.Navarro et al. (1997)), MSPs are not commonly observed there due to difficulties in resolving individual -ray sources.However, MSPs are abundant in GCs, with much more observed number per unit mass compared to the galaxy field (Ransom 2008;Brandt & Kocsis 2015;Ye et al. 2019).This is due to the fact that MSPs are believed to originate from low-mass X-ray binaries (LMXBs), where the neutron star gets spun up through mass transfer from the companion star (Bhattacharya & van den Heuvel 1991).Thus, the high stellar density of GCs provide the desirable environment for both primordial binary formation and dynamical encounter.Galaxy center MSPs can originate from GCs that have migrated in and tidally dissolved.
On the other hand, the central bulge region of galaxies is also a dense stellar environment, although less dense by an order of magnitude than most GCs (Voss & Gilfanov 2007a).Thus, MSPs could in principle form in-situ as well.However, recent studies using scaling relations to probe in-situ MSP luminosity cast doubts on this mechanism as the sole origin of the -ray excess.The Galactic center excess (GCE) has been examined by Cholis et al. (2015) and Haggard et al. (2017), who pointed out that LMXBs are too rarely observed in the MW bulge that in-situ MSPs can account for less than a quarter of the excess luminosity.For M31, where ∼ 10 4 LMXBs are needed to explain the excess, less than 80 were detected within the inner 12 arcmin (∼ 2.7 kpc) (Voss & Gilfanov 2007a).Consequently, it is crucial to explore the contribution of ex-situ MSPs from GCs.
From the previous analysis, we see that the evolution of GCs inevitably contributes to NSC formation and the -ray excess, and it is important to find out the extent of its contributions.However, the formation of GCs still remains highly uncertain.Many previous studies (e.g.Gnedin et al. (2014); Fragione et al. (2019); Leveque et al. (2022)) assumed that all GCs formed at a single redshift, which only serves as a primitive approximation.To achieve a more reliable modeling of GC evolution and mass deposition, a better formulation of GC formation is needed.
In this paper, we use a new semi-analytic model of GC formation and evolution to study its contribution to the NSC formation and galaxy center -ray excess in galaxies similar to the MW and M31.The model adopts the GC formation scenario by Li & Gnedin (2014); Choksi et al. (2018); Choksi & Gnedin (2019), where GC formation was triggered by periods of rapid mass accretion onto the host galaxy across its assembly history, typically triggered by major galactic mergers.To obtain realistic galaxy merging histories, results from the Illustris cosmological simulation are used, and GCs are sampled at qualified simulation snapshots.After formation, the GC population is subject to orbital migration and mass loss depending on their mass and galactocentric distance.In addition, we also model an evolving background potential according to the galaxy assembly history.Through this new model, we hope to enhance our understanding of the connection between GCs,the NSC and galaxy center -ray excess .
We arrange this paper as follows.We introduce our modeling of GC formation and evolution in Section 2. In it we also show how to account for the galaxy center MSP luminosity at  = 0.The calculation of halo parameters from Illustris outputs is introduced in the Appendix A. In Section 3, we present our model predictions of GC properties, the NSC mass and -ray emission by MSPs at  = 0.As we can retrieve from Illustris a collection of halos of similar masses but with different assembly histories, we discuss the moderation effect of assembly history in Section 4. As our -ray luminosity prediction for the M31 falls short to observation, we also discuss alternative explanations.Caveats of our study are listed and discussed as well.Finally in Section 5, we summarize important results and suggest future work.

METHODS
In this section, we present our semi-analytical model for the formation and evolution of GCs in the framework of hierarchical structure formation.Then we describe the calculation of the -ray luminosity at  = 0 from deposited MSPs.

Formation of GCs in cosmological simulations
We introduce our modeling of GC formation in terms of formation times, initial masses and spatial distribution.
In modeling GC formation times, we improve on the simple prescription by previous studies (e.g.Gnedin et al. (2014); Fragione et al. (2019); Leveque et al. (2022)) that GCs formed at a single redshift.While this assumption is partly justified due to the old ages of most GCs, it is recently recognized that GC formation covers a wide range of cosmic time with diverse formation histories (Forbes et al. (2018) and references therein).Therefore, a more physically-motivated GC formation model is desired.Fortunately, over the past decades, our understanding of the origin of GCs in the framework of hierarchical structure formation has been revolutionized.Here, we adopt the GC formation model Choksi et al. (2018) (hereafter CGL), which was built upon earlier works by Muratov & Gnedin (2010) and Li & Gnedin (2014).The CGL model assumes that GC formation was triggered by periods of rapid mass accretion onto the host galaxy, typically during major mergers.This idea was motivated by multiple reasons such as more observed young massive clusters in interacting galaxies (e.g.Wilson et al. (2006); Portegies Zwart et al. (2010)), earlier formation times of GCs than the field stars and that galactic mergers were more frequent at high redshifts (Li & Gnedin 2014), and that galactic mergers are able to induce the high densities and pressures desired for cluster formation (e.g.Li et al. (2017);El-Badry et al. (2019)).
In the CGL model, GC formation is painted onto the halo merger trees of the Illustris simulation, which captures the evolution of halo properties from =47 to 0 (Vogelsberger et al. 2014;Nelson et al. 2015).GC formation is triggered when the specific halo mass accretion rate exceeds a threshold value, which is a tunable parameter.The total GC mass is mapped from the halo mass through observed stellar mass-halo mass relation (SMHM) (Behroozi et al. 2013), the stellar mass-gas mass relation (Magdis et al. 2012;Lilly et al. 2013;Genzel et al. 2015;Tacconi et al. 2018), and GC mass fraction from total gas mass (the second tunable parameter).Once the total GC mass at the epochs of formation is fixed, individual GC masses are sampled from a power −2 mass function observed from young massive clusters (Portegies Zwart et al. 2010;Schulz et al. 2015).It is encouraging to note that this simple two-parameter model successfully reproduced several key observed properties of GC populations, such as metallicity bimodality, GC mass-halo mass relation, etc.This serves as a much improved model of GC formation for the study of NSCs in galaxies with different assembly histories.
As GCs form at multiple epochs across the galaxy assembly history, their initial spatial distribution varies and influences their subsequent orbital evolution.Thus, we need to carefully take care of this issue at each formation epoch.For GCs formed inside the galaxy ('insitu'), we assume that, similar to normal star formation, GCs follow similar spatial distribution of cold gas during their formation.As gas does not exhibit a bulge structure, we adopt the continuous spherical Sérsic density distribution, based on which GCs are assigned to specific galacto-centric radii.On the other hand, for GCs formed inside satellite galaxies and were brought in via galactic mergers ('ex-situ'), we lack detailed knowledge of the stellar dynamics during galactic mergers.As a simplified treatment, we place these ex-situ GCs at half the virial radius of the descendant halo.This shouldn't affect our GC distribution at the galaxy center at  = 0, though, because these massive GCs carry large angular momenta relative to the descendant halo, which are highly unlikely to be sufficiently reduced by dynamical friction from the ambient stellar density.
The Sérsic spatial density distribution was proposed by Prugniel & Simien (1997) to match the well-established Sérsic surface brightness profile (Sérsic 1963).The spherical density distribution is given by: Here  0 is a normalization,   is the effective radius of the galactic disk,   is the concentration index.The term  is a function of   to ensure half of the projected light is contained within   , and can be well approximated analytically by  = 2  − 1/3 + 0.009876/  for 0.5 <   < 10 (Prugniel & Simien 1997).The form of  is adopted from Márquez et al. (2000) as  = 1.0 − 0.6097/  + 0.05563/ 2  for 0.6 <   < 10 and 10 −2 ≤ /  ≤ 10 3 to match the Sérsic surface brightness profile (Terzić & Graham 2005).
The effective radius   is related to the virial radius of the halo,  vir , and halo spin parameter, , by assuming a classical galactic disk formation model (e.g.Mo et al. (1998)): (2) Since we do not have the information of the total energy of the dark matter halo to estimate the traditional spin parameter defined in Peebles (1980), we instead use an alternative definition  B =  sp / √ 2 vir  vir that only requires the specific angular momentum and virial velocity  vir (e.g.Bullock et al. (2001)).Both  vir and  sp are provided at snapshots of the Illustris simulation, and  vir and  vir are calculated accordingly as explained in the Appendix.
Regarding the concentration index   , larger values were conventionally associated with higher concentration, but we found this to be misleading.In Fig. 1, we present the cumulative mass fraction distribution versus normalized radial distance for different   .We see that larger   does exhibit higher concentration in the inner region.However, beyond the crossing point at  ⪆   ( which we refer to as Sérsic crossing hereafter), they reach saturation much slower than curves with smaller   .Consequently, while larger   values indicate a more peaked distribution, they also indicate a greater spread.In our subsequent studies, we investigate the influence of different   values on our results within the range of 0.5 to 4, with increments of 0.5.As a fiducial model, we select   = 2 following the work of Gnedin et al. (2014).
To study the GC system in the MW and M31, we select Illustris halos with similar masses at  = 0.For the MW, we adopt the findings by Deason et al. (2021), who estimated a virial mass of 1.01±0.24×10 12 ⊙ .Among the Illustris halos, 1099 fall within this range.For M31, there are significant uncertainties, therefore we adopt a wider range of 0.7-2.5 ×10 12  ⊙ based on the review by Kafle et al. (2018).Accordingly, 2029 halos were selected, including those selected for the MW.

Dynamical evolution of GCs
After birth, GCs migrate towards the galactic center due to dynamical friction while depositing masses due to stellar evolution and tidal stripping along the way.In situations where they have migrated into the innermost region, the tidal force can be so strong that GCs get completely torn apart (Arca-Sedda & Gualandris 2018; Wang & Lin 2023).Thus, GCs that make all the way into the inner a few parsecs can build up NSC and contribute to the MSP populations there.We directly adopt the analytical prescriptions of Fragione et al. (2019) in modeling: 1) orbital migration, 2) tidal stripping, and 3) direct tidal disruption.These prescriptions include corrections to parameters originally proposed by Gnedin et al. (2014) and Fragione et al. (2018).Below, we describe two improvements we made in this work.
The first improvement is an updated prescription of the mass loss due to stellar evolution, which is obtained from the stellar population synthesis model FSPS (Conroy & Gunn 2010).This allows us to account for the evolving nature of stars within GCs more accurately.
The second improvement is modeling the time-varying gravitational potential by extracting the time-evolution of the dark matter halos from the Illustris simulation.For any cosmic time, we linearly interpolate the halo mass and spin between adjacent Illustris snapshots, and calculate galactic structural parameters accordingly.The overall gravitational potential comprises three components: the dark halo, the stellar component, and modeled GCs.The dark halo is described using the NFW distribution (Navarro et al. 1997), and the calculations for determining halo parameters are introduced in the Appendix A. The stellar component is modeled with a Sérsic distribution, as explained earlier.The mass of modeled GCs is also included in calculating the overall potential.For cosmological parameters, we adopt a flat ΛCDM model with ℎ = 0.704 and Ω m,0 = 0.2726, consistent with the Illustris simulation.
To improve the efficiency of our simulation, we implement subcycling in the evolution of GCs.We divide the entire time span into 100 sections, each characterized by a constant background potential.To evolve individual GCs, we first calculate their evolution timestep   being the smaller of the tidal and orbital evolution timescales multiplied by a fraction, ts m and ts r , respectively.Then, a GC with the smallest value of t i + dt i means that it evolves the fastest, and exerts the strongest influence on other GCs.Thus, we find and evolve such GC at each step, until all GCs cross the current time span.
Meanwhile, we found that it is usually the first step of calculated dt i that is overestimating.To efficiently address this, we introduce a maximum cutoff timestep dt max on top of the ts factors.By testing single GCs with different masses and galacto-centric positions, we optimize the values of dt max , ts m and ts r together.We find that dt max = 0.01 Gyr and ts m = ts r = 0.2 strike a balance between efficiency and accuracy.
For calculating the change in GC mass and galactocentric distance in each step, we employ the Runge-Kutta 4th order method.This method offers a significant speed improvement of approximately 20 times compared to the 1st order Euler method.

MSPs from GCs
MSPs are deposited by GCs due to tidal stripping and disruption, thus they can contribute to the unresolved -ray excess.In calculating the -ray luminosity at  = 0 from MSPs, we follow Fragione et al. (2018Fragione et al. ( , 2019) ) who set the total -ray luminosity of deposited MSPs equal to that of the debris of the GC, using luminosity-mass relation of GCs fitted to observations.To account for uncertainties in this fitting, they also used a constant luminosity-mass relation for comparison, and labeled it as the model 'C'.The model with fitted luminosity-mass relation was labeled 'EQ'.
Then, individual MSP luminosity was sampled according to the observed MSP luminosity function.To account for MSP spin-down due to the loss of rotational energy via magnetic dipole braking, Fragione et al. (2018,2019) use two models of the spin-down timescale: a Gaussian distribution model (GAU) and a log-normal distribution model (LON).In our study, we follow their prescribed methodology and calculate the -ray luminosity at redshift  = 0 using the four models, namely GAU-C, GAU-EQ, LON-C, and LON-EQ.

RESULTS
In this section, we present our results on the properties of GCs at  = 0 for MW and M31-like galaxies, and compare with observations.We then show the predicted mass of the NSC and -ray luminosity distribution from MSPs.We also explore how the above quantities vary with different galaxy assembly histories.

GC-halo mass scaling relation
One of the most striking observations of GCs is a linear correlation between the mass of the GC system and its host halo (Spitler & Forbes 2009;Georgiev et al. 2010;Hudson et al. 2014;Harris et al. 2017).The mass ratio is approximately 10 −5 across 5 orders of magnitude (Harris et al. 2015).Therefore, it is crucial for our model to reproduces this scaling relation.Fig. 2 shows our model prediction of the scaling relation for the MW and M31-like halos combined for different  s .We see that the runs with smaller  s tend to produce lower mass in the GC system at  = 0.This is in line with our findings in Fig. 1, that a smaller   corresponds to a smaller spatial spread of formed GCs, which results in more mass loss due to stronger tidal effect.Nevertheless, the changes of total GC mass in different  s is very small and all choices of  s produce the GC mass in reasonably good agreement with observations.Consequently, our model successfully reproduces the GC-halo mass scaling relation.

Spatial distribution of GC number density
In this section, we closely examine the spatial distribution of GCs in their number density at  = 0. Fig. 3 shows the average trend for MW and M31 combined, as we have checked that they are barely distinguishable.
When we compare the initial and final GC distributions, we notice that the number density outside ≈3 kpc barely decreases.This can be attributed to the fact that beyond this position, most GCs with a mass of approximately 10 5  ⊙ have tidal and migration timescales longer than a Hubble time.In addition, stellar evolution alone cannot exhaust a GC.
When we investigate the effect of different  s values, it turns out that this made minimal difference for both the initial and final GC distributions, except for the innermost region.In the initial distribution, the typical Sérsic crossing is only observed ≈100 pc, since GC formation across cosmic times overlaps in the outer regions as the halo grows, erasing outer crossings.Our inclusion of ex-situ GCs The surface density distribution of GC number with different  s averaged for all selected halos.We do not separate results for MW and M31like halos, because we have checked that they are barely distinguishable, and that the MW mass range is enclosed by that of the M31.In this plot, the redish curves denote all initial GCs formed at different redshift, and the bluish denote final GCs at  = 0. We show in red dotted curve (G14) the initial GC distribution in Gnedin et al. (2014) for comparison.The error-bars (B21, purple; RBC ver.5, black) show the observed GC catalogue for the MW (Baumgardt et al. 2021) and M31 (Galleti et al. 2014).In B21,  was manually removed as it was proved to be the stripped core of a disrupted dwarf galaxy (Noyola et al. 2008).also contributes to this effect.In the final distribution, although the difference between  s values remains minimal, the trend is opposite to that observed in the initial distribution.Larger  s values within the Sérsic crossing lead to stronger tidal disruption, resulting in fewer surviving GCs and a smaller number of GCs in the final distribution.
Comparing our results with observations, we found a minor overshoot within ≈ 5 kpc for the MW.However, individual halos exhibit a spread around the average distribution, and several candidate halos demonstrate conformity to the observation.We present one such candidate halo in Fig. 4. Therefore, our model can be considered an acceptable fit to the observed spatial distribution of MW and M31 GCs.Additionally, we observed no preference for  s , so we will keep  s = 2 as the fiducial choice.
To further analyze the GC population, we utilize the orbital information in the Baumgardt et al. (2021) catalogue to distinguish between the in-situ and ex-situ sub-populations of MW GCS.In order to do so, we follow Massari et al. (2019), who analyzed the dynamical differences of the two branches of GCs on the age-metallicity plot for 151 MW GCs.Based on their distinct features, they assigned 62 GCs as having formed in-situ, which comprise of bulge GCs and disk GCs.The former are defined as having apocenter distance less than 3.5 kpc, while the latter as having orbital altitude less than 5 kpc and circularity greater than 0.5.For our GC catalogue, we calculate the orbital and potential parameters using the ℎ package   by McMillan (2017); Dehnen & Binney (1998).Since our treatment of ex-situ GCs is rudimentary, we focus on in-situ GCs for now.
In Fig. 5, we plot the number density distributions of MW insitu GCs in a similar manner to Fig. 3.We can observe that most of the observations made for the overall GC distribution also apply to the in-situ GC distribution, except for 2 differences.First, in-situ GCs exhibit a more concentrated distribution, with a significant drop in their numbers towards the outskirts.This suggests that in-situ GCs preferentially populate the central regions of the MW-like ha- los.Second, without the contribution from the ex-situ population, the dispersion between different values of   becomes more pronounced in the outskirts, and we can observe the outer portion of the Sérsic crossing.Nevertheless, we can regard our model as satisfactorily reproducing the observed GC number density distribution, both overall and in-situ, and   =2 as a reasonable choice for the concentration parameter.
One advantage of utilizing the halo catalog extracted from Illustris is the ability to investigate the influence of different assembly histories on halos with similar masses.To parameterize the halo assembly histories, we employ the concept of half mass redshift, denoted as  hm , which represents the redshift at which the halo acquired half of its present mass.In Fig. 6, we show the distribution of  hm values.The histogram exhibits a log-normal shape centered around  ∼ 1.2, corresponding to a look-back time of approximately 8 billion years, with a tail extending towards higher redshifts.
We are interested in whether  hm has any discernible effect on the spatial distribution of GCs.Fig. 7 illustrates the GC number density distribution with different  hm for MW and M31-like galaxies.When we look at the influence of  hm for initial GCs, galaxies that 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 formed earlier (referred to as EFGs, earlier formed galaxies) tend to exhibit more concentrated distributions of GCs compared to halos that formed later (referred to as LFGs, later formed galaxies).This trend can be attributed to the higher merger rate of halos at larger redshifts, as indicated by previous studies (Fakhouri & Ma 2008;Fakhouri et al. 2010).Consequently, EFGs experience more GC formation events, starting from smaller halo sizes, while LFGs undergo fewer GC formation events, each resulting in significant growth of halo mass and size and leading to a more spread-out distribution of formed GCs.
In the case of final GCs, however, the impact of  hm is less pronounced.This can be understood, as more concentrated GCs also experience stronger tidal disruption and dynamical friction, which reduces their numbers.Additionally, as these GCs formed earlier, these disruptive processes act over a longer period.Consequently, the final GC number density distribution shows minimal traces of the galaxy assembly history.
When comparing our results with observations, as the influence of  hm on the final GC distribution is not significant, the results are similar to those of Fig. 3.

GC mass function
Besides the spatial distribution of GCs, their mass function at  = 0 is also an important property to compare with observations.Previous studies using the same GC formation model (Muratov & Gnedin 2010;Li & Gnedin 2014;Choksi et al. 2018;Choksi & Gnedin 2019;Chen & Gnedin 2022, 2023) have carried out comprehensive analyses, confirming a transformation from the initial power-law mass function to a log-normal shape at  = 0 that agrees with observations.Therefore, in this study we focus on the in-situ GC population.In Fig. 8 we show the median mass function of candidate MW-like halos for different   , and overplot observation results.We see that our model also shows a log-normal shape, although with its peak shifted towards smaller GC masses compared with observations.This can be due to the incompleteness in observed low mass GCs.Nevertheless, the shaded region corresponding to   = 2 can cover most of the observation trend.Indeed, candidate halos can fit the observation quite well, as illustrated by Fig. 9. Thus, our model provides a good fit to the observed in-situ GC mass function.
We also observe that   results in more noticeable discrepancies in the predicted mass function above ∼ 10 5  ⊙ , while at lower masses, the influence of different   is flipped and less prominent.Larger   leads to slightly fewer light GCs but more heavier ones, and vice versa.This is due to the fact that those light GCs are prone to weak dynamical friction, thus they barely migrate inward.Their disruption is mostly determined by their position at birth.As we have observed in Fig. 1, larger   have peaked density at the galaxy center, which leads to prominent disruption of lighter GCs that initially formed in the region.On the other hand, as the general distribution is more spreaded, heavier GCs have a longer journey to migrate to the galaxy center and are subject to weaker tidal effect.Thus, the distribution leads to slightly stronger disruption of light GCs formed in the innermost region, but weaker disruption of heavier ones.

NSC mass
In this section, We examine the NSC mass contributed by GCs as they spiral in.On average, ∼ 80 and 100 GCs migrated into the NSC of the MW and M31 respectively, with smaller   corresponding to a few more GCs.The standard deviation is ∼ 30 and 40 respectively.In Fig. 10, we compare the deposited GC masses with the observed NSC mass at  = 0. Unlike the number density distribution, the cumulative mass distribution of initial GCs exhibits a prominent Sérsic crossing in the inner region, because the range of the y axis is much smaller in this case.Regarding the deposited mass, we observe that GCs with smaller values of  s contribute more to the total deposited mass, which aligns with the trend we observed in the GC-halo mass scaling relation discussed in Section 3.1.Regardless of the influence of  s , the deposited mass is predominantly confined to the central regions.This contrasts with the results presented by Gnedin et al. (2014); Fragione et al. (2018), where a plateau is established starting from 4 pc outward and continues to rise prominently up to 10 kpc.The disparity arises due to our GC formation occurring throughout the entire assembly history of the halo, resulting in a more concentrated distribution of GCs, as shown in Fig. 3. Consequently, GCs are more susceptible to significant tidal effects and tend to deposit mass towards the galaxy center.
While the average deposited mass exceeds the NSC mass, it is important to note that different galaxy assembly history yield substantial variations, as depicted in Fig. 11.As discussed in the preceding section, EFGs give rise to more GCs, which experience stronger and more prolonged tidal disruption.Consequently, halos with larger  hm have GCs depositing more mass towards the center.Conversely, LFGs exhibit less deposited mass, implying that our MW NSC plausibly originates from a halo with  hm ∼ 1.

Gamma ray luminosity
As our GC fittings turn out consistent with observations, we move forward with  s = 2 to check the spatial distribution of the luminosity of deposited MSPs.As different authors analyzed the -ray excess data in different methods, we present our results in both differential and cumulative distribution.Fig. 12 shows our model prediction of the differential flux distribution.Notably, different models yield similar results.And although they generally underestimate the excess flux, the overall shape is consistent.Consequently, we select the GAU-EQ model as the best choice and plot it for different  hm intervals in Fig. 13.We observe that EFGs exhibit relatively higher flux emission,  as their GCs deposit more MSPs.Halos with  hm ≳ 0.7 provide a good fit to the observations.Note that for clarity purpose, the color code here does not precisely match the colorbar in Fig. 11.With that in mind, we notice that the results are compatible, indicating that a halo with  hm ≳ 0.7 can successfully reproduce both the NSC mass and the spatial distribution of the -ray differential flux emmision.
For the cumulative flux distribution, Fig. 14 illustrates the GAU-EQ model for different  hm .We can observe that our model also shows a consistent shape with the observations, which appears better than the results by Fragione et al. (2018) as we exhibit more flux in the innermost region.This arise from our GC formation occurring throughout the entire assembly history of the halo.As a result, GCs started their evolution closer to the galaxy center.Once again, halos with  hm ≳ 0.7 provide a good fit.And consequently, it is plausible that the GCE arose solely from deposited MSPs.Having examined the GCE, we proceed to the M31 excess.Due to the lack of detailed analysis regarding the excess spatial distribution, we solely compare the cumulative luminosity at 6 kpc.Fig. 15 illustrates the cumulative luminosity distribution for the four MSP luminosity models.The models exhibit notable discrepancies, particularly in the innermost region, which diminish as we move towards the outskirts and nearly vanish at 6 kpc.This discrepancy stems from the fact that the EQ models employ a fitted relation, where log(  /  ) is negatively proportional to log(  ).Consequently, heavier GCs are dimmer compared to the C models that utilize a constant luminosity-mass ratio.Since heavier GCs are more susceptible to dynamical friction and tidal disruption, they primarily contribute to the -ray emission in the innermost regions.As we move farther away from the galaxy center, lighter GCs become more dominant, which reduces the discrepancy.Nevertheless, all models  converge and fall short of matching the excess signal at 6 kpc.Additionally, the highest luminosity among individual halos only reaches approximately 7 × 10 37 erg/s, which is less than one-third of the excess luminosity.Therefore, it is unlikely that MSPs alone account for the M31 excess.This suggests that the two galaxies have distinct origins of the excess emissions despite their similar masses.

Galaxy assembly history
We have observed in preceding sections that galaxy assembly history has clear influences on the GC and NSC properties.We shall take a systematic look in this section.
In Fig. 16 we illustrate the correlations between  hm and four important masses: the halo mass, NSC mass, total initial GC mass and GC mass at  = 0.The halo mass does not correlate with  hm since halos of any mass can assemble early or late.However, the density of data points varies across  hm , consistent with the lognormal distribution presented in Fig. 6.And we have demonstrated previously that EFGs give rise to more GCs, resulting in a positive linear trend between the initial GC masses and  hm .We have also showed that EFGs have GCs distributed closer to the galaxy center, making them susceptible to stronger and prolonged tidal disruption.Consequently, they contribute a greater amount of mass to the NSC, leading to a positive linear trend between the NSC mass and  hm .Thus, the NSC mass serves as a good indicator in breaking the degeneracy to infer the galaxy assembly history.Larger NSC masses indicate relatively earlier accumulation of the halo mass, and vice versa.However, it is intriguing that  hm appears to provide little information about the final GC mass.This seems to suggest a lack of correlation between initial and final GC masses either, contrary to what one might expect.
To verify this observation, we show in Fig. 17 the initial and final GC masses for all MW-like halos, with colors denoting  hm .A positive trend reasonably exist between initial and final GC masses, although there is relatively large dispersion, particularly at smaller final GC masses.If we observe  hm across different initial GC masses, larger GC masses do correspond to larger  hm .However, if we examine final GC masses, each value is associated with a large spread of  hm values.This is due to the fact that the earlier-formed larger GC masses also suffer from stronger disruption, potentially resulting in smaller final masses at  = 0. Consequently, small final GC masses can arise from either small initial masses or earlier-formed large initial masses.This breaks the correlation between final GC masses and galaxy assembly history.
We are also interested in correlations associated with the NSC mass, a highly significant outcome of our GC model.Fig. 18 shows the NSC mass plotted against the halo mass, initial GC mass, and final GC mass.As previously demonstrated in Fig. 16, the NSC mass is positively correlated with  hm , which has no relation to the halo mass.It is not surprising to observe that the NSC mass does not correlate with the halo mass.It is the earlier assembly of the halo that give rise to a heavier NSC, rather than the mass of the halo itself.However, it is worth noting that this observation may change when examining a broader range of halo masses, which will be investigated in our future studies.
In the middle panel of Fig. 18, a strong positive trend is observed between the mass of initial GCs and the NSC mass.Not only do GCs serve as the fuel for the build-up of the NSC, but larger initial GC masses also correlate with larger  hm values, indicating earlier formation and more concentrated distribution.Collectively, these factors contribute to a larger NSC mass.
However, final GC masses exhibit no correlation with the NSC mass.This can be comprehended, as we have demonstrated in Fig. 16 that final GC masses do not correlate with  hm .Each GC mass at  = 0 might arise from either early-formed large   which builds up a large NSC, or later-formed lighter   which contributes little to the NSC mass.Thus, the final GC mass does not serve as an indicator for the NSC mass.
In summary, we presented in this section various correlations associated with the galaxy assembly history.We found that  hm , the total initial GC mass and the NSC mass are correlated, due to the fact that EFGs give rise to an old, heavy and concentrated GC system which contributes a larger amount of mass to the NSC.Thus, either value among  hm , the initial GC mass and the NSC mass serves as an indicator of the other two.And one important observation is that we 0.8 0.9 1.0 1.1 1.2 m halo ( 10 could utilize the NSC mass to infer knowledge on galaxy assembly history.Intriguingly, the final GC mass is not correlated with  hm or the NSC mass.This stems from the degeneracy in the relation between initial and final GC masses.

Possible explanations of the M31 gamma-ray excess
While our results fall short of the M31 excess by an order of magnitude on average, candidate halo can reach approximately one third of the signal at its highest.Nevertheless, even with in-situ MSPs combined, the MSP channel alone is unable to fully explain the M31 excess.Thus, it is evident that a contribution from DM is necessary.
Upon the first report on the detection of the M31 excess, Ackermann et al. (2017a) have brought up the possible explanation of DM.A primitive estimate inferred from a DM-origin GCE results in a flux deficit by 5 times, though the level of uncertainty was high.Subsequent investigations claimed to match the excess luminosity, but commonly identified tensions with observational constraints, such as the under-detection of DM emission in MW dwarf galaxies (Di Mauro et al. 2019) and a lack of DM radio emission for M31 (Mc-Daniel et al. 2018). Additionally, McDaniel et al. (2018) found that the two preferred DM annihilation channels for M31, namely  b and an even mixture of  b/ +  − , favor smaller DM masses compared to those suggested by the GCE.Consequently, it was suggested that DM alone does not explain the M31 excess.Combining these findings with our results, it becomes clear that a combination of MSPs and DM offers a promising and potentially inevitable way for explaining the M31 excess.However, the question remains as to why the MW and M31 have different origins for producing such excess emission.

Caveats and future works
In this section discuss the caveats in our model and improvements to be made in future works.
First, although we adopt the analytical expression of tidal disruption from Fragione et al. (2019), it was based on a static spherical galactic background with circular cluster orbits (Gieles & Baumgardt 2008).Our inclusion of the assembly history of galaxies, however, indicates a more complicated galactic background, especially at large redshifts when galactic mergers were more frequent.This complication is twofold.On one hand, the galactic background keeps varying, although in out treatment of linear interpolation, the variation is steady.Thus, the static approximation is not unreasonable.On the other hand, the process of galactic mergers, especially major mergers, perturbs the galactic environment and GC orbits.However, due to the vibrant nature and lack of knowledge on the process, we leave it to future research.The eccentricity of GC orbits also differs from the circularity approximation.However, to accurately capture realistic GC orbital evolution along galaxy assembly histories can be computationally expensive (e.g.Li et al. (2017Li et al. ( , 2018)); Chen & Gnedin (2022).We could only treat our method as a time averaged approximation to real eccentric GC orbits.However, as very eccentric orbits usually happens for ex-situ GCs, the approximation does not significantly affect our results on the properties of the NSC.With improving powers of N-body simulation on galactic mergers and dynamical friction, more knowledge will enable us to incorporate these effects into a more comprehensive model.
In investigating the formation of NSCs, we didn't take into account the in-situ channel of young stars forming in nuclear regions, despite various observational evidence as mentioned in 1.Our model partially takes care of this channel, as the GC spatial distribution at formation can sometimes sample GCs at the center of the host galaxies.Nevertheless, a more systematic investigation is warranted to obtain a thorough picture of NSC formation.
Besides fitting the NSC mass, future works will be carried out on investigating other NSC properties, such as the age/metallicity distribution and the internal mass profile of the NSC.Besides MW and M31-like galaxies, we will extend our investigation to a broader galaxy mass range and to different galaxy types.
There are also caveats in modeling the MSPs, as their distribution and evolution in GCs remains highly unknown to us.As stellar density increases towards the GC center, LMXBs and MSPs are supposed to peak near the center.There are evidences both from observations and simulation (see Ye et al. (2019) and references therein).Thus, the number of MSPs stripped and deposited to the ambient environment depends on the current size of the GC.Our treatment essentially assumes a uniform distribution of MSPs in GCs, which only serves as a lower limit to the MSP contribution to the galaxy center -ray excess.On the other hand, following Ye et al. (2019), Ye & Fragione (2022) studied the -ray excess problem by depositing MSPs only when the GC is fully disrupted.This methodology goes to another extreme by assuming all MSPs residing at the very center.Further knowledge on the distribution of MSPs inside GCs will enable us to better evaluate their -ray contribution to the galaxy center.Furthermore, we follow Fragione et al. (2018Fragione et al. ( , 2019) ) to assume the effect of new MSP formation canceling MSP spindown in GCs, which is somewhat arbitrary.More knowledge is required to properly evaluate these competing factors.

CONCLUSIONS
In this study, we have presented a comprehensive model of GC formation and evolution, based on the premise that they primarily form following periods of rapid halo mass accretion.Leveraging the results from the Illustris cosmological simulation, we sample GCs across the galaxy assembly history and simulate their subsequent evolution, accounting for the mass loss and radial migration within an evolving galactic background.Our model successfully reproduces key observations of the MW and M31 GC system at  = 0, including the mass scaling of the total GC system with the host halo, and the spatial distribution of GC number density.For the MW, we also reproduced the spatial distribution for the in-situ subpopulation.
With this model at hand, we investigate the spatial distribution of deposited masses of migrated GCs to study its link to the formation of the NSC and galaxy center -ray excess.We find that both NSC masses of the MW and M31 can be reproduced.Detailed spatial distribution of the GCE can also come entirely from deposited MSPs.However, the M31 excess strength is three times as large as our most luminous candidate galaxy.Even factoring in in-situ MSPs born at the galaxy center, the MSP channel still cannot fully account for the excess emission.It becomes evident that DM must play a role in explaining the M31 excess, highlighting a fundamental astrophysical difference between the two galaxies.This constitutes another big difference between them, apart from their galaxy center SMBHs differing in mass by about 50 times.Further investigations are demanded in figuring out the causes to and possible links between these differences.
Another intriguing aspect we discovered is the influence of galaxy assembly history on galaxy properties, which we investigated using halo half mass redshift  hm .Interestingly, we found that it does not correlate with halo mass, but conveys valuable information about the GC system and NSC mass.Specifically, EFGs with large  hm give rise to an old, heavy and concentrated GC population as they formed, and vice versa.This results in more deposited mass from GCs and a heavier NSC, which in turn serves as an informative indicator of galaxy assembly history.
In conclusion, our comprehensive model of GC formation and evolution provides a robust framework for understanding the properties of the GC system, the NSC and galaxy center -ray emissions.Additionally, our study unveils the significance of galaxy assembly history in shaping galaxy properties.However, the need to invoke DM to explain the M31 excess emphasizes the distinct astrophysical origins of these high energy emissions from the two similar galaxies.Further investigations are warranted to unravel the precise mechanisms that drive these differences and establish a comprehensive understanding of galaxy formation and evolution.

Figure 1 .
Figure 1.Distribution of the cumulative mass fraction (   ) for the Sérsic volume density profile (Eqn.1), plotted against the radial distance  normalized by the effective radius  e .Different concentration index  s are plotted in different colors according to the colorbar at the right.

Figure 2 .
Figure 2. The ratio between the total GC system mass ( GC ) and the host halo mass ( halo ) at  = 0 for different Sérsic index  s for all selected MW and M31-like halos.The black diamonds show the median ratio, while the error bars show the 25-75 interquartile range.For reference, we also overplot horizontal dashed lines for a compilation of recent observations from Spitler & Forbes (2009) (S09, blue), Georgiev et al. (2010) (G10, red), Hudson et al. (2014) (H14, cyan), and Harris et al. (2017) (H17, green).
Figure3.The surface density distribution of GC number with different  s averaged for all selected halos.We do not separate results for MW and M31like halos, because we have checked that they are barely distinguishable, and that the MW mass range is enclosed by that of the M31.In this plot, the redish curves denote all initial GCs formed at different redshift, and the bluish denote final GCs at  = 0. We show in red dotted curve (G14) the initial GC distribution inGnedin et al. (2014) for comparison.The error-bars (B21, purple; RBC ver.5, black) show the observed GC catalogue for the MW(Baumgardt et al. 2021) and M31(Galleti et al. 2014).In B21,  was manually removed as it was proved to be the stripped core of a disrupted dwarf galaxy(Noyola et al. 2008).

Figure 4 .Figure 5 .
Figure 4. GC number density distribution for one candidate MW halo with Illustris subhalo ID 00972.The red line shows its initial GCs, and blue line for GCs at  = 0.The error bars are again the Baumgardt et al. (2021) observed GC catalogue.

Figure 6 .
Figure 6.Histogram showing the half-mass redshift  hm for all candidate halos.Look-back time is shown on the upper horizontal axis as well.

Figure 7 .Figure 8 .
Figure 7. Number density distribution for all candidate halos with different  hm .For clarity purposes, we only show the average distribution for three  hm ranges.Dashed lines correspond to initial GCs and solid line to final ones.

Figure 9 .
Figure 9.The same as Fig. 8, but for one candidate halo with Illustris subhalo ID 17238.

Figure 10 .
Figure 10.Spatial distribution of the cumulative masses of initial GCs (summed over all formation epochs) and their deposition at  = 0, averaged for MW and M31-like halos with different  s .The color code is identical to Fig. 3.The MW SMBH and NSC observations are adopted from Do et al. (2019); Neumayer et al. (2020), and M31 from Bender et al. (2005); Leveque et al. (2022) respectively.

Figure 11 .
Figure 11.Cumulative spatial distribution of masses of the initial GCs and their deposition with different  hm for all candidate halos.Three dashed curves show the average for halos falling under different  hm spans.Observed positions and masses of the MW and M31 NSC and SMBH are plotted with uncertainties.

Figure 13 .
Figure 13.Similar to Fig. 12, but for the GAU-EQ model with three  hm ranges.

Figure 15 .
Figure 15.Cumulative luminosity distribution for the four MSP luminosity models for all 2029 M31-like halos.The average distributions for each model are also plotted.The Fermi data was adopted from Fragione et al. (2019) in red errorbar.

Figure 16 .
Figure16.Correlations of  hm with (from the top panel to the bottom) the mass of the halo, initial GCs, NSC and final GCs for all MW-like halos.Pearson correlation coefficients are shown at the upper left corner of each panel.Solid lines and shaded regions show the best-fit trend and 1- dispersion.The standard deviations of the four masses are 1.38×10 11 , 1.30×10 8 , 4.82×10 7 and 3.33×10 7  ⊙ respectively, which corresponds to 0.14, 0.54, 0.76 and 0.66 of the respective mean.

Figure 17 .
Figure17.GC initial and final masses for all MW-like halos, with colors denoting  hm .The best-fit trend with 1- band is plotted in black solid line and shaded region.

Figure 18 .
Figure 18.Correlations of the NSC mass with the halo mass (upper panel), total mass of initial GCs (middle panel) and that of final GCs (lower panel) for all MW-like halos.Pearson correlation coefficients are shown at the lower right corner of each panel, and best-fit trend and 1- range was shown as well.Only the initial GC mass exhibits a positive trend with the NSC mass.