Stellar Collisions in the Galactic Center: Massive Stars, Collision Remnants, and Missing Red Giants

Like most galaxies, the Milky Way harbors a supermassive black hole (SMBH) at its center (e.g., Genzel et al. 2003; Ferrarese & Ford 2005; Ghez et al. 2005; Kormendy & Ho 2013). A dense region of stars and stellar remnants, known as the nuclear star cluster, surrounds the SMBH (e.g., Schödel et al. 2003; Ghez et al. 2005, 2008; Gillessen et al. 2009, 2017). The proximity of the Milky Way's Galactic center (GC) presents a unique opportunity to observe the consequences of this dense, dynamic environment. One such consequence is that direct collisions can occur between stars (e.g., Dale & Davies 2006; Dale et al. 2009; Rose et al. 2020, 2022; Mastrobuono-Battisti et al. 2021). As we show below, any 1 M_⊙ star within about 0.1 pc of the SMBH experiences at least one collision over its lifetime. At 0.1 pc from the SMBH, the collision timescale is comparable to the main-sequence lifetime of a 1 M_⊙ star. Collisions are therefore essential to understanding the properties and demographics of the nuclear star cluster.

A number of studies have explored the effects of this dynamical process on the nuclear star cluster. For example, stellar collisions can modify the stellar density profile near the SMBH, where they deplete the stellar population (e.g., Duncan & Shapiro 1983; Murphy et al. 1991; David et al. 1987a, 1987b; Rauch 1999; Freitag & Benz 2002). Collisions have also been invoked to explain the dearth of 1–3 M_⊙ red giants (RGs) in the inner 0.1 pc of the GC, a persistent observational puzzle (e.g., Genzel et al. 1996; Bailey & Davies 1999; Buchholz et al. 2009; Do et al. 2009; Baumgardt et al. 2018; Gallego-Cano et al. 2018; Schödel et al. 2018; Habibi et al. 2019). Genzel et al. (1996) first suggested that RG collisions with main-sequence stars may destroy the RGs. However, theoretical studies diverge on whether such a collision can divest a RG of its envelope, with more recent work indicating insufficient mass loss post-collision (e.g., Alexander 1999; Bailey & Davies 1999; Dale et al. 2009; see the latter for discussion). Other mechanisms, such as black hole–RG collisions, RG encounters with stellar binaries, and a nuclear jet from the SMBH during its active phase, can only partly explain the apparent underabundance of RGs within 0.1 pc of the SMBH (e.g., Davies et al. 1998; Dale et al. 2009; Zajaček et al. 2020), though RG interactions with a gaseous disk in the GC represent a viable alternative explanation for their underabudance (e.g., Amaro-Seoane & Chen 2014; Amaro-Seoane et al. 2020). The cause of the missing RGs remains an open question, one which main-sequence stellar collisions may yet address (see Mastrobuono-Battisti et al. 2021, and our discussion in Section 7).

