Collisional Shaping of Nuclear Star Cluster Density Profiles

The properties of the surrounding fixed cluster govern the dynamical processes that shape our evolving sample of stars. In this section, we describe the properties of our model nuclear star cluster. The stellar mass density of the cluster is described by a power law:

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

where α is the slope and r_• denotes the distance from the SMBH. The density profile within the sphere of influence of the SMBH is normalized using ρ₀ = 1.35 × 10⁶ M_⊙ pc^‒3 at r₀ = 0.25 pc, based on observations of this region (Genzel et al. 2010). In our simulations, we test three values of α, 1.25, 1.5, and 1.75, to encapsulate the range of theoretical predictions and observed density profiles.

The velocity dispersion within the cluster is also a function of distance from the SMBH:

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

where α is the slope of the density profile and M_• is the mass of the SMBH (Alexander 1999; Alexander & Pfuhl 2014). We set M_• equal to 4 × 10⁶ M_⊙, the mass of the Milky Way's SMBH (e.g., Ghez et al. 2003). For a uniform mass cluster of 1 M_⊙ stars, the number density n is simply $\tfrac{\rho ({r}_{\bullet })}{1\,{M}_{\odot }}$. Together, the density and velocity dispersion in the nuclear star cluster determine the frequency of interactions between stars and set key timescales for various dynamical processes.

2.2.1. Collision Rate

In dense star clusters, direct collisions between objects become possible (e.g., Dale & Davies 2006; Dale et al. 2009; Rubin & Loeb 2011; Mastrobuono-Battisti et al. 2021; Rose et al. 2022, 2023). The timescale for a direct collision can be estimated as ${t}_{\mathrm{coll}}^{-1}=n\sigma A$, where A is the cross section of interaction, n is the number density, and σ is the velocity dispersion, here a measure of the relative velocity between objects. For a direct collision, the cross section of interaction A is the physical cross section enhanced by gravitational focusing. The eccentricity of a star's orbit about the SMBH, e_•, can also affect the collision timescale; more eccentric orbits have shorter collision timescales compared to circular orbits with the same semimajor axes (e.g., Rose et al. 2020). Including the eccentricity dependence, the collision timescale can be written as

$$\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{ 3 }$$

where f₁(e_•) and f₂(e_•) are given by Rose et al. (2020; Equations (20) and (21)), G is the gravitational constant, a_• is the semimajor axis of the star's orbit, and r_c is the sum of the radii of the colliding stars. The second term in the parentheses comes from gravitational focusing (e.g., Binney & Tremaine 2008). We plot the collision timescale in red in Figure 1 for the range of density profiles considered in this study. α = 1.75 is shown in the solid line, while α = 1.25 is shown in the dashed line. This timescale assumes a uniform population of solar mass stars. Therefore, r_c = 2R_⊙.

Zoom In Zoom Out Reset image size

2.2.2. Two-body Relaxation

Even more frequent than direct physical collisions are weak gravitational interactions between passing stars. These interactions cause the orbital parameters to change slowly over time in a diffusion process. Over the so-called relaxation timescale, the star's orbital energy and angular momentum change by order of themselves. Like the collision timescale, the relaxation timescale t_rlx depends on properties of the cluster such as the density and velocity dispersion, both functions of r_•. The timescale can be expressed as

$$\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{ 4 }$$

where 〈M_*〉 is the average mass of the objects in the cluster, here taken to be 1 M_⊙, and $\mathrm{ln}{{\rm{\Lambda }}}_{\mathrm{rlx}}$ is the coulomb logarithm (e.g., Binney & Tremaine 2008; Merritt 2013). We plot this timescale in blue in Figure 1 for a range of density profiles, spanning α = 1.25 (dashed)–1.75 (solid). The timescale is less than or comparable to the duration of our simulations, 10 Gyr, shown in gray in Figure 1. Additionally, while outside of 0.1 pc, the collision timescale is long compared to the total simulation time, and relaxation processes can change the orbits of the stars and move them into regimes where collisions become likely. The reverse is also true: stars that begin in regions where collisions are common can migrate further from the SMBH. It is therefore important to account for this diffusion process in our semianalytic models.

