Close stellar encounters at the Galactic Centre I: The effect on the observed stellar populations

We model the effects of stellar collisions and close encounters on the stellar populations observed in the Milky Way nuclear stellar cluster. Our analysis is based on N -body simulations in which the nuclear stellar cluster forms by accretion of massive stellar clusters in the presence of a supermassive black hole. We attach stellar populations to our N -body particles and follow the evolution of their stars and the rate of collisions and close encounters between different types of stars. We find that the most common encounters are collisions between pairs of main-sequence stars, which lead to merger; destructive collisions of stars with compact objects are relatively rare. We find that the effects of collisions on the stellar populations are small for three reasons. First, the core-like density profile in our N -body models limits the stellar density. Secondly, the velocity dispersion in the nuclear stellar cluster is similar to the surface escape velocities of the stars, which minimises the collision rate. Finally, whilst collisions between main-sequence stars can destroy bright giants by accelerating the evolution of their progenitors, they can also create them by accelerating the evolution of lower-mass stars. To a good approximation these two effects cancel out. We also investigate whether the G2 cloud that made a close pass by Sgr A* in 2014 could be a fuzzball: a compact stellar core which has accreted a tenuous envelope in a close encounter with a red giant. We conclude that fuzzballs with cores below about 2M have thermal times-scales that are too short to be compatible with observations of G2. A fuzzball with a black-hole core could reproduce the surface properties of G2 but the production rate of such objects in our model is low.