Collisions and mergers in the nuclear star cluster can also yield a wide range of observables (e.g., Amaro Seoane 2023). Stellar collisions may be associated with several electromagnetic signatures, including variability of active galactic nuclei (e.g., Murphy et al. 1991; Torricelli-Ciamponi et al. 2000; Freitag & Benz 2002), the presence of blue stragglers (e.g., Sills et al. 1997, 2001; Lombardi et al. 2002), and supernova-like explosions (Dale & Davies 2006; Balberg et al. 2013). Furthermore, stellar mergers represent a possible explanation for the origin of G objects, stellar objects that are enshrouded by gas and dust observed in the nuclear star cluster (albeit secular binary mergers rather than impulsive collisions; Antonini et al. 2010; Witzel et al. 2014, 2017; Stephan et al. 2016, 2019; Ciurlo et al. 2020). Lastly, collisions between compact objects can release gravitational-wave emission, potentially detectable by the LIGO-Virgo-KARAGA collaboration (e.g., O'Leary et al. 2009; Hoang et al. 2018, 2020; Arca Sedda 2020).

Also with implications for gravitational-wave sources, repeated collisions between stars and black holes (BHs) can cause the latter to grow in mass, producing more massive BHs than predicted by stellar-evolution models (see Rose et al. 2022 for BH mass predictions; see also Heger et al. 2003; Spera & Mapelli 2017; Woosley 2017; Limongi & Chieffi 2018; Belczynski et al. 2020; Renzo et al. 2020). However, for the BHs to grow significantly, there must be a sufficiently large population of stars to sustain the collisions. Star–star collisions may deplete this supply, especially in the region closest to the SMBH, ≲0.01 pc, where BH–star collisions can act most efficiently to produce larger BHs (Rose et al. 2022). Therefore, understanding the impact of collisions on the stellar population has important implications for the compact-object–star interactions that give rise to massive BHs and interesting electromagnetic signatures (for the latter, see Kremer et al. 2022).

In this study, we consider the effect of stellar collisions on the masses and survival of stars in the GC, with implications for the topics described above. We show that collisions between low-mass stars can give rise to a population of young-seeming, massive stars in the vicinity of the SMBH as well as peculiar low-mass stars that have been divested of their outer layers. Additionally, our results suggest that stellar collisions may explain observations of dust-enshrouded gaseous stellar objects and a missing RG cusp in the GC. This paper is organized as follows. In Section 2, we describe our methodology. Many studies rely on smooth particle hydrodynamics (SPH) to determine the outcome of a collision, a complex problem that depends on assumptions about the stellar structure, impact velocity, and mass ratio of the stars (e.g., Benz & Hills 1987; Lai et al. 1993; Rauch 1999; Fregeau & Rasio 2007; Kremer et al. 2020; Rodriguez et al. 2022). Given these uncertainties, we test different treatments of the mass loss and collision outcomes, documented in Section 2.3 and compare their results in Section 3. Sections 4 and 5 present and discuss our findings in the context of the stellar demographics and mass function. Sections 6 and 7 discuss the implications for the G-object and RG populations, respectively. In Section 8, we summarize our findings.

2.1. Initial Conditions

We follow a sample population of 1000 stars embedded in the nuclear star cluster. We draw the initial masses of our sample population using a Kroupa initial mass function (IMF). We limit the stars to masses between 0.5 and 100 M_⊙. However, in practice, the most massive stars in our simulations tend to be around 30 M_⊙ because of the steep IMF. The semimajor axes of the stars are drawn from a uniform distribution in log distance, i.e., one in which ${dN}/d(\mathrm{log}r)$ is constant. This distribution ensures adequate sampling at all distances from the SMBH, necessary to build comprehensive understanding of stellar collisions in the GC.

We draw the orbital eccentricities of the stars from a thermal distribution. While the young stellar population has lower eccentricities on average than a thermal distribution (e.g., Yelda et al. 2014; von Fellenberg et al. 2022), the nuclear star cluster is dominated by an old population of stars, with ages on the order of a billion years (e.g., Chen et al. 2023). We therefore assume that the eccentricity distribution for most of the stars has been thermalized; the stars have had sufficient time to interact with one another and exchange energy. Additionally, the S-star cluster exhibits a thermal eccentricity distribution, which may have been achieved through relaxation processes within a 10⁷ year time frame (e.g., Gillessen et al. 2017; Generozov & Madigan 2020). We explore deviations from an initial thermal eccentricity distribution in Appendix B and find that it has only a minor impact on our final results.

Each star in our sample population is influenced by the surrounding stellar cluster, which we treat as a reservoir. For simplicity, we assume that the surrounding stars have an average mass of 1 M_⊙, which remains constant over the duration of our simulation. This value is selected as a good approximation for the average stellar mass in a star cluster with a top-heavy IMF, as has been observed in the nuclear star cluster (e.g., Paumard et al. 2006; Bartko et al. 2010; Lu et al. 2009, 2013; Alexander & Pfuhl 2014). ⁵ We note that if the IMF is not top heavy, perhaps with an average stellar mass closer to 0.33 M_⊙, then we are overestimating the amount of mass added to a star in a collision. However, over time, collisions will also change the average mass within the population. Heuristically, we expect that the average stellar mass doubles every collision timescale. The average stellar mass is therefore a function of both time and distance from the SMBH. Outside of 0.1 pc, it should be approximately constant.

We describe the density of the cluster as a function of r_•, the distance from the SMBH, using a power law with the form

$$\begin{eqnarray}&&\rho ({r}_{\bullet })={\rho }_{0}{\left(\displaystyle \frac{{r}_{\bullet }}{{r}_{0}}\right)}^{-\alpha },\end{eqnarray} \tag{ 1 }$$

where α determines the slope of the density profile. For the inner parsec of the GC, the density is normalized using the observationally constrained values ρ₀ = 1.35 × 10⁶ M_⊙ pc⁻³ at r₀ = 0.25 pc (Genzel et al. 2010). For a uniform population of solar-mass stars, the number density becomes

$$\begin{eqnarray}&&n({r}_{\bullet })=\displaystyle \frac{\rho ({r}_{\bullet })}{1\,{M}_{\odot }}.\end{eqnarray} \tag{ 2 }$$

Unless otherwise noted, we assume that the surrounding star cluster is dynamically relaxed, resulting in a cuspy Bahcall–Wolf profile (i.e., α = 1.75 Bahcall & Wolf 1976; Alexander & Hopman 2009; Keshet et al. 2009; Bar-Or et al. 2013; Aharon & Perets 2016). The true density profile in the GC may be shallower, with an index closer to 1.5 or 1.25 (e.g., Buchholz et al. 2009; Do et al. 2009; Bartko et al. 2010; Schödel et al. 2014, 2018, 2020; Gallego-Cano et al. 2018, 2020; Linial & Sari 2022). Appendix D shows results that use α = 1.25. The most realistic picture of the GC may lie somewhere between the α = 1.25 and 1.75 results.

The mass of the SMBH, M_•, sets the velocity dispersion of the cluster:

$$\begin{eqnarray}&&\sigma ({r}_{\bullet })=\sqrt{\displaystyle \frac{{{GM}}_{\bullet }}{{r}_{\bullet }(1+\alpha )}},\end{eqnarray} \tag{ 3 }$$

where α is the slope of the density profile (Alexander 1999; Alexander & Pfuhl 2014), and σ is approximately the orbital velocity. Together, the velocity dispersion and number density govern key physical processes like collisions because these properties determine the frequency of interactions between a star and its surrounding objects.

2.2. Direct Collisions

A star in the GC is expected to experience a collision over a characteristic timescale, t_coll. This timescale can be estimated using an (n σ A)⁻¹ rate calculation, where A is the cross section of interaction and n and σ are determined by Equations (1) and (3). The collision timescale also has a weak dependence on eccentricity:

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{coll}}^{-1} & = & \pi n({a}_{\bullet })\sigma ({a}_{\bullet })\\ & & \times \left({f}_{1}({e}_{\bullet }){r}_{c}^{2}+{f}_{2}({e}_{\bullet }){r}_{c}\displaystyle \frac{2G({M}_{\odot }+M)}{\sigma {\left({a}_{\bullet }\right)}^{2}}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where f₁(e_•) and f₂(e_•) are given by Rose et al. (2020), G is the gravitational constant, and r_c is the sum of the radii of the colliding stars.

We plot the collision timescale for different density profiles, spanning α = 1.25 (dashed) to α = 1.75 (solid), in red in the bottom-left panel of Figure 1. In the upper-left panel, we plot the number of collisions expected per star as a function of distance from the SMBH. The number of collisions per star are also plotted in red, with dashed and solid lines corresponding to the α = 1.25 and α = 1.75 density profiles, respectively. These curves assume a uniform population of solar-mass stars. On the plot, we use gray vertical lines to mark the distance within which each star will experience at least one collision. This distance corresponds to where the collision timescale is equivalent to the main-sequence lifetime of a solar-mass star. Outside of this distance, the number of collisions per star is less than one, indicating the probability that a single star will experience a collision or, equivalently, the fraction of the population that collides. In the upper panel, we also show the velocity dispersion as a function of distance from the SMBH in units of escape velocity from the star. Since we take the velocity dispersion to be the relative velocity of the two stars as they collide, this gray curve is a measure of how energetic each collision is. Where σ is greater than the escape speed from a solar-mass star, we expect the collisions to be destructive, leading to mass loss. We discuss collision outcomes in greater detail in the next section.

Zoom In Zoom Out Reset image size

In our simulations, we account for collisions using a statistical approach following Rose et al. (2022). Over a time step, Δt, the probability that a given BH will experience a collision is approximately Δt/t_coll. We draw a random number between 0 and 1. If that number is less than or equal to the probability, we assume a collision has occurred. We adjust the mass of the star accordingly and repeat this process until the simulation has reached the desired total time, 10 Gyr. As in Rose et al. (2022), we use a time step of 1 × 10⁶ yr.

2.3. Collision Outcomes

2.3.1. General Overview

The outcome of a direct collision between two stars depends on many factors, including the velocity dispersion, the impact parameter, and assumptions about the stellar structure (e.g., Lai et al. 1993; Freitag & Benz 2005). Given the uncertain outcome of stellar collisions, especially at high velocities, we explore a range of possibilities, varying our treatment of the mass-loss and merger products. Broadly, there are two possible outcomes of a collision:

• 1.  
  Merger: in this case, the star in question gains a solar mass. We reiterate that in our toy model, the surrounding cluster contains only 1 M_⊙ stars, while the masses of our sample are drawn from a Kroupa IMF. From the merger product, some mass may be ejected depending on the nature of the collision. We estimate the final mass, M_f , of the star as M_f = M_i + 1 − f_ml (M_i + 1), or M_f = (1 − f_ml )(M_i + 1), where M_i is the stellar mass before the collision took place and f_ml is the fractional mass loss caused by the collision, discussed in greater detail below.
• 2.  
  Mass-loss only: in the second case, there is no merger—the primary star fails to capture the collider—but the star still loses some fraction of its mass following the impact. Its final mass equals (1 − f_ml )M_i .

In the following subsections, we detail different physically motivated approaches to estimating the fractional mass loss. Following a collision, we recalculate the radius of the star using its updated mass. The thermal timescale is short enough (10⁷ yr or shorter) for us to take the radius of a main-sequence star, except for the low-mass (<1 M_⊙) stars. At 0.01 pc, approximately the innermost radius where we expect mergers to occur (see Figure 2), the collision timescale is more than 10 times longer than a thermal timescale of 10⁷ yr. As the collision timescale is inversely proportional to R² and the Kelvin–Helmholtz timescale depends on 1/R, the radius of the merged star must be inflated to 10 times its initial value for these timescales to be comparable. Our default merger conditions require the relative velocity between the colliders to be less than the escape velocity from the star (see Section 2.3.3). For this reason, the degree of inflation post-collision should be not so extreme, and most of the star's mass should remain concentrated toward the center (e.g., Lai et al. 1993). We therefore conclude that the stars will have time to recollapse before they experience another collision. For the low-mass stars, however, the thermal timescale may be longer than the collision time, meaning the stars remain distended when the next collision occurs. The larger radius than accounted for in our model may hasten their destruction.

Zoom In Zoom Out Reset image size

Irrespective of the specific prescription adopted for f_ml , we can make predictions about where in the GC certain outcomes are likely to occur. Figure 1, right panel, illustrates the parameter space in which different outcomes are expected. In particular, the right plot shows the collision timescale as a function of the star's mass at different distances from the SMBH. Each collision timescale curve is color-coded according to distance from the SMBH, as indicated by the color bar on the left side of the plot. We also include the main-sequence lifetime as a function of stellar mass (solid black line) and the simulation time of 10 Gyr (dashed black line). The region highlighted in red represents a regime in which most stars will evolve off the main sequence before a collision can occur. On the other hand, when the collision timescale is less than the main-sequence lifetime, stars may experience one or more collisions before evolving off the main sequence. The regime highlighted in dark blue in the right-hand plot represents the scenario where the dispersion velocity on the onset of the collision is larger than the escape velocity of the two colliding stars (i.e., σ > v_esc). Thus, we expect collisions to result in significant mass loss from the stars and no mergers to occur. However, in the small light blue wedge, the dispersion velocity is smaller than the escape velocity, and collisions may result in mergers and mass growth. As mentioned above, the details of the mass loss and growth are uncertain, and we consider several scenarios below.

2.3.2. Energy Ratio Estimate of the Mass Loss

The fractional mass loss can be approximated as the ratio of the kinetic energy to the binding energy of the two colliding stars:

$$\begin{eqnarray}&&{f}_{{ml}}=\displaystyle \frac{\mu {\sigma }^{2}}{{{GM}}_{i}^{2}/{R}_{i}+{{GM}}_{\odot }^{2}/{R}_{\odot }},\end{eqnarray} \tag{ 5 }$$

where μ is the reduced mass of the two colliding stars and M_i and R_i are the mass and radius of the star from our sample population. This equation is roughly consistent with equations presented in the qualitative overview of collisions in Lai et al. (1993), with some adjustments to make the fractional mass-loss unitless. The aforementioned study includes a parameter β in the denominator of f_ml , a variable that accounts for the gas pressure in the star. While β is order unity for solar-mass stars, radiation pressure becomes important for more massive stars; it is easier to unbind material from them during a collision (Lai et al. 1993). The main-sequence lifetime for these massive stars, however, is short compared to the collision timescale. Their low probability of experiencing a collision means that our assumption that β = 1 does not significantly impact our results.

Stars are not uniform-density spheres: in fact, most of a star's mass is concentrated near its center, leading to a higher binding energy than the estimate in the denominator of Equation (5) (e.g., Christensen-Dalsgaard et al. 1996). To illustrate this point, we plot the mass enclosed as a function of radius within the Sun based on Christensen-Dalsgaard et al. (1996) in Figure 10 (see Appendix A). A spherical object with this mass profile has binding energy 1.61GM²/R. For a slightly more accurate estimate of the fractional mass loss, we include this prefactor, 1.61, in Equation (5). In Figure 2, we plot the fractional mass loss with and without the prefactor in green and blue, respectively, for a collision between two 1 M_⊙ stars. The fractional mass loss decreases with distance from the SMBH because it depends on the velocity dispersion.

In this prescription, if $\sigma \lt \sqrt{2{GM}/R}$, roughly the escape velocity from the stars, we assume that a merger occurs (see Figure 1). Otherwise, we calculate the final mass of the star as M_post = M(1 − f_ml ). This prescription overpredicts the amount of mass lost for reasons discussed in Section 2.3.4.

2.3.3. Toy Core+Envelope Model

The purpose of this toy model is to provide physical insight into the outcome of stellar collisions. Most of a star's mass is concentrated at its center. As a result, the escape velocity from the surface of the star may be much lower than that from the core, a subtlety not fully encapsulated by the approximations in Section 2.3.2. Here, we construct a more nuanced toy model based on the internal structure of the Sun (e.g., Christensen-Dalsgaard et al. 1996). Based on the mass enclosed, M_enc, within a radius r in the Sun, we calculate the escape velocity from a sphere with equivalent mass and radius, ${v}_{\mathrm{esc}}=\sqrt{2{{GM}}_{\mathrm{enc}}/r}$, as a function of r. We determine that this v_esc reaches a maximum of about 800 km s⁻¹ at R = 0.33 R_⊙ (see Figure 10 and Appendix A). This radius encloses almost exactly two-thirds of the Sun's mass. We define this region as the "core," though we stress that it is not representative of the true core of the Sun; in our simple model, there is a tightly bound inner part of the star and a more loosely bound outer "envelope."

We calculate the fractional mass loss as follows. We begin by calculating the ratio of the kinetic energy to the binding energy of the stars as in Equation (5). We then take the maximum between the corresponding change in mass and the mass of the "envelope," 0.33 M_⊙. If the kinetic energy is greater than the binding energy of the envelope, we assume the remaining kinetic energy can go into ejecting material from the more tightly bound "core." In other words, if the collision has enough energy to unbind more than the outer third of the star's mass, the total mass lost equals

$$\begin{eqnarray}&&{M}_{\mathrm{lost}}={M}_{\mathrm{env}}+\displaystyle \frac{\mu {\sigma }^{2}/2-{U}_{\mathrm{env}}}{{{GM}}_{\mathrm{core}}^{2}/{R}_{\mathrm{core}}},\end{eqnarray} \tag{ 6 }$$

for a collision between two 1 M_⊙ stars. As noted above, in this equation M_core is about 0.66 M_⊙ and R_core is about 0.33 R_⊙. In Figure 2, we plot the fractional mass loss in orange.

The main advantage of this prescription is that it can be expanded to account for an impact parameter. From a geometric standpoint, head-on collisions are unlikely. Collisions with a nonzero impact parameter result in less mass loss because only the outer envelopes of the stars interact, while most of the mass is concentrated in the core. Accounting for the impact parameter, b, is therefore crucial in this work. We draw an impact parameter between zero (a head-on collision) and the sum of the two stellar radii in question, R + R_⊙ (a grazing collision) using a probability density that goes as ${b}^{2}/{(R+{R}_{\odot })}^{2}$.

In our simple model, during a collision with ${R}_{\mathrm{core}}\lt b$, only the envelopes interact and there cannot be any additional mass loss from the core. On the other hand, any collision with an impact parameter smaller than the core radius can result in mass loss from this region. Additionally, in this prescription, there are two sets of conditions that, when met, result in a merger: $b\lt 2{R}_{\mathrm{core}}$ and $\sigma \lt {v}_{\mathrm{esc},\mathrm{core}}\sim 800$ km s⁻¹, or $2{R}_{{\rm{core}}}\lt b\lt 2\,{R}_{\odot }$ and σ < v_esc,surface ∼ 600 km s⁻¹. Note that this formula for drawing the impact parameter does not account for gravitational focusing, which can make head-on collisions more likely and lead to slightly higher fractional mass loss.

In order to expand this approach to stars of different masses, we assume that all stars have a similar structure to our Sun, allowing us to scale the quantities described above. In other words, we calculate the binding energy of the star as 1.61GM²/R, and assume that about two-thirds of the star's mass is concentrated in the inner third of the radius. While the true structure of a star depends on its mass, which determines the convective regions, stars with M ≳ 1 M_⊙ have a similar structure. Massive stars (M ≳ 10 M_⊙) have convective cores. However, their shorter main-sequence lifetime coupled with their low abundances in our sample means that very few of them collide. Future work may consider a more accurate approach to stars of other masses, perhaps by using different polytropic indices as was done by Rubin & Loeb (2011), but for the purposes of this study the approach detailed above provides a simple mass-loss recipe that complements the other prescriptions used herein (see Figure 2).

2.3.4. Fitting Formulae

In other simulations, we rely on fitting formulae from SPH studies of collisions at high velocities. Considering collisions in galactic nuclei, Lai et al. (1993) explores a wide variety of relative velocities, mass ratios, and impact parameters. This study provides a variety of fitting formulae, physically reasoned and based on their numerical results, that calculate the fractional mass loss and determine whether or not a collision results in a merger. Specifically, we use their Equations (4.4) and (4.8)–(4.10) and the fitting parameters in their Table 1. We label simulations that use these fitting formulae as "Lai+93."

As an alternative, we also use the fitting formulae in Appendix B of Rauch (1999). Based on a set of SPH simulations for high-impact velocities by M. Davies (1994, private communication), Equation (B1) relates the fractional mass loss of the two stars in question using their mass ratio, and Equation (B2) calculates the total fractional mass loss from the system in terms of the impact parameter, mass ratio, and relative velocity. The label "Rauch99" denotes simulations in this paper that use these equations.

In Figure 2, we plot the fractional mass loss expected from these fitting formulae for a head-on collision (b = 0) between two 1 M_⊙ stars as a function of distance from the SMBH. The equations from Rauch (1999) generally give a lower fractional mass loss than the other models used in this work. ⁶ We note that both Rauch (1999) and Lai et al. (1993) assume that stars are polytropes, which may provide an incomplete picture of the full range of collision outcomes (Freitag & Benz 2005). Additionally, Rauch (1999) note the difficulties in determining the post-collision radius of the star and approximate it as being equal to that of a main-sequence star with the same mass, an approach that is consistent with ours.

In the bottom panel of Figure 2, we also plot the maximum impact parameter, b_max, that will result in a merger as a function of distance from the SMBH, which acts as a proxy for the velocity dispersion. The red line uses the fitting formala provided in Lai et al. (1993), while the green dashed curve represents our simple toy-model approach, used as the default in all of our other simulations. Unlike Lai et al. (1993), whose equations give the specific impact parameter needed to achieve a merger at a given impact velocity, Rauch (1999) use an escape-velocity argument to determine whether or not a merger occurs, much like our own toy model. We therefore adopt the conditions described in Section 2.3.3 (green dashed line in the bottom panel of Figure 2) to determine whether or not a collision results in a merger in Rauch99 simulations. These conditions lead to more mergers overall than the Lai et al. (1993) fitting formula. However, as noted in Section 2.3.3, our formula for drawing the impact parameter does not account for gravitational focusing, which makes head-on collisions more likely. Our Lai+93 simulations may therefore underestimate the number of mergers because of the stricter b_max condition. Our Rauch99 simulations, which use the less restrictive b_max merger conditions, may provide a more accurate picture of mergers outside of 0.02 pc, where gravitational focusing becomes increasingly important.

2.4. Two-body Relaxation

As a star orbits the SMBH, it will experience frequent weak gravitational interactions with passing objects. Over time, the small effects of these interactions accumulate, altering the star's orbit. This process is known as relaxation and acts over an associated timescale, t_rlx, the amount of time needed for the star's orbital energy to change by order of itself (e.g., Binney & Tremaine 2008):

$$\begin{eqnarray}&&{t}_{\mathrm{rlx}}=0.34\displaystyle \frac{{\sigma }^{3}}{{G}^{2}\rho \langle {M}_{* }\rangle \mathrm{ln}{{\rm{\Lambda }}}_{\mathrm{rlx}}}.\end{eqnarray} \tag{ 7 }$$

In Figure 1, the blue lines show this timescale as a function of distance from the SMBH for a range of stellar density profiles.

We simulate the effects of relaxation in our code. Over each time step, we apply a small instantaneous velocity kick to the star, from which we calculate the new, slightly altered orbital parameters (see Lu & Naoz 2019 and Naoz et al. 2022 for the full set of equations). The velocity kick is drawn from a Gaussian distribution with a standard deviation ${\rm{\Delta }}{v}_{{\rm{rlx}}}/\sqrt{3}$, where ${\rm{\Delta }}{v}_{{\rm{rlx}}}={v}_{\mathrm{orbital}}\sqrt{{P}_{\mathrm{orbital}}/{t}_{{\rm{rlx}}}}$ (e.g., Bradnick et al. 2017; Naoz et al. 2022; Rose et al. 2022). Here, v_orbital and P_orbital are the orbital speed and period of the star in question. This prescription allows us to account for changes in the orbital parameters over time from interactions with the surrounding stars.

We note that we do not include resonant relaxation processes, which change the magnitude and orientation of an orbit's angular momentum vector (e.g., Rauch & Tremaine 1996; Hopman & Alexander 2006; Kocsis & Tremaine 2011). However, we do not consider any orbital inclinations as our models assume a spherically symmetric nuclear star cluster. Additionally, while scalar resonant relaxation may change the eccentricity of the orbit faster than standard two-body relaxation (e.g., Bar-Or & Fouvry 2018), we show in Appendix B that the orbital eccentricities have a weak effect on the final simulation results.

2.5. Stopping Conditions

At each time step in our code, we calculate the star's main-sequence lifetime, t_MS, according to its present mass: ${t}_{\mathrm{MS}}={10}^{10}{(M/{M}_{\odot })}^{-2.5}$. When a star's age exceeds its main-sequence lifetime, we terminate the simulation. In the case of a merger, we calculate the new age of the merged product using a "sticky-sphere model" (see, e.g., Kremer et al. 2020). Mixing during the collision and merger can introduce a fresh supply of hydrogen to the core, extending the lifetime of the merged star. The dimensionless parameter f_rej, ranging from 0.1 to 1, quantifies the degree of mixing and therefore the degree of rejuvenation (e.g., Hurley et al. 2002; Breivik et al. 2020). While the true value of f_rej is challenging to ascertain since it can depend on the details of each collision and stellar structure, we adopt f_rej = 1, the most conservative value. The new age of the star, t_age,f , can be calculated as follows:

$$\begin{eqnarray}&&{t}_{\mathrm{age},f}={f}_{\mathrm{rej}}\displaystyle \frac{{t}_{\mathrm{MS},M+{M}_{\odot }}}{M+{M}_{\odot }}\left(\displaystyle \frac{M\times {t}_{\mathrm{age},i}}{{t}_{\mathrm{MS},M}}+\displaystyle \frac{{M}_{\odot }\times {t}_{\mathrm{sim}}}{{t}_{\mathrm{MS},\odot }}\right).\end{eqnarray} \tag{ 8 }$$

Above, t_age,i is the age of the star in question before it collides with a solar-mass star—recall that we assume that solar-mass stars comprise the surrounding cluster—and t_MS,M is its main-sequence lifetime. The star with which it collides has a main-sequence lifetime t_MS,⊙ and age equal to the time elapsed in the simulation, t_sim. ${t}_{\mathrm{MS},{\rm{M}}+{{\rm{M}}}_{\odot }-{{\rm{M}}}_{\mathrm{lost}}}$ denotes the main-sequence lifetime of the merger product. Our simulations consider only one population of stars, whose age is set by the time elapsed in the simulation. However, in the actual nuclear star cluster, star formation episodes may replenish the stellar population. Observations suggest two major star formation episodes in the GC, separated by billions of years (e.g., Chen et al. 2023). If a star from the older, evolving population collides with a younger star, as opposed to one from its own formation episode, the merger product will be even more rejuvenated. We may therefore be underestimating the degree of rejuvenation and new stellar lifetime in a small fraction of the collisions. These longer-lived stars may go on to experience more collisions.

In addition to the stellar age, we include two other stopping conditions in our code. Relaxation processes may drive a star's orbital eccentricity to high values, shrinking its pericenter passage to fall within the tidal radius of the SMBH. We include a stopping condition to account for these tidally disrupted stars:

$$\begin{eqnarray}&&a{(1-e)\lt R(2{M}_{\bullet }/M)}^{1/3}.\end{eqnarray} \tag{ 9 }$$

Lastly, some stars undergo significant mass loss through one or more collisions. The code terminates once the mass of a star falls below 0.1 M_⊙.

We begin by comparing the different mass-loss prescriptions detailed in Sections 2.3.2, 2.3.3, and 2.3.4. We therefore run several simulations with a sample population of only 1 M_⊙ stars and identical initial conditions, constraints we later relax in Section 5. At this stage we do not include relaxation in order to build a clear picture of the differences between our mass-loss prescriptions. Figure 3 juxtaposes the results. The initial conditions, 1000 solar-mass stars distributed uniformly in log distance from the SMBH, are shown in gray, while red dots show the final masses versus distance from the SMBH. The blue dots show the total change in each star's mass over the simulation. Each row corresponds to a different mass-loss prescription, indicated by the label in the left corner.

Zoom In Zoom Out Reset image size

All of the mass-loss prescriptions yield qualitatively similar results, consistent with the expectations from Figure 1. Close to the SMBH, collisions work to destroy the stars, while ∼1 pc from the SMBH few collisions occur. In between these two extremes, one or more collisions occur, and they result in mergers. Generally, our results suggest that multiple collisions between 1 M_⊙ stars can produce more massive stars. These massive products are localized between 0.01 and 0.1 pc from the SMBH. In this region, a few conditions must align to support the formation of massive stars. First, the collision timescale is shorter, in some cases by an order of magnitude, than the main-sequence lifetime and several collisions can take place. Second, the velocity dispersion is small enough that stars can capture each other and the factional mass loss is low. The final mass of the stars formed through collisions therefore depends on the specific mass-loss prescription adopted.

For these initial conditions, all of our mass-loss prescriptions can produce a >8 M_⊙ star. However, the Rauch99 prescription was the only one to do so consistently. This outcome is not surprising given that the Rauch99 equations result in less mass loss compared to the others (see Figure 2), producing more massive stars. Over several Rauch99 simulations with 1000 solar-mass stars each, we find that a few stars out of the 1000 will exceed 10 M_⊙ through mergers with other low-mass stars. We scale this number using a realistic density profile normalized by the M–σ relation (Bahcall & Wolf 1976; Tremaine et al. 2002). Accounting for the number of stars residing between 0.01 and 0.1 pc in the GC, we estimate that there may be more than ∼100 stars over 10 M_⊙ formed through collisions between lower-mass stars. This estimate assumes that collisions tend to produce less mass loss, as in the Rauch99 prescription, and result in only moderate rejuvenation (f_rej = 1 in Equation (8)). If mergers in the GC are more effective at rejuvenating stars, then collisional processes should lead to even more massive stars and in greater abundance.

Another difference between mass-loss prescriptions is their ability to produce low-mass remnants. Often, a collision within ∼0.01 pc of the SMBH, where the velocity dispersion is high, will only remove the outer layers of a star. This outcome is common because most collisions have a nonzero impact parameter (see also discussion in Rauch 1999). The contrast between the upper row and last two rows in Figure 3 illustrates the importance of the impact parameter. The upper row does not include an impact parameter dependence, and as a result it overestimates the number of stars that are destroyed. Rauch99, on the other hand, predicts many low-mass collision remnants, a result discussed in their original study (Rauch 1999). The evolution of these remnants post-collision is uncertain. However, from a dynamical standpoint, these objects may be long-lived. Their small cross section leads to a longer collision timescale, making it less likely that they will undergo another destructive encounter.

In the previous section, we demonstrate that the general collision outcome trends described in Figure 1 hold. Relaxation has the potential to complicate the physical picture, however, by allowing stars to diffuse into different regions from the SMBH. We show the effects of this dynamical process in Figure 4. The simulations shown in the Figure use identical initial conditions as those in the last two rows of Figure 3, except with the addition of relaxation. One main difference is that, due to two-body relaxation, low-mass remnants can migrate out of the region in which they form, around 0.01 pc, to larger distances from the SMBH. We therefore predict a diffuse population of peculiar low-mass stars throughout the GC. In contrast, stars that experience the most mass growth are still localized to 0.01–0.1 pc from the SMBH. These stars must form in this region because otherwise either the collision timescale is too long to achieve multiple mergers or the collisions are destructive in nature. Additionally, as a star grows in mass, its main-sequence lifetime becomes short compared to the relaxation timescale (Figure 1). As a result, the most massive stars have insufficient time to diffuse out of the region in which they form. We therefore expect massive stars formed through stellar collisions to be limited within ≲0.1 pc from the SMBH.

Zoom In Zoom Out Reset image size

Thus far, we have considered only uniform populations of 1 M_⊙ stars. In this section, we use a Kroupa IMF to generate our initial conditions and explore general trends in the stellar mass distributions. We consider two mass-loss prescriptions, Lai+93 and Rauch99, because they encapsulate the range of results in Figure 3. Recall that Rauch99 favors mergers with minimal mass loss, while the Lai+93 equations lead to greater mass loss and have a stricter merger condition.

Figure 5 shows the results of two simulations, Rauch99 on the left and Lai+93 on the right, with the aforementioned initial conditions. In the bottom row, we plot the final masses and semimajor axes of the stars in red and the initial masses and semimajor axes in gray. The bottom row also shows the change in mass of each star versus its initial position. For effective visualization, we divide the nuclear star cluster into three regions that roughly correspond to different collision outcomes based on Figure 1. The regimes are as follows:

• 1.  
  r < 0.01 pc (top row, light blue border), where the velocity dispersion is high and we expect repeated collisions lead to mostly destruction of stars.
• 2.  
  0.01 < r < 0.1 pc (middle row, dark blue border), where the velocity dispersion and relatively short collision timescale allows mergers to take place.
• 3.  
  0.1 < r < 2 pc (third row, red border), where in this regime, while the velocity dispersion allows mergers to take place, the collision timescale is longer than the stellar lifetimes. The number of collisions is therefore low.

These three rows in Figure 5 depict the the number of stars within a certain mass range as a function of time. The solid lines represent the expectation in the absence of collisions. The solid lines decrease when stars within each mass range reach the end of their main-sequence lifetimes. The symbols show stellar abundances when collisions are included in the simulation. In the bottom row of the figure, we also include plots in the style of Figure 4 with the different regimes highlighted to guide the eye.

Zoom In Zoom Out Reset image size

In the first region (top row), the populations of all the low-mass stars decline over time. Collisions destroy these stars before the can evolve off the main sequence. However, multiple collisions are required in order to destroy the stars. As can be observed in the figure, it takes about a billion years to halve the population of low-mass stars in both the Rauch99 and Lai+93 simulations. This half-life is about an order of magnitude longer than the collision timescale. We attribute this longevity to the fact that most collisions are not head-on and therefore less destructive. Additionally, stars losing mass through collisions can temporarily increase the number of stars in a lower mass range. For example, as 1–1.5 M_⊙ stars lose mass, they will temporarily increase stars in the 0.5–1 M_⊙ mass range.

Between 0.01 and 0.1 pc, the number of stars in the ranges M < 1 M_⊙ and 1–1.5 M_⊙ still decrease. However, while some of these stars may undergo destructive collisions, it is apparent from the final masses, as shown in the last row of the figure, that their changes in number are due mostly to mergers. This can be seen most clearly in the Rauch99 simulation, where there is a visible increase in the number of 1.5–3 M_⊙ stars such that they outnumber the 1–1.5 M_⊙ stars. Additionally, increasing the mass of the stars even slightly causes the stars to evolve off the main sequence more quickly. As a result, the red symbols, for example, decrease more quickly compared to the expectation based on the original population (solid line). The Lai+93 prescription leads to more moderate mass growth than the Rauch99 one. As a result, the faster decline of the 1–1.5 M_⊙ stars than expected is likely due to accelerated evolution, not stars leaving this mass range entirely. Lastly, in the 0.1 < r < 2 pc, similar effects from mergers are present, but to a far lesser degree.

We suggest that collisions may result in more massive, young-seeming stars. In particular, because collisions need time to act, a massive stellar population may emerge after about a billion years. This emergent population skews toward higher masses, yielding a mass function that resembles a top-heavy one. This behavior is depicted in Figure 6, which shows stellar mass distribution for a few representative snapshots in time. In this figure, we again show our two main simulations, also depicted in Figures 4 and 5, in the left and right columns. The simulation shown in the middle column is identical to the Rauch99 simulation on the left, except the initial distances from the SMBH of the 1000 stars are drawn from a Bahcall–Wolf distribution (α = 1.75). This initial condition does not adequately sample stars close to the SMBH. However, this simulation provides a check on mass functions from the other two simulations, which overrepresent stars close to the SMBH in population demographics because most of the stars in the nuclear star cluster are located further from the SMBH. Each subplot shows the mass function expected at a specific age for the population based on a Kroupa IMF in gray. The red dashed line shows the mass function at a given time for the entire cluster. The blue dashed–dotted and green dashed lines show the mass distributions for the stars internal to 0.1–0.01 pc, respectively. We do not show the latter for the middle simulation because stars within this region are not adequately sampled for any meaningful result. While Figure 4 shows snapshots at only four different times, the expanded Figures 12 and 13 in the Appendix display the full evolution of the mass function.

Zoom In Zoom Out Reset image size

As shown in Figure 4, the Rauch99 model, which favors mergers with minimal mass loss, leads to an excess of more massive stars in the inner 0.1 pc (see bottom two rows of the left and middle columns). The Lai+93 column exhibits a similar effect, though to a lesser degree. These stars may appear rejuvenated and younger, especially compared to the surrounding cluster. Interestingly, the IMF of the young stars within ∼0.5 pc of the SMBH has been suggested to be top heavy, with an excess of massive stars (e.g., Paumard et al. 2006; Bartko et al. 2010; Lu et al. 2009, 2013). Furthermore, the origins of young massive stars in the nuclear star cluster remain an open question; it remains uncertain whether they could have formed in situ in the extreme gravitational potential of the SMBH (e.g., Jackson et al. 1993; Morris 1993; Gerhard 2001; Ghez et al. 2003; Genzel et al. 2003; Levin & Beloborodov 2003; Christopher et al. 2005; Davies & King 2005; Paumard et al. 2006; Montero-Castaño et al. 2009; see Section 6 of Genzel et al. 2010 and Section 7 of Alexander 2005 for a review). It has been suggested that stellar collisions and mergers may instead create these massive stars, explaining the so-called "paradox of youth" of the S-star cluster (e.g., Genzel et al. 2003; Ghez et al. 2003; Antonini et al. 2010; Stephan et al. 2016, 2019). Our findings are consistent with a collisional formation scenario.

Given the abundance of collisions within 0.1 pc, here we consider what the merger products might look like post-collision, before the stellar object, a newly merged star, has had time to recollapse over a Kelven–Helmholtz timescale. In particular, these merged stars may have connections to gas- and dust-enshrouded stellar objects, or G objects, observed in the GC (see Ciurlo et al. 2020). Several works in the literature have addressed the mysterious origin of the G-object population (e.g., Burkert et al. 2012; Murray-Clay & Loeb 2012; Schartmann et al. 2012; Madigan et al. 2017; Zajaček et al. 2017; Owen & Lin 2023). For example, stellar mergers from binary systems may explain the presence of these distended stellar objects (e.g., Witzel et al. 2014; Prodan et al. 2015; Stephan et al. 2016, 2019). In the nuclear star cluster, gravitational perturbations from the SMBH on a stellar binary's orbit can drive the binary to merge (e.g., Antonini & Perets 2012; Prodan et al. 2015; Stephan et al. 2016, 2019). Collisions between other species, like compact objects and RGs, have also been invoked to explain the formation of G objects (Mastrobuono-Battisti et al. 2021).

We estimate the number of G objects from main-sequence stellar collisions in the GC in Figure 7. We begin by calculating the number of collisions that occur per year at each distance from the SMBH. A post-collision merger product remains puffy and extended for a thermal timescale (see Stephan et al. 2016, 2019). ⁷ Therefore, the number of G-type objects visible at a given distance from the SMBH is equal to the collision rate times 10⁷ yr, plotted in the second row of the Figure. The last row shows the cumulative number of G objects as a function of distance from the SMBH for three different stellar density profiles, α = 1.25, 1.5, and 1.75. This analytic calculation assumes a uniform population of solar-mass stars. We do not include collisions internal to about 0.008 pc, based on Figures 4 and 5, because collisions in this region are largely destructive and do not result in mergers. For the least cuspy of the profiles considered, α = 1.25, we expect approximately 10 G objects, a number consistent with the current G-object population (Ciurlo et al. 2020). For the other two steeper density profiles, α = 1.5 and 1.75, G objects should number in the tens and hundreds, respectively. The density profile in the GC in the present day appears to be shallower than a Bahcall–Wolf profile, lying between about α = 1.25 and 1.5 (e.g., Genzel et al. 2003; Gallego-Cano et al. 2018).

Zoom In Zoom Out Reset image size

Observations of the nuclear star cluster indicate that the surface number density of RGs plateaus within about 0.3 pc of the SMBH (e.g., Genzel et al. 1996; Buchholz et al. 2009; Do et al. 2009; Baumgardt et al. 2018; Gallego-Cano et al. 2018; Habibi et al. 2019). This plateau diverges from the cusp-like density profile expected for an older, evolved stellar population (e.g., Bahcall & Wolf 1976; Linial & Sari 2022; see the latter for the distribution expected from steady-state scatterings). The observed missing RG cusp has motivated several proposed mechanisms to destroy RGs near the SMBH, many of which have been unable to completely explain the phenomenon (Davies et al. 1998; Alexander 1999; Bailey & Davies 1999; Dale et al. 2009; Zajaček et al. 2020). Of the proposed mechanisms, the most promising relies on interactions between the RGs and a gaseous disk in the GC, which strip the RGs of their outer layers (Amaro-Seoane & Chen 2014; Amaro-Seoane et al. 2020), though we also note that the missing cusp may have a dynamical explanation, without requiring the destruction of the individual RGs (e.g., Merritt 2010; Antonini 2014). Here, we explore the possibility that main-sequence stellar collisions contribute to the dearth of RGs. As shown in Figure 1, 1–3 M_⊙ stars within ∼0.1 pc of the SMBH are likely to experience a collision before evolving off the main sequence. To ensure a deficit in RGs, these stars must lose significant mass—or be completely destroyed—before they can evolve off the main sequence, or mergers may accelerate their evolution off the main sequence.

Recently, Mastrobuono-Battisti et al. (2021) undertook a comprehensive N-body study of collisions between many different species in the GC and comment on a variety of applications, including the missing RGs. Their results suggest that collisions between main-sequence stars cannot explain the underabundance of RGs in the GC. They assume that the stars merge with no mass loss when the velocity dispersion is between 100 and 300 km s⁻¹, as is the case in their model. However, a velocity dispersion of 300 km s⁻¹ corresponds to a distance of about 0.07 pc from the SMBH. Furthermore, a 1 M_⊙ star may diffuse into regions of larger velocity dispersion over its main-sequence lifetime, which is longer than the ∼1 Gyr relaxation timescale. In this section, motivated by Figures 5 and 6, we assess whether collisions closer to the SMBH can affect the RG count in the GC.

We generate RG surface-density profiles at different snapshots in time from our fiducial Rauch99 simulation, shown for example in the left column of Figure 5. For adequate resolution, we increase the number of stars in our sample from 1000 to 4000. Only a subset of this population will be RGs at any given time. The RGs used to generate the surface-density profiles shown in Figure 8 number at about 100 each. We bin the number of RGs by distance from the SMBH. We then take the number of RGs formed per star and scale it by the number of stars expected at that distance, normalized using the M–σ relation (Tremaine et al. 2002). Figure 8 shows the results at three different ages for the stellar population, 3, 4, and 8 Gyr in gray, peach, and red, respectively. We compare these results to data from Figure 9 of Gallego-Cano et al. (2018), represented by the blue star-shaped dots. Extinction may explain the dip in the RG count data at 0.2 pc (Buchholz et al. 2009; Gallego-Cano et al. 2018). We have scaled our simulated results to coincide roughly with the number of observed RGs at 1 pc. In each panel, we also include a Bahcall–Wolf profile to guide the eye.

Zoom In Zoom Out Reset image size

Both our simulated results and the observed profile lack a cusp internal to about 0.1 pc. In this region, mergers occur. As a result, these stars evolve off the main sequence more rapidly. This process can boost the number of RGs at earlier times. The timing at which these stars evolve off the main sequence has a radial dependence, based on the collision timescale. A population over a few billion years old may have already experienced this accelerated boost of RGs, leading to a deficit in present time. Additionally, in this region, where multiple collisions are possible, many stars may undergo enough mergers to exceed 3 M_⊙ (see also the last panel of Figure 4). Lastly, more massive stars spend less time as RGs, so merging 1 M_⊙ stars into 2 M_⊙ stars, for example, reduces the chances of catching and observing them as RGs.

The discussion above pertains to stars outside ∼0.01 pc, the innermost distance from the SMBH included in Figure 8. However, as seen in Figure 5, stellar collisions destroy most low-mass (≲3 M_⊙) stars within ∼0.01 pc of the SMBH before they can evolve off the main sequence. Thus, we also expect a clear RG deficit in this region. We note that it remains uncertain whether stellar collisions, even at extremely high velocities, can destroy stellar cores. The first one or few collisions in this ≲0.01 pc region can easily remove that outer layers of a star. In our simulations, later collisions can go on to destroy the low-mass collision remnants. However, if these stellar cores prove more resilient against collisions, our results may need to be revisited.

We stress that the simulated curves in Figure 8 should not be used as predictions for the number of RGs at a given distance from the SMBH. This exercise is meant to draw attention to the general shape of the RG surface-density profile based on our results. The actual RG density is likely a composite of several of these theoretical curves. The GC contains populations of different ages and its star formation history remains an open area of research (e.g., Pfuhl et al. 2011; Schödel et al. 2020; Nogueras-Lara et al. 2021; Chen et al. 2023). Additionally, our simulations consider an evolving sample of stars embedded in a fixed, unchanging cluster. In an actual stellar cluster, a merger between two stars results in a single object. Each collision should therefore reduce the number of particles in the system. In our simulations, however, the number of particles, or stars, remains the constant in our sample. Our final numbers may therefore overcount by a factor of 2. Lastly, in addition to the caveats mentioned above, we also assume an initial Bahcall–Wolf profile (α = 1.75) for the stars. Two-body scattering in the presence of a BH population in the nuclear star cluster may instead cause the stars to settle on a slightly shallower (α = 1.5) profile within a critical radius from the SMBH (Linial & Sari 2022). A shallower initial profile leads to a longer collision timescale, lessening the impact of this dynamical process on the density profile within 0.1 pc (S. C. Rose & M. MacLeod 2023, in preparation).

We also note that the number of RGs depends on the interplay between distance from the SMBH and time. The profile's shape is determined by when mergers occur at a given distance from the SMBH and how many mergers are possible before the star evolves off the main sequence; it is sensitive to both the destruction and formation of RG progenitors via collisions. For example, in Figure 5 mergers produce an abundance of RG progenitors at about 1 Gyr. Figure 9 presents two examples of RG profiles that lack an apparent plateau within 0.1 pc. In conjunction, Figures 8 and 9 show that collisions can produce both cuspy and core-like RG density profiles depending on the age of the population. Interestingly, the density profiles in the latter figure still lack RGs within ∼0.03 pc. We speculate that this dearth occurs for two reasons: collisions begin to become destructive for the lowest-mass stars in this region, while stars that are not destroyed instead experience multiple mergers. In other words, RG progenitors are not being replenished at the rate they are being eliminated.

Zoom In Zoom Out Reset image size

Gallego-Cano et al. (2018) and Schödel et al. (2018) find that only the brightest giants exhibit a flattening of their profile within ∼0.1 pc, while observations of the less bright main-sequence stars and subgiants are consistent with a cusp. Schödel et al. (2018) note that the type of stars being observed and the degree of contamination from young, unrelaxed stellar populations are uncertain, especially as they rely on diffuse light to probe the stellar profile. However, their results imply that the mechanism destroying the RGs must act after the stars have ascended the giant branch. The mechanism we discuss here, main-sequence stellar collisions, would deplete the RG progenitors before they evolve off the main sequence. Our results do not necessarily contradict the observations for the reasons discussed above: the RG density profile depends on the age of the parent population(s) and the initial density profile assumed for the stars. We reserve a more thorough examination of the surface-density profile for future work. However, our results suggest that main-sequence stellar collisions merit further consideration in the context of the evolved population.

Collisions that take place in the dense nuclear star cluster surrounding our Galaxy's central SMBH can shape the stellar population. Within about 0.1 pc, every star undergoes at least one collision during its lifetime. More generally, within the inner 1 pc, we estimate that 10% of the stellar population will experience at least one collision while on the main sequence. In this study, we characterize the outcomes of these collisions as a function of distance from the SMBH. In particular, we consider three regimes: (1) within ∼0.01 pc of the SMBH, where the velocity dispersion is large and collisions are destructive; (2) 0.01–0.1 pc, where collisions are common and result in stellar mergers; and (3) outside of ∼0.1 pc, where collisions do not occur very frequently but always result in mergers. The details of each collision, specifically the amount of mass ejected during the collision, are more difficult to determine. In particular, the mass lost from the stars depends on their relative velocity, mass ratio, structure, and impact parameter (e.g., Lai et al. 1993; Freitag & Benz 2002). We adopt different prescriptions to treat the mass loss, including a few based on fitting formulae from SPH studies (Lai et al. 1993; Rauch 1999). Our models also account for two-body relaxation, which works to change the semimajor axes and eccentricities of the stellar orbits in the GC over time. We find the following:

• (i)  
  The formation of massive stars. The mass-loss prescriptions all lead to qualitatively similar results, consistent with our expectations for the different regimes. In particular, in the 0.01–0.1 pc regime, stars with mass greater than ∼8 M_⊙ can form from mergers between 1 M_⊙ stars. Furthermore, in the case that collisions result in only moderate mass loss even at high impact velocities, as is the case for the Rauch99 equations, repeated collisions between low-mass stars can consistently form massive (≳10 M_⊙) stars. These massive products are mainly localized within the ∼0.01–0.1 pc region, as there is insufficient time for them to diffuse to distances ∼1 pc from the SMBH. However, as can be seen in Figures 4 and 5, significant mass growth through mergers remains possible, albeit less probable, at ≲1 pc from the SMBH. We predict that there should be at least 100 massive stars formed from lower-mass stars in this region. Generally, mergers may cause the mass function of the population to appear top heavy (see Figure 6, in keeping with observations of this region; e.g., Paumard et al. 2006; Bartko et al. 2010; Lu et al. 2009, 2013).
• (ii)  
  A diffuse population of low-mass remnants. Collisions can divest many main-sequence stars of their outer layers, creating a low-mass stellar object. This result is consistent with previous studies (Rauch 1999). These collision remnants can migrate through relaxation processes from their place of formation, ≲0.01 pc, to greater distances from the SMBH. We therefore predict a population of low-mass peculiar stars throughout the GC. Additionally, the mass ejected through grazing or off-axis collisions may produce gas and dust features in the GC, like the observed X7 (e.g., Ciurlo et al. 2023), whose distance from the SMBH is comparable to where σ becomes ≲v_esc (see Figure 1).
• (iii)  
  The survival of low-mass stars within 0.01 pc. While collisions within 0.01 pc tend to be very destructive for the stars, our results suggest that their populations may persist for longer than expected based on the collision timescale alone, similar to findings in other studies (e.g., Rauch 1999). In fact, it takes ∼1 Gyr to halve the populations of low-mass (≲3 M_⊙) stars. This timescale is the same order of magnitude as the relaxation timescale in this region (Figure 1). Episodes of star formation, not included in this study, may therefore help replenish the stellar population in this region. Depending on the age of the population, there may be an undetected population of low-mass stars and collision remnants in the innermost region of the nuclear star cluster.
• (iv)  
  G objects. Collisions between main-sequence stars may have connections to G objects, stellar objects enshrouded in gas and dust observed in the GC (e.g., Ciurlo et al. 2020). A merger between two collided stars may produce a puffy G-like object that recollapses over the thermal timescale of a solar-mass star, 10⁷ yr (similarly, for binary mergers, see e.g., Stephan et al. 2019). We provide predictions for the number of G-type objects expected in the GC from stellar collisions given different stellar density profiles. For a stellar density profile with index α between 1.25 and 1.5, there may be ∼10 G objects observable in the inner parsec of the GC.
• (v)  
  Connections to the missing red giants. Observations of the GC indicate a deficit of RGs within about ∼0.3 pc of the SMBH (e.g., Genzel et al. 1996; Buchholz et al. 2009; Do et al. 2009; Gallego-Cano et al. 2018; Habibi et al. 2019). We consider whether main-sequence stellar collisions may help explain this observational puzzle. We find that within ∼0.01 pc of the SMBH, stellar collisions destroy most low-mass stars before they can evolve off the main sequence (see Figure 5). Thus, we expect a lack of RGs in this region. Further from the SMBH, between ∼0.01 and 0.1 pc, mergers can result in fewer RGs. As seen in Figure 8, the RG surface-density profile plateaus in this region, an effect similar to what is observed in the GC. Both the profiles from our models and observations diverge from the density cusp expected for an old, dynamically relaxed population: α = 1.75 in Equation (1).

Collisions play a crucial role in shaping the stellar populations in the inner parsec of the GC. We show that this dynamical channel leads to a number of interesting observables, including rejuvenated, massive stars near the SMBH, G objects, and missing evolved stars. Our simulations focus on the dynamics of this dense environment using a statistical approach. We leverage previous SPH studies (e.g., Lai et al. 1993) and intuitive toy models. However, our results highlight the need for future work to explore the rich interplay between hydrodynamics, stellar evolution, and collisions in these dense, dynamic environments.

We thank Anna Ciurlo, Cheyanne Shariat, Jessica Lu, Mark Morris, Morgan MacLeod, Ylva Gotberg, and Enrico Ramirez-Ruiz for helpful discussion. S.R. thanks the Charles E. Young Fellowship, the Dissertation Year Fellowship, and the Thacher Fellowship for support. S.R., S.N., and R.S. acknowledge the partial support of the NSF-AST 2206428 grant. S.R. and S.N. acknowledge the partial support from NASA ATP grant No. 80NSSC20K0505. S.N. thanks Howard and Astrid Preston for their generous support.

We base our toy "core+envelope" model on the structure of the Sun as described in Christensen-Dalsgaard et al. (1996). ⁸ The lower panel of Figure 10 shows the mass enclosed within a given radius. Based on the mass enclosed within a given radius r, we calculate the escape velocity from a sphere with equivalent mass and radius, ${v}_{\mathrm{esc}}=\sqrt{2{{GM}}_{\mathrm{enc}}/r}$. Note that this quantity is not the true escape velocity from some point r within the Sun as that would require integrating from r to infinity to account for the outer layers. We plot this quantity as a function of distance from the center of the Sun all the way to the surface in the upper panel of Figure 10. Because most of the Sun's mass is concentrated at its center, this quantity peaks at 0.33 R_⊙, a radius which encloses 66% of the Sun's mass. We define the "core" of the star such that its radius coincides with this point. The "envelope" comprises the material external to this point, which is more loosely bound. Assuming that all stars have a similar structure, we scale this model to all stars within our considered mass range; we take 66% of the star's mass to be within 0.33 × R_star of the star's center and define this region as the "core." We then calculate collision outcomes as detailed in Section 2.3.3. This approach gives a straightforward, intuitive way of calculating the fractional mass loss and determining the outcome of a collision.

Zoom In Zoom Out Reset image size

In this appendix, we relax the assumption that the orbital eccentricities of the stars are on a thermal distribution, in which the probability density f(e) is proportional to e. We test two other initial eccentricity distributions. First, we test a uniform initial distribution, in which all eccentricities are equally likely. The average eccentricity for this population is 0.5, as opposed to ∼0.67 for a thermal distribution. We also test a distribution with the form f(e) ∝ e^0.5, which favors low eccentricities. The results of these simulations are shown in Figure 11. For ease of comparison to Figure 3, they do not include relaxation. Other than the eccentricity distribution, the initial conditions are identical to the simulation shown in the last row of Figure 3.

Zoom In Zoom Out Reset image size

The first and second rows of Figure 11 show simulations that use f(e) ∝ e^0.5 and a uniform eccentricity distribution, respectively. The results are qualitatively similar to Figure 3. While the less eccentric initial conditions result in fewer massive stars, their formation is still likely. This result is expected based on Rose et al. (2020), who find that the collision timescale of a very eccentric orbit (e > 0.95) and a circular orbit differ by a factor of 2.

In this section, we show expanded versions of Figure 6, including more time steps to detail the evolution. Figure 12 shows stellar mass distributions before 1 Gyr. Figure 13 shows snapshots from 1 Gyr onward. Together, they provide a more detailed evolution of the mass function. Mergers skew the mass function toward higher masses, giving it a slightly top-heavy appearance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our simulations assume a background population of stars with some density profile that remains fixed. However, over time we expect stellar collisions to modify the distribution, in particular by eroding the cusp near the SMBH. In the absence of star formation, the new, shallower density profile can only lead to a lower collision rate. As a result, our predictions of the rate of destruction of stars may be overpredicting their decline in population. Additionally, two-body scattering with a population of heavier BHs may cause the stars to lie on a profile with index α = 1.25 (Linial & Sari 2022), and observations of the GC's nuclear star cluster suggest that the profile may be shallower (e.g., Genzel et al. 2003; Gallego-Cano et al. 2018). Therefore, we are including a version of the Rauch99 simulation with a α = 1.25 density profile for the surrounding star cluster, as opposed to our default α = 1.75. The most realistic picture probably lies between the α = 1.75 and α = 1.25 cases.

We show results in the same style as Figures 3 and 5. In Figure 14, we plot the results of a Rauch99, α = 1.25 that does not include the effects of relaxation. There are fewer mergers due to the shorter collision timescale. The merger products are therefore less massive compared to the α = 1.75 case. Additionally, the stars are not destroyed as efficiently. As can be seen in the upper panel of Figure 14, the low-mass collision remnants persist for longer closer to the SMBH.

Zoom In Zoom Out Reset image size

In a similar fashion, Figure 15 is identical to the left column of Figure 5 except for the assumed density profiles. More massive stars are still formed through collisions, but they are generally limited to less than 3 M_⊙ because fewer collisions occur. For the same reason, the decline of the low-mass stars occurs more gradually than in the α = 1.75 case.

Zoom In Zoom Out Reset image size