We use a toy model developed by Rose et al. (2022, 2023) to simulate the effects of collisions on the star cluster. This model follows a subset of 1000 stars of varying masses, drawn from a Kroupa initial mass function (IMF) with an assumed range of 0.5–100 M_⊙, embedded in a fixed cluster of 1 M_⊙ stars. The index α of the surrounding star cluster's density profile is treated as a free parameter with a default value of 1.75, the expectation for an old, dynamically relaxed population (e.g., Bahcall & Wolf 1976). While the value of α sets the spatial distribution of stars in the nuclear star cluster, and therefore the probability and frequency of collisions, we draw the semimajor axes of our sample stars from a uniform distribution in log distance from the SMBH. We sample stars using this distribution because it allows us to understand the effects of collisions at all distances from the SMBH. However, calculations related to each star's dynamical evolution still rely on the spatial distribution of the surrounding stellar cluster. The orbital eccentricities of the stars in our sample are drawn from a thermal distribution.

The code takes a statistical approach to collisions. It computes the probability of a collision occurring over a timestep Δt, taken to be 10⁶ yr, as Δt/t_coll. We then draw a random number between 0 and 1. If the number is less than or equal to the collision probability, we treat the star as having collided and update its properties given the models described in Section 2.4. This prescription repeats until the code has reached the desired simulation time, 10 Gyr, or the star has reached the end of its main-sequence lifetime, whichever occurs first. We note that collisions are a Poisson process, and two collisions occurring within one timestep is within the realm of possibility, though highly improbable. We estimate that one star out of 20,000 in the inner 0.1 pc may experience two collisions over a 10⁶ yr timestep. The effects of this stochastic process are likely already captured by our statistical approach to collisions.

We also simulate the effects of relaxation in our code. Over each timestep, we apply a small instantaneous velocity kick to the star, from which we calculate the new, slightly altered orbital parameters. The velocity kick is drawn from a Gaussian distribution with a standard deviation that depends on the ratio of the orbital period to the relaxation timescale for the star in question (e.g., Bradnick et al. 2017; Lu & Naoz 2019; Naoz et al. 2022; Rose et al. 2022, 2023, see Naoz et al. 2022 for the full set of equations). This prescription allows us to account for diffusion in the orbital parameters over time from interactions with the surrounding stars. Additionally, relaxation processes can drive stars onto near-radial orbits such that they become tidally disrupted by the SMBH (e.g., Magorrian & Tremaine 1999; Wang & Merritt 2004; Hopman & Alexander 2005). If a star's periapsis is driven to a value less than the tidal disruption radius, it triggers a stopping condition in the code, effectively terminating the lifetime of the star. ⁵

If a collision occurs, we must adjust the mass of the star in our sample accordingly. The collision outcome depends in part on the speed of the impact. The velocity dispersion in the nuclear star cluster exceeds 100 km s⁻¹ within about 0.1 pc of the SMBH. Heuristically, high-velocity collisions should result in mass loss from the colliding stars. In extreme cases, when the relative velocity is larger than the escape speed from the star, one might expect the collision to fully unbind the star.

Rose et al. (2023) find that the collision outcomes can generally be understood in three regimes. Near the SMBH, the velocity dispersion exceeds the escape velocity from the stars, leading to destructive collisions with high mass loss. Between about 0.01 and 0.1 pc, collisions are common and can lead to mergers. These merger products, more massive than their progenitor stars, evolve off the main sequence over a shorter lifetime. Outside of 0.1 pc, collisions are less frequent, but always lead to mergers with little mass loss.

In more detail, collision outcomes depend on the gas dynamics of the collisions themselves and are dependent on several other conditions beyond the impact velocity, such as the impact parameter and the internal structure of the stars (e.g., Lai et al. 1993; Freitag & Benz 2002; Rubin & Loeb 2011). Previous studies have leveraged hydrodynamic simulations to understand collision outcomes at high velocities for different impact parameters (e.g., Lai et al. 1993; Rauch 1999; Freitag & Benz 2002). Lai et al. (1993) and Rauch (1999), in particular, provide fitting formulae based on their results. These formulae can be easily implemented in a code such as ours to estimate the mass lost from the stars and whether or not a merger occurred (Rose et al. 2023).