INTRODUCTION
Nuclear star clusters (NSCs) are among the densest stellar systems in the Universe (see Neumayer et al. 2020, for a recent and comprehensive review of the properties and evolution of NSCs). With half-light radii of a few parsecs and masses of 10 6 -10 7 M , they are located at the dynamical centres of the majority of observed galaxies (Böker et al. 2004;Côté et al. 2006;Böker 2010;Neumayer et al. 2011;Turner et al. 2012;Georgiev & Böker 2014;den Brok et al. 2014;Neumayer et al. 2020). Many of these NSCs coexist with a central supermassive black hole (SMBH, see e.g. Neumayer & Walcher 2012;Nguyen et al. 2019). This is the case for our Galaxy, whose NSC hosts at its centre SgrA*, an SMBH of 4.3 × 10 6 M (Ghez et al. 1998;Eisenhauer et al. 2005;Gillessen et al. 2009;Boehle et al. 2016;Gillessen et al. 2017). The Galactic NSC has a total mass of ∼ 2.5 × 10 7 M and a projected half-light radius of 4.2 pc (Schödel et al. 2014b,a). Thanks to multi-epoch high-quality adaptive optics observations, Schödel et al. (2020) found that 80% of the total stellar mass of the Galactic NSC formed 10 Gyr ago. After a quies-E-mail: mastrobuono@astro.lu.se cent phase that lasted around 5 Gyr, there has been another star formation episode that generated 15% of the observed stellar mass. Finally, a few percent of the mass of the NSC formed 100 Myr ago.
Although the NSC properties and evolutionary history seem to be correlated to those of their host galaxies, there is not yet a clear consensus on their formation mechanism. One hypothesis is that NSCs formed in situ, from gas that fragmented after flowing into the centres of the galaxies (Loose et al. 1982;Milosavljević 2004;Schinnerer et al. 2006Schinnerer et al. , 2008. This scenario has been tested using hydrodynamical and Nbody simulations (see e.g. Levin & Beloborodov 2003;Nayakshin & Cuadra 2005;Paumard et al. 2006;Hobbs & Nayakshin 2009;Mapelli et al. 2012;Mastrobuono-Battisti et al. 2019) that mostly focused on reproducing the disc of young stars observed within the central 0.5 pc of the Galactic NSC (see e.g. Genzel et al. 2010).
Another alternative mechanism, proposed by Tremaine et al. (1975), involves the merger of massive, dense stellar clusters, similar to globular clusters, that spiralled into the centre of the galaxy because of the action of dynamical friction (see also Capuzzo-Dolcetta 1993). This mechanism has been studied in detail using N -body simulations (see e.g. An-tonini et al. 2012;Gnedin et al. 2014;Arca-Sedda et al. 2015;Antonini et al. 2015). Some of the simulations were tailored to the Milky Way and could produce an NSC structurally and dynamically compatible with that observed at the Galactic centre (Antonini et al. 2012;Mastrobuono-Battisti et al. 2014;Tsatsi et al. 2017;Abbate et al. 2018). The two proposed scenarios could both partake in the formation of NSCs. It has also been suggested that clusters could be the potential source of gas from which new stars form in galactic nuclei (Guillard et al. 2016).
Using recent observations of the Galactic centre, Schödel et al. (2020) seem to constrain the contribution of classical globular clusters to a small fraction of the total mass. Recently, Feldmeier et al. (2014), Feldmeier-Krause et al. (2017) and Do et al. (2020) found a low metallicity substructure with kinematical properties that seem to suggest that it could have been generated by a recent stellar cluster merger at the Galactic centre (Tsatsi et al. 2017;Arca Sedda et al. 2020). However, a clear consensus on the origin of NSCs is still missing.
Following their dynamical evolution, stars around an SMBH should relax and reach a cusp-like distribution in density space (Bahcall & Wolf 1976). However, observations have shown that red giant (RG) stars show a core-like distribution in the inner 0.5 pc (12") of the Galaxy, with a flat or even radially decreasing projected profile (see e.g. Sellgren et al. 1990;Genzel et al. 1996;Buchholz et al. 2009;Do et al. 2009;Yusef-Zadeh et al. 2012). Several dynamical processes have been suggested to explain this so called "missing RG" problem. For example, Genzel et al. (1996) suggested that the core could be the result of collisions between RGs and main sequence (MS) stars. These kind of collisions should be frequent in the extremely dense central regions of the NSC. This scenario has been explored theoretically (Davies et al. 1998;Alexander 1999;Bailey & Davies 1999;Dale et al. 2009); these studies found that collisions would produce a NSC core too small with respect to the observations. Other mechanisms invoked to explain the discrepancy between models and observations include the action of an intermediate mass black hole or the infall of a star cluster (see e.g. Gualandris & Merritt 2012;Antonini et al. 2012). These scenarios imply that the Galaxy and the NSC went through processes not supported by other observables (e.g. a major merger or frequent and recent stellar cluster infalls). More recent works Schödel et al. 2018;Habibi et al. 2019;Schödel et al. 2020), using improved observational techniques, revised our knowledge of the stellar distribution in the central region of the NSC, finding a steeper density profile than seen in previous studies. In particular, Gallego-Cano et al. (2018) showed that fainter late-type stars have a spatial distribution well described by a power-law with index −1.43. However, brighter RG stars show a shallower cusp, with power-law index −1.2. This suggests that ∼ 100 RGs could be missing within 0.3 pc from Sgr A*. Similarly, Habibi et al. (2019) found a core-like structure only for the brightest RG stars. However, the number of missing RGs in their case seems to be smaller than 100. Based on results by Amaro-Seoane & Chen (2014), Amaro-Seoane et al. (2020) illustrated how the interaction between the stars in the cusp and the clumpy gaseous discs that generated new stellar populations at the Galactic centre could have stripped the envelope of the RGs making them invisible. Zajaček et al. (2020) recently suggested that the lack of RGs could be the result of the interaction between these stars and the nuclear jet that left the observed γ-ray Fermi bubbles as an imprint of its past action.
Another open question is the nature of G2, a cloud-like object orbiting the Galactic Centre. This cloud has been observed several times in various wavelengths (e.g. Gillessen et al. 2012;Phifer et al. 2013;Hora et al. 2014;Witzel et al. 2014;Bower et al. 2015;Pfuhl et al. 2015;Tsuboi et al. 2015;Valencia-S. et al. 2015;Plewa et al. 2017). There are two main classes of proposed origins for this object. In the former G2 is considered as purely a gas cloud. In the latter, G2 is the result of the ejection of gas from a compact source, like a planet or a star: see Calderón et al. (2018) and references therein for a review of different possibilities. In this paper we consider a formation channel where a close encounter between a red giant and a more compact star leads the compact star to accrete a tenuous envelope.
The remainder of this article is organised as follows. Sections 2 and 3 describe the contents of our models: N -body simulations of the formation of the NSC and the treatment of stellar collisions and close encounters. Section 4 describes the effects of encounters on the stellar populations in the NSC, and Section 5 the behaviour of compact stars that have accreted a tenuous envelope. Following a discussion in Section 5.2.1 we conclude in Section 6.

N -BODY SIMULATIONS
The N -body run used in this paper is described in Antonini et al. (2012), Mastrobuono-Battisti et al. (2014),  and Tsatsi et al. (2017). In these works we modelled the formation of a Milky Way-like NSC, through the inspiral and merger of twelve dense and massive stellar clusters, similar to globular clusters (GCs), in the centre of a nuclear bulge of mass 10 8 M . The Galactic centre initially hosts a supermassive black hole (SMBH) of 4×10 6 M , similar to Sgr A * (Ghez et al. 1998;Eisenhauer et al. 2005;Gillessen et al. 2009;Boehle et al. 2016;Gillessen et al. 2017). All clusters are described by a tidally truncated King (1966) model with total mass of 1.1 × 10 6 M , core radius rc = 0.5 pc and concentration parameter W0 = 5.8. The adopted simulation corresponds to 'Simulation 1' in Tsatsi et al. (2017). While the nuclear bulge is modelled using N b = 227523 equal mass (400 M ) N -body particles, every GC is modelled using NGC = 5715 super-particles of 200 M . We assume that the star clusters have already migrated to the inner Galaxy and that each cluster is initially located 20 pc from the central SMBH. All GCs are on randomly inclined circular orbits.
Our GCs decay consecutively at regular intervals of time (0.85 Gyr), over ∼ 10.2 Gyr. The system is then left to evolve without any additional GC infall until it reaches an age of 12.4 Gyr. At the end of the simulation, we obtain a central NSC that kinematically (including its rotation and velocity dispersion) and morphologically (e.g. mass, shape and density profile) resembles the one observed at the centre of the Milky Way.

Scaling
To obtain a more accurate match between the final mass and central density of the simulated and observed NSCs, we scaled our system using the observed values of the half-mass radius and total mass for the Galactic nucleus (Schödel et al. 2014b,a).
We defined the scaling factor for the distances as the ratio between the observed projected half-light radius (R h,obs = 4.2 pc) and the simulated projected half-mass radius R h,sim . Following Schödel et al. (2014b) we use the projected radial co-ordinate which accounts for the oblateness q. The observed NSC has q = 0.64 whilst our model has q = 0.72 (Tsatsi et al. 2017). We obtain The mass scale is the ratio between the observed NSC mass, M obs = 2.5 × 10 7 M , and the simulated NSC mass Msim: Finally, the velocities have been scaled so as to keep the virial radius Q of the simulated NSC constant after the other two scalings have been applied. Since where mi and vi are the superparticle masses and velocities, and ri their positions, we obtain the velocity scaling v scale = m scale l scale = 2.02.
The scaled number density is then defined as and the spatial and projected mass densities are and It is also necessary to consider how to treat the evolution of the particle masses. In the N -body code, the superparticles have a constant mass, but the mass in the coagulation calculations changes as the stellar population evolves. Essentially it is necessary to choose a time to set the N -body superparticle masses and the stellar mass in the coagulation code equal. Figure 1 shows the range of options available. The black line, plotted against the left-hand ordinate axis, shows the turnoff age in our models as a function of turnoff mass, whereas the blue curves (plotted against the right-hand ordinate axis) show the stellar mass at that age as a fraction of the initial population mass. Grey dotted lines show that the fractions of mass remaining at 10 8 and 10 9 years are 0.815 and 0.725 respectively. We choose to set the masses equal at 10 8 yr, since although that time precedes most of the dynamical evolution in the N -body code, a large fraction of the mass loss has already happened. Hence we define a second mass scale as the ratio between the initial stellar mass in each superparticle and the simulated mass, m scale,0 = 2.23.

Dynamical interactions
In order to make a rough estimate of the degree of dynamical interaction between stars we introduce the interaction parameter for N -body superparticle i, Γ ,i, defined as where ni(t) is the number density of stars in the particle i and σi(t) is the velocity dispersion local to the particle. The integral is carried out over all times from the introduction of the particle to the simulation until the end of the simulations at 12.4 Gyr. Γ can be understood as the collision probability for an individual star within the N -body particle, if it had at all times the mass and radius of the Sun. Whilst the actual collision probability for any given star must be calculated by a full integral over the time-varying stellar population, as described below, this gives a simple way of assessing how much interaction the stars that make up a given N -body particle are likely to experience. The final positions of the N -body superparticles, projected onto the y − z plane, are shown in Figure 2. The colour scheme tracks the value of Γ , as can be seen in the bottom panel. Whilst the particles are plotted in an order that systematically highlights the most interactive particles, it can be seen that these all lie towards the centre of the cluster, with the outer parts comprising entirely non-interactive particles. Hence we expect that the effects of close encounters and collisions will be most pronounced at the centre of the cluster. the interaction parameter Γ is plotted as a function of projected galactrocentric radius R. An offset of 0.1 pc is added to make the plot more compact. The colour scale is the same as in the top panel.

TREATMENT OF STELLAR COLLISIONS AND CLOSE ENCOUNTERS
In order to make a quantitative calculation of the effects of stellar evolution, collisions and close encounters on the stellar population in the Milky Way NSC, we attach a population of stars to each superparticle, calculate the rates of collisions between the various stellar species, and model the effects of the collisions that occur. The populations are coeval and each have initial mass where the superparticle mass in the N -body simulations, M spcl = 200 M (see Section 2) and m scale,0 is described in Section 2.1. The additional scale f overwight allows us to over-sample the IMF in order to avoid introducing problems from stochastic sampling. Empirically we have found that f overwight = 10 is sufficient, which corresponds to one or two black holes per superparticle. The stars are all taken to be single, and their birth masses of the stars are distributed according to the IMF of Kroupa et al. (1993).
Having built a stellar population, we evolve it using the dynamical history of that superparticle from the N -body simulations. We take the number density and velocity dispersion of the 100 superparticles closest to the particle being modelled and use that to compute the stellar interactions using the stellar population attached to the given superparticle. This effectively requires us to assume that the stellar population is only weakly affected by stellar encounters, an assumption that turns out to be quite good.
The number density ni of species i is then given by where n spcls is the number density of superparticles local to the superparticle being modelled and Ni is the number of stars of species i in the simulation. The velocity dispersions are taken to be Maxwellian, and the velocity dispersion σ is taken directly from an average over the 100 nearest superparticles in the N -body simulations. The instantaneous collision rate of an individual star i with mass Mi and radius Ri with a star of type j at velocity v∞ is then given by (Binney & Tremaine 1987, Eqn. 8-116) where and Rij = Ri + Rj. For cases where the outcome is largely insensitive to v∞, such as the collision of two main-sequence stars, we normalise over the Maxwellian velocity distribution to obtain For collisions involving red giants, where in general both the velocity and pericentre separation of the collision make a large difference to the outcome we first sample v∞ from the velocity distribution. There is a physical collision in timestep ∆t when X < R v ij ∆t, where X is uniformly chosen between 0 and 1, and R v ij comes from Eqn. 13. We obtain the collision radius as with Resc = 2G(Mi + Mj)/v 2 ∞ and R 2 rand = X/(πnjv∆t). In order to accelerate the calculation we bin the stars by mass, and treat all the main-sequence stars in a given bin as being equivalent. This is possible since the stellar radius and mass do not change very much over the main sequence for the low-mass stars that are relevant for our analysis. The main-sequence lifetime of each star is calculated individually, however, depending on its precise mass, and we treat giants individually since their masses and radii change significantly. Once the stars have evolved into compact remnants we again bin them, and have only a single compact object species of each type. For convenience the bin masses are the same as the masses for which stellar evolution tracks are calculated, as described below. End of HB Core He mass fraction reaches zero 5a Exhaustion of H and He burning b CO formation Code non-convergence during C burning c Nuclear reactions beyond C burning Table 1. Key points in the evolution of the stars and how they are defined

Stellar evolution
We follow the evolution of single stars using the stellar evolution code stars, originally written by Eggleton (1971), using the physical prescriptions described by Pols et al. (1995). Details of modifications to the code to make it run autonomously can be found in Church (2006), but are not particularly significant for this work. Our evolution tracks are computed at solar metallicity, with wind mass loss on the first giant branch and during core helium burning given by the formula of Kudritzki & Reimers (1978), with η = 0.4. On the asymptotic giant branch we follow the formula of Vassiliadis & Wood (1993). Our mass-loss law is not particularly suitable for massive stars, but our N -body tracks start at a time of ≈ 400 Myr, by which time all of the massive stars have already evolved into compact objects. Instead of calculating the evolution of each star individually, we make models of single stars at 76 different masses, chosen carefully to resolve the physics relevant to the questions that we are trying to answer. For each star of initial mass mi,0 we define a set of times tj(mi,0) which delimit key phases in the evolution of the star (end of core H burning, base of giant branch, etc.). The list of key times and their definitions is given in Table 1 To compute the state of a star of initial mass m0 at time t we then use the following procedure: (i) Identify the two tracks whose initial masses mi,0 and mi+1,0 bracket m0.
(ii) Interpolate the log of the key times, log tj, linearly 1 in mass between the two tracks to obtain a set of key times for the star in question, tj(m0).
(iii) Identify the pair of key times tj(m0) and tj+1(m0) that bracket the current time t and calculate the fractional time (iv) Identify the corresponding times in the two bracketing tracks; i.e.
and similarly for mi+1,0, and interpolate each stellar property (mass, core mass, log T eff , etc.) to that time.
CO species Initial mass range CO mass Table 2. Initial mass ranges and compact object masses for each compact object (CO) species in our stellar evolution model.
(v) Finally interpolate each stellar property linearly in mass between the two bracketing tracks at the corresponding values of τ .
This complicated procedure allows us to accurately reproduce stellar behaviour over the full range of masses using 76 tracks; a spot check of interpolated tracks against generated ones shows essentially identical HR diagrams with typical time-scale errors of ≈ 0.02 per cent. We choose this process in preference to using the fitting formulae of Hurley et al. (2000) because we wish to use additional derived quantities from the full models in the case of stellar mergers. To obtain K-band magnitudes we use the bolometric corrections and empirically calibrated T eff -colour relations of Lejeune et al. (1998).
At the end of its nuclear burning lifetime each remaining star is converted into a compact remnant. For performance reasons, we consider only four species of compact remnant; white dwarfs produced by normal single-star evolution (WDs), white dwarfs that have been produced by stellar collisions stripping the envelope of giant stars (cWDs), neutron stars (NSs), and black holes (BHs). Initial stellar mass ranges and compact object masses for these species are listed in Table 2.

Treatment of collisions and close encounters
We model collisions of both types of evolving star (MS and RG) with all types of compact object (cWD, WD, BH and NS), as well as collisions of MS stars with other MS stars and RGs. A summary can be found in Table 3. In the following subsections we describe and justify the collision prescriptions that we implement.

Main-sequence -main-sequence collisions
We assume that any collision between two main-sequence stars leads to merger with no loss of mass. Velocity dispersions in our model of the Galactic Centre are typically between 100 and 300 km s −1 , which is comfortably less than the surface escape velocities of main-sequence stars (≈ 600 km s −1 for the Sun). This treatment ignores the small fraction of mass that is expected to be lost in such encounters, and the high-velocity tail of potentially destructive encounters; neither are expected to be significant for our results.

Main-sequence -red giant collisions
For red giants, the analysis of the previous section is reversed. For a 100 R , 1 M giant, the surface escape velocity vesc,RG ≈ 60 km s −1 , less than the typical encounter velocity at the Galactic centre. The dwarf is also very much  Table 3. Effect of impactor (columns) on species collided with (rows). The table is not symmetric since a collision between two stars of different types will have different effects on each. See Section 3.2 for details more dense than the outer layers of the giant. As a result, the most likely outcome is for the dwarf to pass unscathed through the outer layers of the giant's atmosphere. Some of the giant's envelope is removed by the bow shock driven by the impacting main-sequence star, and some of that lost mass is accreted by the impactor. We parameterise these effects following the simulations presented by Bailey & Davies (1999) (BD99 hereafter). They studied collisions between a 2 M giant and a 1 M mainsequence star. 2 They present measurements of the mass lost from the system and the mass accreted by the impactor as a function of the encounter pericentre Rmin and velocity at infinity v∞ (figure 6 of BD99). We find that, when v∞ > vesc,RG, their results are well fit by the empirical formulae and where facc(r,ṽ) = max(−3, a acc r = Rmin/RRG, andṽ = V∞/Vesc,RG. The coefficients are given in Table 4 and the fits are plotted graphically in Figure 3. As we only have data for one giant mass, it is necessary to scale the results to other stellar masses. Physical reasoning suggests that, for homologous giant structures, the total mass and surface escape velocity are reasonable quantities to scale the mass loss with. It is also necessary to generalise the results from BD99 to different impactor masses. A more massive impactor will affect the atmosphere of a giant star from further away, and will both eject and accrete more mass as it passes through. In the absence of a comprehensive set of simulations we make two assumptions. Firstly, since the mass loss from a giant in an encounter is somewhat analogous to tidal stripping, instead of scaling Rmin with the giant 2 Whilst BD99 in fact studied a much larger range of impactor masses and two different target giants, they only presented detailed results for this one case, and the original simulation outputs are no longer available.

Mass loss (RGB
following Eggleton (1983). Secondly, we assume that the accretion of mass by the impactor as it passes through the giant's envelope largely follows the mechanism of Bondi & Hoyle (1944) and hence, following their equation 1, should be scaled with m 2 i . To test these assumptions we utilise the results of Davies et al. (1991) and Davies et al. (1992). They made SPH models of interactions of a 0.8 M giant with 0.4, 0.6 and 1.4 M point masses. Whilst they took V∞ = 10 km s −1 , appropriate to globular clusters but not to the majority of encounters at the Galactic Centre, the impactor velocities inside the giant envelope are still well in excess of the local sound speed, and so the relevant physics should not be very different. The rightmost panels of Figure 3 show their results scaled as described in the paragraph above; the curves for the three different impactor masses show rough agreement other than at very low masses where the simulations are underesolved. We consider that this validates our assumptions in the absence of further hydrodynamic simulations. Our final formulae for the mass lost and accreted by the impactor in giant-dwarf encounters are therefore   Davies et al. (1991) and Davies et al. (1992). In panel (e) we plot the mass accreted by the impactor, ∆Macc, divided by M env,RGB m 2 i , where M env,RGB is the mass of the giant envelope and m i is the mass of the point-mass impactor. The three sets of points lie on a single trend with min /R L,RGB , which suggests that this is the correct scaling of the accreted mass with mass of the giant and impactor, and closest approach distance R min . Similarly, in panel (f ) the mass lost by the giant, ∆M loss is divided by M env,RGB and again shows a single trend.
with f loss and facc given by Equations 19 and 21, and For each giant in our simulations, we accumulate the mass lost in successive collisions with main-sequence stars and compact objects. At the point when the total mass lost in collisions is greater than the remaining envelope mass the giant is converted into a compact object; i.e. even weak collisions that do not remove the entire envelope can shorten the giant's lifetime. The second-order effect of the change in stellar evolution induced by the lower envelope mass is ignored.
For the main-sequence star, we assume that the (relatively compact) main-sequence star acquires a hot, low-density halo of accreted material, as seen in hydrodynamic simulations of the encounters of compact stars with giants. We refer to such objects as fuzzballs. The amounts of mass accreted are typically rather small and we do not expect them to significantly affect the long-term evolution of the main-sequence stars, but we record the formation of these objects and consider their formation rates and expected properties further in Section 5.

Compact object -red giant collisions
For low-mass compact objects (WDs, cWDs and NSs) colliding with giants, our treatment follows that for MS-RG collisions in Section 3.2.2 above. The density contrast between a giant envelope and a main-sequence star is so large that the impactor can be treated as a point mass (and indeed is in the SPH simulations that we base our approach on). The same treatment then obviously applies for compact objects. Such collisions will also produce fuzzballs by accretion from the giant envelope. For collisions of black holes with giants, the mass ratio, and hence behaviour, is typically rather different. We base our treatment on Dale et al. (2009) and in particular their figure 11, where they present the envelope fraction remaining after collisions of a 1 M red giant with a 10 M black hole. We fit their results as where and, as before,ṽ = V∞/Vesc,RG,r = Rmin/RRG r L (q) r L (10) analogously with Equation 25, and the fit coefficients are given in Table 4. A comparison of the SPH results and fit is plotted in Fig. 4: whilst the fit is not perfect it provides a conservative estimate of the effect of such collisions on the population of red giants.
Most collisions that red giants suffer with other stars remove only part of their envelope. Dale et al. (2009) investigated the post-collision evolution of giants and found that it was necessary to remove almost all of the envelope to significantly affect the visibility of the giants. As with MS-RG collisions, for encounters that do not remove all of the envelope, we keep a running sum of the mass removed, which we add to the mass loss through stellar winds that the giant would have experienced in any case. When this total exceeds the mass of the envelope we take the giant to have evolved into a compact object. In this manner we are able to calculate a history for each giant without resorting to carrying out prohibitively expensive hydrodynamic or stellar evolution calculations on a star-by-star basis.

Main sequence -compact object collisions
We expect collisions between MS stars and compact objects to usually result in capture of the compact object by the star, following the same reasoning as in Section 3.2.1. Following capture, the accreted compact object will be more dense than the rest of the star, and hence will sink into the star's core. This will lead to the formation of a bright giant. If the compact object is a white dwarf the giant will be very similar to an ordinary star at the tip of the AGB, if it is a neutron star a Thorne-Żytkow object (Thorne & Zytkow 1975, 1977 will form, and if the compact object is a black hole then the result will be a quasistar (Begelman et al. 2006(Begelman et al. , 2008Volonteri & Begelman 2010). The details of the evolution of these objects are complicated and uncertain, but for us the most significant property is that they are expected to be short-lived. The MS star is thus effectively destroyed, and the compact object recovered as the final product of the evolution of the bright giant. It is possible that e.g. a neutron star could be converted into a black hole by this process, but it is rare enough that we do not need to take this pathway into account.

Collisions between pairs of giants
Collisions between pairs of giants could be as destructive as collisions between giants and main-sequence stars of similar mass, or perhaps even more so given the lower density contrast. However they are rare. Making a simple comparison between the rate of collision of two identical giants of mass M and radius Rg, and the rate of collision of a giant with a main-sequence star of equal mass M , we find that an enhancement of a factor of between two and four, depending on the velocity dispersion, from the ratio of number densities. But the giant phase is short, so the number density ratio is always small, particularly when lower-mass main-sequence stars are taken into account. Therefore we choose to ignore giant-giant collisions.

Collisions between pairs of compact objects
Collisions between pairs of compact objects are very rare owing to the small radii of the objects. Their only effect on our work would be if they removed a significant fraction of the compact objects, but this is not expected. We ignore them entirely (see e.g. Quinlan & Shapiro 1990;Benz et al. 1989, for an analysis of collisions between compact objects in galactic nuclei). ings. For any mass of main-sequence star the most common main-sequence impactor is from the lowest-mass bin, simply because those stars are most numerous, but the mass which they bring with them is relatively small. Plotting d 2 N coll /d ln M ,1 d ln M ,2 shows which mass impactors have the most effect on a given species. The plot shows that, for a given target star of mass M1, the fractional change in mass as a result of collisions is roughly equal for all impactors with mass M2 < M1. The effect of impactors with M2 > M1 is negligible. The peak effect is seen around 0.8 to 1 M , which is where the effects of collisions between main-sequence stars is having the most effect on the mass function: these are the highest-mass stars to be on the main sequence for the entire simulation duration. Figure 6 shows the final state of stars that initially had masses between 0.9375 M and 0.9875 M , binned by interaction parameter Γ . For a non-interacting population all of these stars would be giants. For an interacting superparticle, three factors affect the population. First, collisions with other main-sequence stars remove stars from the mass bin (grey long-dashed line). Secondly, the population is replenished at almost the same rate by the creation of stars in collisions. Many of these stars, however, do not manage to evolve onto the giant branch by the end of the simulation and are still on the main sequence at the end (green thin solid line). Finally, some of the initial population is destroyed by collisions with compact objects (brown dot-dashed line). It can be seen that the rate of destruction by compact objects is small compared to the rate of stellar mergers, but that the population removed by mergers is largely replenished. Indeed, if red giants at masses other than the bins around the turnoff mass are included (red dotted line) then the number of visible giants is actually slightly increased by collisions. Figure 7 gives an alternative view of how mergers and collisions with compact objects sculpt the stellar population. This figure shows the mass distributions of different types of stars and of the effects of stellar mergers and collisions at the end of the simulation. The thin dotted purple line shows the initial stellar population: for a non-interacting population this would be equal to the main sequence population (thin solid green line) below the turnoff, giants (thick solid red line) around the turnoff and compact objects (black solid line) above the turnoff. The short-dashed blue line shows the stars created by mergers, and the long-dashed grey lines the stars thus destroyed. The main effect of collisions is to convert low-mass stars into higher-mass stars, with the largest relative increase in population being around 2−3 M . Almost all of these stars have evolved by the end of the simulation, however, which leads to an increase in the number of white dwarfs descended from such stars (the bulge in the black line). At the present-day turnoff mass, more stars are destroyed than created by collisions, which leads to a significant reduction in the number of red giants around the turnoff; however, this is compensated by higher-mass red giants. Finally, it can be seen that the effect of collisions between stars and compact objects is negligible for all stellar masses. Figures 6 and 7 show that, even for the most interactive superparticles, the effects on the visible population of stars are modest. At no values of Γ do we see more than a 10% change to the mass function, and the total number of giants is barely changed by stellar encounters. Part of the reason for this is because the surface escape speeds of the stars are similar to the relative velocities of the stars in our simulations. This means that, with respect to the stellar velocity dispersion, the encounter rates are at a minimum. Hence, despite the very high number densities of stars, direct encounters between single stars are relatively inefficient at modifying the stellar population.  ible within the errors. Fainter, and probably older and more relaxed, stars in the Galactic NSC seem to show a steep profile (Schödel et al. see figure 12 in 2020 and figure 4 in Habibi et al. 2019). Given the slight collision-related changes in the different magnitude bins, we conclude that the luminosity function is not significantly affected by the collisions.

FORMATION AND EVOLUTION OF FUZZBALLS
One prediction of our results is that red giant stars in the Galactic Centre should undergo periodic collisions with objects that are much more compact than their envelopes: mainsequence stars, white dwarfs and black holes. These collisions are typically neither destructive of the giant nor, at the relative velocities found in the Galactic NSC, sufficiently dissipative for the compact star to be captured by the giant. Instead the giant suffers some mass loss, and some of that mass emerges bound to the compact star. Here we discuss the structure and observability of compact stars with tenuous accreted envelopes, which we refer to as fuzzballs.

Structure and evolution of fuzzballs
We assume that the captured gas (which we refer to as the envelope) reaches hydrostatic equilibrium within a freefall time-scale. Thereafter its properties are governed by the equations of stellar structure. We solve these equations to connect its observable surface properties -the radius R, effective temperature T eff and luminosity L -with the mass of the compact star Mc, the mass of the envelope Menv and the thermal time-scale of the fuzzball envelope, τ th,env . Our analysis is based on Section 7.3.4 of Eldridge & Tout (2019), with appropriate modifications. We make the assumption that the envelope is unstable to convective motions, and hence isentropic with γ = 5/3: this allows us to write the pressure P in terms of the density ρ as where K is a constant to be determined. We further assume that Menv Mc, as seen in SPH simulations of compact stars colliding with giant envelopes: hence we can take the enclosed mass m = Mc in the equation of hydrostatic equilibrium, which becomes dP dr where r is the radial co-ordinate, G is the gravitational constant and we have substituted Eq. 29. Rearranging and integrating with respect to r gives the density as where ρs is the surface density. This allows us to obtain the envelope mass as where Rc is the radius of the base of the envelope. The value of Rc is largely immaterial because, close to the core, the enclosed mass scales as r 3/2 : hence the inner envelope's mass is too small to have much effect on the fuzzball's bulk properties. It remains to define boundary conditions. In addition to the isentropic relation we apply an ideal gas equation of state at the surface which allows us to obtain an expression for K in terms of the surface density ρs: where R = 8.314 kJ kg −1 K −1 is the gas constant and the mean molecular weight µ ≈ 0.7 for a typical population I star envelope. Finally, we apply a grey surface boundary condition for the surface pressure Ps and opacity κ: We first choose R and T eff , and then solve this boundary condition iteratively using a bicubic spline fit to the opacities of Ferguson et al. (2005) to obtain ρs. The envelope mass is then obtained by integrating Eq. 32. Finally, we obtain the internal energy of the envelope as where we have again used the ideal gas law, and the gravitational potential energy of the envelope by The thermal time-scale for the envelope is then given by where the luminosity L is given by the usual surface boundary condition L = 4πR 2 σ sb T 4 eff and σ sb is the Stefan-Boltzmann constant. Our calculation of the time-scale implicitly neglects any luminosity provided by the central compact star or by nuclear burning at the base of the envelope; the sole source of energy is the contraction of the envelope.

Results
We start by considering the properties of fuzzballs with white dwarf cores, since these should be one of the more common cases. Figure 9 shows the envelope mass Menv as a function of T eff and R. Black, red, and white dashed contour lines show envelope masses of 10 −3 , 3 × 10 −3 and 10 −2 times the white dwarf mass; these represent the upper end of the distribution of envelope masses found in collisions in our models. A given fuzzball has a fixed envelope mass and hence will follow a contour line, radiating away its internal energy as it shrinks. The structure in the diagram is largely set by the complex dependence of the surface opacity on temperature and density. It shows three branches in the figure -the branch at around 1500 K corresponds very roughly to the classical Hayashi track. Our somewhat cavalier treatment of the inner boundary condition frees us from the Hayashi limit, however, and permits an additional, cooler, set of solutions. These intersect nicely with the properties of the G2 cloud inferred by Witzel et al. (2014) (T eff,G2 = 560 K and RG2 = 570 R ). We infer that if G2 is a WD-cored fuzzball it must have an envelope mass of 1.5 × 10 −3 M .
While the surface properties of our model of G2 as a WDcore fuzzball are compatible with the observations, the model is excluded by the inferred thermal time-scale. Figure 10 shows the thermal time-scale as a function of T eff and R, with the star again marking the location of G2; our model at  Figure 10. Thermal time-scale τ env,th as a function of T eff and R for fuzzballs with a white dwarf core. Black regions of the plot have total envelope energies greater than zero; i.e. the envelope is not bound. The white, red and black dashed contour lines show envelope mass ratios Menv/Mc of 10 −3 , 3 × 10 −3 and 10 −2 respectively. The black star shows the location of G2 as measured by Witzel et al. (2014). this point in the diagram has a thermal time-scale of 4.4 yr. This implies that G2 would have shown significant evolution of its surface properties over the last few years, which has not been observed, and that it was still in close proximity to the site of the close encounter. In addition, the total mass is probably insufficient for G2 to have survived tidal disruption during its close passage by Sgr A*. Figure 11 shows how the envelope mass and thermal timescale of fuzzball models for G2 vary as a function of Mc. As can be seen, models start to become compatible with observational constraints on the thermal time-scale from Mc ≈ 2 M . This rules out WDs and probably all but the most massive NSs as cores of a putative G2 fuzzball. Black holes are still compatible, as are intermediate-mass main-sequencestar cores, but in the latter case more careful modelling would be required to take into account the effect of the star's own luminosity.
For completeness we include plots of the envelope mass ( Figure 12) and thermal time-scale ( Figure 13) for fuzzballs with central 5 M black holes. In this case a larger ratio of envelope to core mass is required to match the G2 surface boundary conditions. As a result of the larger envelope masses the thermal time-scales are greatly increased. A BH fuzzball with the surface properties of G2 requires an envelope of 0.34 M and has a thermal time-scale of 9400 yr, compatible with the properties of G2.

Fuzzball formation rate
Encounters between giant and compact stars are fairly common in our model. We look at the average number of encounters in the whole central projected 2 pc during the last Gyr of the simulation, after the NSC is fully assembled. We find a total of 2100 WD-RG encounters, 65 NS-RG encounters and 15 BH-RG encounters. Encounters with low-mass mainsequence stars are even more common, mostly with MS stars in the lowest mass bin. However, most of the encounters lead to very little mass accretion. Figure 14 shows the mass accreted by white dwarfs in encounters with giants in the final Gyr of our model. The vast majority of encounters are fast and wide, since the the target presented by a giant is larger at larger radii, and the velocity dispersion in the central cluster is significantly larger than the surface escape speeds of the giants. Only the rare, slow, close encounters accrete sufficient mass to form G2-like clouds, with a total rate of about 60 per Gyr. Given the very short lifetimes of WD fuzzballs it is clear that G2 cannot be explained in this way.

Frequency of harassment of giants
A corollary to our previous statement is that the fraction of giant stars that have had an encounter with a compact star since leaving the main sequence is potentially rather large. Figure 15 shows the fraction of giant stars that have had at least one encounter with a compact star, split into bins by apparent K-band magnitude. The figure shows that, as stars evolve up the first giant branch, the collision fraction increases steadily up to around one quarter of giants by the top of the RGB. By the tip of the AGB about 40% of giants will have had an encounter. It is beyond the scope of this paper to study in detail the effect of these encounters on the giants, but they will include a reduction in the giant lifetime, and a temporary increase in luminosity from deposition of thermal energy by the passing stars.
NSCs are among the densest stellar systems in the Universe and, as such, their stars should undergo frequent interactions. Stellar collisions have been addressed in a number of studies, relying on analytical and observational based models of the Galactic nucleus (Davies et al. 1998;Alexander 1999;Bailey & Davies 1999;Dale et al. 2009). In this work we tackle this problem using the phase-space distribution of stars in a Milky Way-like NSC, built through the infall and merger of GC-like massive clusters (see e.g. Tsatsi et al. 2017;Abbate et al. 2018). The NSC formed in this N -body simulations show morphological (e.g. flattening) and kinematic (e.g. rotation, velocity dispersion) properties comparable to those of the Galactic NSC. After converting each N -body particle into a stellar population, we applied a stellar evolution routine and collision prescriptions derived from the literature to study the cumulative effects of the stellar encounter history at the Galactic Centre. We took into consideration collisions of main sequence and red giant stars with compact objects (WDs, BHs and NSs), as well as collisions of main sequence stars with other main sequence stars and red giants, while we neglected rare events such as collisions between red giant stars and between pairs of compact objects.
The most common encounters in our models are between pairs of main-sequence stars. Stars of a given mass are only significantly affected by collisions with less massive impactors, and the peak effect on the mass function is observed for stars with masses between 0.8 and 1 M , i.e. for the highest-mass main sequence stars whose lifetime is not less than the simulation length. While the rate of stellar mergers is greater than the rate of destruction by compact objects, stars are removed and created in mergers almost at the same rate, so at the present day the red giants whose progenitors were destroyed by collisions are replaced by higher-mass red giants born in mergers.
Furthermore, we find that the rate of collisions between stars and compact objects is so low that their effect on the mass function is negligible. The mass function of the visible stars in the Galactic NSC should, therefore, only be slightly affected by stellar encounters: we see changes smaller than 10%. The final core-like density profile of our NSC is not significantly altered by the collisions. Stellar encounters marginally affect the more luminous stars but leave the luminosity function of the stars in the NSC mostly unchanged.
Our NSC formation scenario leads to a central core at the end of the N -body simulation. We, therefore, do not probe the role that collisions might play in disrupting a cusp (Davies et al. 1998;Alexander 1999;Bailey & Davies 1999;Dale et al. 2009;Amaro-Seoane & Chen 2014;Amaro-Seoane et al. 2020). Analysis of a static cuspy NSC model with properties similar to those of the Galactic nucleus suggest that, in this case, stellar encounters could have a non negligible impact on the luminosity function and spatial distribution of the observable stars. However, a more sophisticated analysis of a time-evolved MW-like NSC showing a central stellar cusp Schödel et al. 2020;Habibi et al. 2019) will be necessary to test this scenario in detail.
We follow the evolution of fuzzballs -compact stars that have accreted a tenuous envelope in collisions with red giants -and conclude that the properties of the enigmatic G2 object observed at the Galactic centre could be explained as the result of an encounter between a red giant and a black hole. On the other hand, we exclude the possibility that G2 could be a WD fuzzball, given the short lifetime of these objects. Our collision rates make the BH-fuzzball model for G2 very unlikely, as our model has a core number density of black holes in the central 0.5 pc of only 1400 BHs pc −3 . The collision rate would be much higher if a cusp of black holes is present in the central NSC (see e.g. Miralda-Escudé & Gould 2000;Freitag et al. 2006;Alexander & Hopman 2009;Hailey et al. 2018;Davies & Lin 2020, for theoretical and observational studies on the presence of massive stars and SBHs at the Galactic centre). The number of black holes in the galactic centre could also be enhanced by local star formation. The observed population of massive young stars seen today will likely produce up to 100 black holes. If this mode of star formation occurs every 100 Myr or so (Amaro-Seoane & Chen 2014), then in total, this mode of star formation will add roughly 10 4 black holes in the central regions thus enhancing the population by a factor of roughly ten. In turn if some of these black holes are more massive (around 20 M ), then the thermal relaxation times for any fuzzball produced will be commensurately longer (as shown in Figure 9). Thus with a population of black holes enhanced in this way, the chance of seeing G2 today could be a few per cent or higher.
Another possibility is fuzzballs produced by collisions of giants with intermediate-mass stars. In this case the luminosity provided by the central star might be expected to compensate for the loss of energy from the photosphere of the tenuous envelope, since G2's luminosity is roughly the same of a 2 M main-sequence star. However, it is by no means clear that the structure of a central core and tenuous envelope would persist over multiple thermal timescales, and modelling the evolution of such objects would require a solution of the full set of equations of stellar evolution which is beyond the scope of this paper.
Another possible route to form G2 is collisions between main-sequence stars. Collisions deposit large quantities of kinetic energy which produce a distended star, possibly with surface properties similar to G2. We see seven MS-MS collisions per Myr in the final 10 9 yr of simulations, which suggests that this is an interesting avenue to explore in further work.
In our model, about 40% of the giants undergo at least one encounter by the end of their AGB phase. Such encounters may shorten their lifetimes, and temporarily increase their luminosity due to the energy deposited by the interacting star. Further SPH simulations and longer-term follow-up with stellar evolution calculations will be necessary to follow the details of such events, as well as to to better understand the formation and evolution of fuzzballs.
Our results might suggest that stellar encounters are unlikely to be of any significance to the stellar population at the Galactic Centre, and that further work on this complex edifice of encounter simulations promises to be of little interest. However, as aluded to above, various assumptions that we have made conspire to make this a very pessimistic scenario. The core-like structure of our NSC reduces the central density of the system and the presence of the supermassive black hole throughout the evolution in effect enforces the high velocity dispersion that reduces the stellar encounter rate. Secondly, we inject the stars into the simulation when the cluster enters the field, but ignore their previous evolution. Globular clusters undergo dynamical encounters which not only modify their stellar populations by collisions and coagulations, but also cause segregation of the more massive stars into the centres of the clusters. Failing to include these effects means that our results are inherently a pessimistic estimate of the effects on the stellar population from encounters.
We plan to refine our models to investigate the likely effects of some of these assumptions in a subsequent paper. In particular, we will combine our NSC formation simulations with realistic GC models (e.g. obtained with MOCCA Giersz et al. 2008, 2013 to provide more accurate predictions on the characteristics of the stellar populations in the Galactic nucleus. We will also explore the collisional evolution of a cusp-like NSC with the same properties of the Galactic NSC to once and for all check if the collisions can flatten the density profile.

CONCLUSIONS
In this work we attempt to quantify the effect that close encounters between stars has had on the stellar populations in the Milky Way nuclear star cluster (NSC). NSCs should witness significant rates of stellar encounters owing to their high number densities. We model the stellar evolution and collision history of a Milky Way-like NSC built through the infall of stellar systems similar to GCs, basing our simulations on N -body models of NSC formation through the merger of multiple massive stellar clusters. We find that the most common encounters to be expected in the NSC are collisions between pairs of main-sequence stars.
Stellar collisions and mergers have little effect on the observable stellar populations, both in terms of their mass and luminosity functions. There are three reasons for this. First, our N -body models form a core rather than a cusp, which limits the central number density that we reach and hence the collision rate. Second, the presence of a supermassive black hole from the beginning of the simulation leads to a high velocity dispersion, similar to the surface escape speeds of the interacting stars. This coincidence minimises the rate of stellar encounters. Third, while collisions between main-sequence stars can deplete the population of bright giants by accelerating the evolution of their progenitors, it can also enhance the population by accelerating the evolution of lower-mass mainsequence stars. For modest collision rates these effects cancel and the population of visible giants is almost unaffected.
Additionally we investigate the possibility that a mainsequence star or compact remnant could accrete a tenuous envelope when passing through the outer layers of a red giant. We test whether such objects, which we refer to as fuzzballs, could explain G2, a cloud observed to have made a close passage past the MW SMBH and whose nature is still debated. We find that fuzzballs where the accretor is a white dwarf or low-mass main-sequence star are ruled out on the grounds of their thermal time-scales being too short. A blackhole fuzzball has a compatible thermal timescale but in our model the formation rates of such objects are too low to convincingly explain G2.
We plan to overcome the limitations of our N -body simulations and analysis by considering more realistic GC models obtained e.g. in Monte Carlo simulations (e.g. MOCCA Giersz et al. 2008Giersz et al. , 2013. We will also consider models in which the underlying density profile of the NSC is steeper, which could lead towards a higher rate of collisions and to a more significant effect on the stellar populations at the Galactic centre.