In this study, we compare three different recipes for determining the collision outcome. We begin by considering a limiting case in which collisions are fully destructive. In this prescription, the occurrence of a single collision terminates the code. Once a star in our sample experiences a collision, its mass is set to zero and it is removed from the population. While unphysical, this treatment of collisions allows us to build an intuition for the underlying physics of our models. We then proceed to include a more complicated treatment of the collision outcomes. These simulations utilize either fitting formulae from Rauch (1999) or Lai et al. (1993), discussed in more detail in Rose et al. (2023). Henceforth, we refer to simulations that use fitting formulae from these studies as "Rauch99" and "Lai+93," respectively. These fitting formulae allow us to calculate the mass loss from a given collision given the mass ratio of the colliding stars, the impact parameter, which is drawn statistically, and the relative velocity. Our simulations always assume that the relative velocity is equal to the velocity dispersion, a function of distance from the SMBH given by Equation (2). As described in Rose et al. (2023), our Rauch99 simulations favor mergers with little mass loss, while collisions lead to higher mass loss in the Lai+93 simulations and mergers are less likely. Together, these prescriptions are selected to span the range of possible collision outcomes.

Collisions play a crucial role in shaping the stellar demographics where the collision timescale, t_coll, is less than the age of the population. Setting t_coll equal to t_age gives a critical radius, r_coll, within which the vast majority of the stars have collided. Outside of r_coll, a fraction of the stellar population will still experience collisions. This fraction can be estimated using t_coll/t_age (see also Figure 1 in Rose et al. 2023). In addition to population age, r_coll also depends on the steepness of the density profile. For example, the age of an old 10 Gyr population and collision timescale intersect closer to 0.04 pc for α = 1.25, compared to ∼0.1 pc for α = 1.75.

This analysis informs where in the nuclear star cluster we expect collisions to be an important process in shaping the cluster properties. Because collisions modify or destroy stars at high enough velocity, we predict that r_coll will correspond to a break in the stellar density profile. Within r_coll, collisions are an important process in determining the stellar density profile. Outside of this critical radius, collisions are rare, and the density profile is not shaped by collisions. Over time, as the age of the population increases, the inflection point will move further from the SMBH. For an old 10 Gyr population, with properties similar to the Milky Way GN, r_coll occurs at about 0.1 pc, shown in Figure 1.

We can perform a similar analysis to determine where in the cluster collisions can effectively deplete the entire supply of stars. Within about 0.01 pc of the SMBH, the velocity dispersion, given by Equation (2), exceeds the escape velocity from a Sun-like star. In this region, collisions have the potential to destroy the stars. About two-thirds of the Sun's mass is concentrated in the inner third of its radius (e.g., Christensen-Dalsgaard et al. 1996). A collision will result in high mass loss when the impact parameter is small enough that the dense cores interact (e.g., Lai et al. 1993; Rauch 1999; Rose et al. 2023). We define a characteristic timescale over which the stellar cores will collide by setting r_c from Equation (3) equal to 2 × 0.33 R_⊙. Within about 0.01 pc of the SMBH, this timescale represents the characteristic timescale over which a star will be destroyed. While we account for gravitational focusing in our calculations, the high relative velocity ensures that gravitational focusing is very weak in this regime. Figure 1 shows this timescale in green. This timescale is consistent with Rose et al. (2023), who find the time needed to deplete the stellar population within 0.01 pc to be about 1 Gyr. Similar to r_coll, we define r_dest as the radius at which the destruction timescale equals the population age. We stress that this definition is only valid in regions where the collision velocity is high enough to destroy the stars.

The break radius for any GN can be found by equating the population age and collision timescale: t_age = (n σ A)⁻¹. Similar to the Milky Way's GN, the SMBH of mass M_• dominates the gravitational potential within the sphere of influence, and the relative velocity σ can be calculated using Equation (2). The initial stellar density profile of the cluster must be calibrated to the mass of the central SMBH. Using the M–σ relation, the stellar density profile for a GN with an SMBH of arbitrary mass can be written as a power law:

$$\begin{eqnarray}&&\rho ({r}_{\bullet })=\displaystyle \frac{3-\alpha }{2\pi }\displaystyle \frac{{M}_{\bullet }}{{r}_{\bullet }^{3}}{\left(\displaystyle \frac{G{\left({M}_{0}{M}_{\bullet }\right)}^{1/2}}{{\sigma }_{0}^{2}{r}_{\bullet }}\right)}^{-3+\alpha },\end{eqnarray} \tag{ 5 }$$

where M₀ = 1.3 × 10⁸ M_⊙, σ₀ = 200 km s⁻¹, and r_• is the distance from the SMBH (Tremaine et al. 2002).

We can derive a scaling relation using r_coll = 0.09 pc for the Milky Way's GN at 10 Gyr for a Bahcall–Wolf profile (α = 1.75). For a GN with arbitrary stellar density profile slope α, the break radius can be calculated using the following:

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\mathrm{coll}} & = & {\left[1.49\times {10}^{-6}{\left(6.04\right)}^{\alpha }\displaystyle \frac{{\left(3-\alpha \right)}^{2}}{1+\alpha }\right]}^{\tfrac{1}{1\,+\,2\alpha }}\\ & & \times \ {\left(\displaystyle \frac{{t}_{\mathrm{age}}}{{10}^{10}\,\mathrm{yr}}\right)}^{\tfrac{2}{1\,+\,2\alpha }}{\left(\displaystyle \frac{{M}_{\bullet }}{4\times {10}^{6}{M}_{\odot }}\right)}^{\tfrac{\alpha }{1\,+\,2\alpha }}\mathrm{pc}.\end{array}\end{eqnarray} \tag{ 6 }$$

For values of α between 1 and 2, the first term in brackets changes by a factor ∼5.4. The following equation can serve as an approximation in this range of α:

$$\begin{eqnarray}&&{r}_{\mathrm{coll}}\sim 0.07{\left(\displaystyle \frac{{t}_{\mathrm{age}}}{{10}^{10}\,\mathrm{yr}}\right)}^{\tfrac{2}{1+2\alpha }}{\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{\odot }}\right)}^{\tfrac{\alpha }{1+2\alpha }}\,\mathrm{pc}.\end{eqnarray} \tag{ 7 }$$

If the GN also has a Bahcall–Wolf profile, i.e., α = 1.75, the full relation can instead be simplified to

$$\begin{eqnarray}&&{r}_{\mathrm{coll}}=0.09{\left(\displaystyle \frac{{t}_{\mathrm{age}}}{{10}^{10}\,\mathrm{yr}}\right)}^{0.44}{\left(\displaystyle \frac{{M}_{\bullet }}{4\times {10}^{6}\,{M}_{\odot }}\right)}^{0.39}\,\mathrm{pc}.\end{eqnarray} \tag{ 8 }$$

We use this last relation to calculate r_coll in our fiducial models, setting M_• = 4 × 10⁶ M_⊙.

Figure 2 shows the time evolution of the model population of "tracer" stars. We show four snapshots for two different simulations, both of which adopt α = 1.75 for the density profile. The top row assumes that collisions are fully destructive, while the bottom row uses a Rauch99 prescription for collision outcomes. The gray points in the plots show the stellar masses at a given time due to stellar evolution alone. As can be seen in both rows of the figure, the vertical extent of the gray points decreases with time because the time elapsed has exceeded the main-sequence lifetime of stars above a mass threshold.

Zoom In Zoom Out Reset image size

The red points show the stellar masses from our simulations, which account for the effects of collisions. We mark the location of r_coll with a dashed red vertical line in each snapshot, the distance within which the majority of the stars have collided. Calculated from Equation (8), r_coll sweeps out over time.

Depending on the collision prescription and the impact velocity, collisions can either add or remove mass from the impacting stars (as described in detail by Rose et al. 2023). In the top row, collisions carve out a starless region over time. The bottom row, on the other hand, incorporates a more complete and complex treatment of collision outcomes. Thus, the region inside r_coll contains both collisionally merged, more massive stars and collisionally stripped, less massive stars as compared to their progenitors.

We consider the evolution of the stellar density profile over time due to collisions, which can work to either destroy or merge stars. To generate a density profile, we slice up each snapshot, including those depicted in Figure 2, in distance. In each annulus with width δ r that lies distance r from the SMBH, we divide the total mass in the red points by the total mass in the gray points, which tells us the fractional change in mass at that distance at a particular time. We then convolve those fractions with the original density profile, a cusp with index α, of the stars to obtain a collisionally modified profile.

We juxtapose the evolving density profiles of three simulations in Figure 3. Each simulation uses the same initial value for α, 1.75. They use different collision outcome prescriptions noted in the plot titles. The upper panel of each plot shows the stellar density as a function of distance from the SMBH at a particular time. As indicated by the color bar on the right, the redder curves correspond to older populations. The bottom panel shows the fractional change in stellar density due to collisions compared to the expected value. The black line shows the original, reference density profile.

Zoom In Zoom Out Reset image size

In each case, collisions deplete the stellar mass near the SMBH. This process causes the density profile to flatten within r_coll. Over time, this r_coll moves further from the SMBH, shifting the break radius in the density profile. The fully destructive case in the left column of Figure 3 provides the clearest example of the break radius sweeping outward over time: the bluest density profiles diverge from the unmodified profile at smaller radii. Equating timescales suggests a break radius of ∼0.1 pc for a ∼10 Gyr population with these initial conditions (Equation (8)). The break radius of the fully destructive model falls slightly outside of this radius, closer to 0.3 pc for a 7 Gyr old population, because collisions still affect a fraction of the population outside of r_coll leading to mass loss for some percentage of the stars (see Figure 1 in Rose et al. 2023), causing break radii to be gradual rollovers rather than abrupt transitions.

A more realistic mass loss and merger prescription gives a less distinct break radius, as seen in the second two columns of Figure 3. In these prescriptions, a single collision results in comparatively smaller fractional changes in total stellar mass. Compared to the Rauch99 model, the Lai+93 density curves show greater flattening within 0.05 pc. The greater flattening occurs because fitting formulae from the corresponding study (Lai et al. 1993) give a higher fractional mass loss per collision in this region compared to Rauch (1999). Mergers can also contribute to a break in the density profile. While mergers result in fractional mass loss, generally between 5% and 10% (see Figure 2 in Rose et al. 2023), they also hasten the evolution of the stars off the main sequence. As a result, there is less mass in main-sequence stars at late times within 0.1 pc compared to the outer region of the nuclear star cluster.

We also examine the dependence of the break radius on the initial density profile assumed for the population. We test three values of α: 1.25, 1.5, and 1.75. Figure 4 juxtaposes the evolving density profiles for three simulations with these different values of α. These simulations use the same collision outcome prescription, Rauch99. We find that a break in the density profile is ubiquitous, but its location depends on α, shifting from ∼0.1 pc for α = 1.75 to ∼0.04 pc for α = 1.25. The steeper the initial profile, the greater the fractional change in density, as seen in the bottom panels of the figure. To illustrate these trends, we indicate the point at which the density at 7 Gyr is halved compared to the expected value using a gray vertical line in each panel. It shifts inward with increasing α. Regardless of initial conditions, we find that the density profile is always preserved, i.e., unmodified by collisions, outside about 0.2 pc. Observations of this region can therefore be used to constrain the underlying, original density distribution.

Zoom In Zoom Out Reset image size

Our models assume the nuclear star cluster holds a single stellar population, originating from the same star formation episode, that has relaxed onto some equilibrium density distribution, to which we apply the effects of collisions. In this section, we address these assumptions. Recent star formation episodes, for example, can contaminate the cusp as the younger population may not have had sufficient time to relax onto a cusp (e.g., Baumgardt et al. 2018; Schödel et al. 2018). Young stars are present in the Galactic center (e.g., Paumard et al. 2006; Lu et al. 2009; Schödel et al. 2018, 2020), though observations suggest that an older (> few Gyr) population dominates the central cluster (e.g., Nogueras-Lara et al. 2019; Schödel et al. 2020; Chen et al. 2023; Schödel et al. 2023). Chen et al. (2023) find that about 90% of the stars are 5 Gyr old. Future work will take into account a multiple component population for the Milky Way's nuclear star cluster. However, when applying Equation (6) to other GN, it must be assumed that their stellar populations are similarly dominated by a single, relaxed stellar population. Additionally, we have presented examples of stellar density profiles in accordance with various theoretical predictions and observational studies for the Milky Way (e.g., Gallego-Cano et al. 2018; Schödel et al. 2018; Linial & Sari 2022). However, the density profile may be even shallower and evolving over ≳Gyr timescales (e.g., Antonini 2014; Panamarev et al. 2019). The break radius may still be estimated using Equation (6). r_coll varies slowly with changing α (see Section 3.3 and Equation (7)), so the equation may still be used for evolving stellar density profiles so long as α does not change by more than 1 over the time interval, t_age, selected.

As we discuss in the previous section, mass loss due to collisions together with the accelerated evolution of merger products can lead to a discernible flattening of the stellar density profile over time. In this section, we briefly address the effect of collisions on the stellar luminosity. Mergers are the typical outcome of collisions outside of about 0.01 pc (Rose et al. 2023). Consistent with our previous treatment of the merger products, we assume that the properties of the merged star are similar to those of a typical main-sequence star of the same mass (e.g., Leiner et al. 2019). In this case, the luminosity of the stars should follow a simple mass–luminosity relation: L ∝ M^2.3 for M < 0.43 M_⊙, ∝M⁴ for Sun-like stars, and ∝M^3.5 for M > 2 M_⊙ (e.g., Duric 2003; Salaris & Cassisi 2005). Additionally, we can calculate the peak wavelength for the stars in our sample with Wien's displacement law, in which we use the peak wavelength as a proxy for color. The general assumption here is that we consider the properties of these objects after a period of thermal relaxation, not when they are still cooling following collisional shock-heating (see discussion of thermal versus collision timescale in Rose et al. 2023).

We illustrate the effects of collisions on color and luminosity in Figure 5. This figure presents a snapshot of our sample stars at 4.64 × 10⁹ yr. The left panel shows the expected population in the absence of collisions, equivalent to the gray points in Figure 2, while the right side shows the simulated population. We plot the peak wavelength of the stars versus their distance from the SMBH, both on a logarithmic scale. The y-axis is inverted so that more massive (bluer) stars are higher, facilitating a direct comparison with the mass versus radius plot in Figure 2. Additionally, we color-code the points by bolometric luminosity, calculated using the relations above.

Zoom In Zoom Out Reset image size

As can be seen in the figure, mergers produce brighter and bluer stars than expected for a population of that age, like "blue stragglers" in a star cluster population (Leiner et al. 2019). We note that in addition to younger-seeming collision products, there may also be young stars from recent star formation episodes residing in the nuclear star cluster (e.g., Lu et al. 2009; Neumayer et al. 2020). Collisionally stripped stars are also present in Figure 5. These stars have undergone one or more high-speed collisions, leading to mass loss. Their new, lower masses would suggest that they are redder and less luminous. However, while we have plotted them here based on this assumption, their appearance is highly uncertain. Models of stripped stars in binary systems may provide clues as to their appearance, suggesting they are in fact more luminous than their progenitor stars (Götberg et al. 2018). We therefore distinguish them with gold outlines in the figure.

Previously, we showed that collisions always produce a break in the stellar mass density profile. The break occurs because collisions can only ever reduce the mass contained in the stars; it cannot increase the mass. Luminosity, however, scales as mass to the 3.5 power. We therefore do not expect collisions to necessarily produce a decrease in the bolometric luminosity profile of the cluster within r_coll. The bluer, brighter merger products may in fact outshine the other stars in their vicinity. Figure 5 shows that the colors and luminosities of stars are affected by collisions in a way that varies systematically with radius. Luminous, blue merger remnants exist from 0.01 to 1 pc, around r_coll ∼ 0.1 pc. Stripped stars from the highest-velocity collisions preferentially lie within 0.1 pc. Future work may examine the evolution of the luminosity profile, both bolometric and in specific bands, by leveraging comprehensive hydrodynamic and stellar evolution simulations in order to understand the dynamical and thermal evolution of merger products.