Stellar Collisions in Galactic Nuclei: Impact on Destructive Events Near a Supermassive Black Hole

Centers of galaxies host both a supermassive black hole and a dense stellar cluster. Such an environment should lead to stellar collisions, possibly at very high velocities so that the total energy involved is of the same order as supernovae explosions. We present a simplified numerical analysis of the destructive stellar collision rate in a cluster similar to that of the Milky Way. The analysis includes an effective average two-body relaxation Monte-Carlo scheme and general relativistic effects, as used by Sari and Fragione (2019), to which we added explicit tracking of local probabilities for stellar collisions. We also consider stars which are injected into the stellar cluster after being disrupted from a binary system by the supermassive black hole. Such stars are captured in the vicinity of the black hole and enhance the expected collision rate. In our results we examine the rate and energetic distribution function of high velocity stellar collisions, and compare them self-consistently with the other destructive processes which occur in the galactic center, namely tidal disruptions and extreme mass ratio inspirals.


INTRODUCTION
The presence of supermassive black holes in the centers of galaxies appears to be ubiquitous in the universe (Kormendy & Ho 2013).Such black holes regulate galaxy formation (Heckmn & Best 2014), and drive radiative emission in active galactic nuclei (Rees 1984;Fabian 2012).A central black hole also determines many of the properties of the dense stellar cluster that surrounds it near the galactic center.The structure and dynamics of the cluster are a combined result of the black hole dominating the potential, along with a multitude of gravitational interactions amongst the stars (Rauch & Tremaine 1996;Hopman & Alexander 2005;Merritt 2013;Alexander 2017;Fouvry et al. 2022).This combination tends to generate a particular spatial density profile in the stellar cluster, which evolves dynamically through gravitational relaxation processes.
Corresponding author: Shmuel Balberg shmuel.balberg@mail.huji.ac.ilA unique aspect of the dense stellar cluster is the destructive events which its stars experience.Much attention has been devoted to tidal disruption events (TDEs), when a star approaches the supermassive black hole (SMBH) at a distance R T , where the tidal force of the black hole overpowers the gravitational binding of the star.Such an event is possible when R T lies outside of the event horizon, and possibly leads to a highly luminous transient, as the stellar matter interacts with itself while accreting onto the black hole (Rees 1988;Gezari 2021).
As noted before by several authors (see, e.g.Alexander (2017) for a review), there are actually two channels that lead a star to tidal disruption.First, stars throughout the cluster can be gravitationally scattered by other stars into highly eccentric orbits, whose periapse lies below R T .Such scatterings can occur at any distance from the SMBH, including stars at the radius of influence, R h , which are usually the most abundant.Hereafter we will refer to this channel as a TDE.The second channel is due to general relativistic effects, in which gravitational wave (GW) losses cause single stars to inspiral gradually toward the SMBH.Initially GW losses tend to circularize the orbit while keeping an almost fixed periapse, followed by a gradual inspiral toward R T along roughly circular orbits.While technically this inspiral also terminates with the star being tidally disrupted, the process occurs over multiple orbits causing a more "gentle" disruption.This latter channel is referred to as a main-sequence extreme mass ratio inspiral (EMRI), and is mostly studied in the context of potential sources of observable GW emission (along with EMRIs of compact objects which cannot be tidally disrupted and therefore spiral all the way to the event horizon).However, given that a significant fraction of the star's mass may shed prior to disruption, optical transients are also possible (Linial & Sari 2023), especially if more than one star can undergo an EMRI simultaneously (Metzger & Stone 2017;Metzger et al. 2022).
The rates of TDEs and EMRIs depend on the properties of the SMBH and the stellar cluster.In both cases, the rate is determined by gravitational scattering among the stars, which repopulate orbits that terminate in disruption.In general, the typical estimated TDE rate for a Milky Way-like SMBH and star cluster is of order 10 −5 yr −1 .The EMRI rate (of main-sequence stars) in a fully relaxed cluster with Milky Way properties is estimated to be about 1% of the TDE rate, but as shown by Sari & Fargione (2019), this rate is significantly enhanced when disrupted binaries are taken into account (see below).
Relatively little attention has been given to destructive stellar collisions, which are also possible in a dense cluster surrounding an SMBH.Physical collisions in the dense stellar cluster should actually be common, but mostly at large distances from the black hole, with velocities which are smaller than the stars' escape velocities.Such collisions are more likely to result in one or both stars surviving or merging (Freitag & Benz 2005).Mass loss, mass transfer, and mergers in such collisions should reshape the mass function of the stellar population (Rose et al. 2023).In particular, collisions are suggested as a possible formation mechanism for young, high-mass stars, such as S-stars observed surrounding Sagittarius A (Fragione & Antonini 2019).However, closer to the SMBH, where the orbital velocities reach thousands of kilometers per second, collisions are completely destructive, and could led to a variety of potential observational consequences (Amaro-Seoane 2023).
A fast, head-on collision of this sort will result in hot expanding gas, implying that it should produce some form of an optical flare.In fact, at the higher end of the velocity range, collisions between sun-like stars at 10,000 km S −1 contain kinetic energies of 10 51 erg or more, which conform with supernova-like energy output.Accordingly, Balberg et al. (2013) suggested that such hypervelocity collisions should be considered as another, rare, form of supernova.A rough estimate suggests that the resulting optical flare should indeed be "supernovalike" over a time scale of order a few days (and decay quickly given that no radioactive isotopes are produced).A second luminous transient might be possible if there is significant fast accretion of the debris onto the SMBH Clearly destructive collisions (DCs) are inevitable, but at what rate?How do they compare, and possibly compete, with TDEs and EMRIs as stellar destruction scenarios near the SMBH?The answers depend on the steady-state distribution of stars in the cluster surrounding the SMBH.In turn, this distribution strongly depends on the mechanisms that supply stars to orbits which are relevant for destruction in the different scenarios.
In this work, we present a simplified numerical study of the expected rates of DCs in a "Milky Way-like" galactic center.We approximate the stellar cluster as being composed of identical mass m stars orbiting the SMBH, and track the evolution of the cluster with a simplified Monte Carlo simulation.This scheme was originally presented by Sari & Fargione (2019) for TDEs and EMRIs, to which we added DCs.In the simulations we include stars which originated from disrupted tight binaries by the SMBH (Hills 1988), and captured in highly eccentric orbits with relatively small (1 − 100R T ) periapses.The existence of such captured stars is inferred from observations of hypervelocity stars ejected from the Milky Way (Brown 2015), which are believed to be companions that became unbound during the binary disruption.As shown by Balberg et al. (2013), these captured stars are prime candidates for DCs, and in fact could make their rate comparable to that of TDEs.Stars from disrupted binaries obviously also affect the TDE and EMRI rates (Sari & Fargione 2019), and the main goal of the present work is to study all three in a selfconsistent calculation.
The structure of this manuscript is as follows.In Section 2, we review the physics of two-body relaxation, GW emission, and collisions and evaluate their time scales in a schematic (power-law density) cluster.The equations and setup of our Monte Carlo simulation, following Sari & Fargione (2019) and the introduction of the possibility for collisions, are presented in Section 3. In Section 4 we examine the results of the simulations, and most notably the rates of the different destruction scenarios when no stars from disrupted binaries are included.Such stars are added to the simulations in Sec-tion 5, where we demonstrate their dramatic enhancement of the DC rate.Finally, in Section 6 , we discuss the implications of the results and our corresponding conclusions.

PRINCIPAL PROCESSES IN THE DENSE STELLAR CUSP
A stellar cluster surrounding an SMBH is expected to form a dense cusp, with the density monotonically increasing toward the center.The conventional approach is to identify the SMBH radius of influence, R h , at which the total enclosed mass of the stars is equal to the SMBH mass, M • .Inside of R h we approximate that the SMBH dominates the gravitational potential, and that at every instant the motion of each individual star is consistent with a Keplerian orbit in this potential.In this work we set all the stars in the cusp to have identical mass m, so the number of stars enclosed by the radius of influence is simply The orbits of the stars do evolve over time with respect to a stationary Keplerian orbit.The two main driving processes are gravitational scattering among the stars and GW emission.Both processes can be estimated by approximating the stars as point masses.
Stellar collisions introduce an additional mechanism which can change both the mass and the number of stars.Low-velocity collisions (and even near misses) far from the SMBH can result in mass transfer, mass loss, and mergers, while hypervelocity collisions at small distances from the SMBH remove stars through complete destruction in a singe event.The prospects and physics of collisions depend on the sizes of the stars, which we treat here as identical with radius R ⋆ .

Gravitational scattering and two-body relaxation
While the SMBH dominates the gravitational potential, stars also interact with each other.Any given star experiences mostly weak, small-angle scatterings with other stars, but rare strong scatterings at small impact parameters occur as well, and have a similar effect on the total accumulated changes in the specific energy, E, and the specific angular momentum, J (Lightman & Shapiro 1977).Under the simplest assumption that all scatterings are random, the star's orbital parameters change diffusively in phase space.A star orbiting the SMBH with a semi-major-axis (sma) r will change its specific angular momentum by order unity over a two-body relaxation time, given by (Bahcall & Wolf 1976, 1977) (1) Here P (r) is the orbital period of a Keplerian orbit around the SMBH, N (r) is the number of stars in the proximity of the radius r and ln Λ is the resulting Coulomb logarithm.Assuming that most stars have roughly circular orbits, N (r) is approximately the number of stars enclosed by radius r.The time scale for order-of-unity diffusion in energy space is somewhat longer than T 2B (r), so stars generally change their periapse prior to changing their smas significantly.This two-body relaxation time allows us to derive a rough estimate for the TDE rate (Sari & Fargione 2019).Assuming that the cluster is dominated by stars with an sma r ∼ R h and that at R h the velocity distribution is isotropic, TDEs reflect the rate at which stars at the radius of influence are scattered into orbits with angular momentum A system with Milky Way-like parameters will tend to an "empty-losscone" scenario, in which there is on average less than one star in an orbit that terminates in tidal disruption.Hence the TDE rate is set by stars evolving from a circular orbit to a highly eccentric one that will "plunge" below R T .This occurs as a diffusive process in twodimensional angular momentum space, originating from the circular angular momentum J c (R h ) = √ GM • R h and terminating at the loss-cone value J LC .The TDE rate should be of order since by construction N (R h ) = M • /m, and we assume that ln(J c (R h )/J LC ) and ln Λ have similar values (of order 10).For the Milky Way P (R h ) = 2π(R 3 h /GM • ) 1/2 ≈ 1.3 × 10 5 yr, implying a TDE rate of about 7.6 × 10 −6 yr −1 .
For our purpose of studying the steady-state rates of all destructive events which occur close to the SMBH, we also consider the evolution of orbits.A star on an eccentric orbit with sma r and periapse distance r p will change its specific angular momentum by order unity over a time scale smaller than T 2B (r) by a factor of r p /r ≈ (J/J c (r)) 2 .The eccentric orbital relaxation time is therefore (Binney & Tremaine 1987) which is appropriate if the density profile is not too steep, so that stars at r dominate the scatterings (this assumption applies as long as the density of the stars, n(r), is shallower than n(r) ∼ r −4 , which we shall assume henceforth).
In this work we follow the standard low-order approximation of gravitational dynamics between stars, and limit ourselves to random two-body scatterings.In reality the dynamics between stars are subject to other relaxation mechanisms, commonly referred to as resonant relaxation (Rauch & Tremaine 1996;Alexander 2017).These arise from residual torques between the stars, which lead to nondiffusive evolution of the angular momenta, and over time scales which are significantly shorter than T 2B (Kocsis & Tremaine 2011, 2015).On the other hand, the effect of resonant relaxation is restricted by the timescales over which a star experiences coherent fluctuations in the gravitational potential (Hopman & Alexander 2006).Detailed analysis (Bar-Or & Alexander 2016) shows that random two-body scattering are sufficient to suppress coherent effects far from the SMBH, and the same occurs close to the SMBH due to general relativistic effects.As a result, resonant relaxation in a Milky Way-like cusp is limited to intermediate distances from the SMBH.Since TDEs are mostly the result of stars being scattered into eccentric orbits close to R h , and EMRIs arise from general relativistic effects close to the SMBH, resonant relaxation was found to have little impact on their rates.Our focus is on collisions in regions close to the SMBH, so we expect that neglecting resonant relaxation is a legitimate approximation for our analysis as well.

Orbital evolution through GW emission
Gravitational wave emission by the stars in the background field of the SMBH is a dissipation mechanism.It drives an orbit toward a smaller sma, along with a slower decay of the periapse distance (Peters 1964).As a result, GW emission will tend to circularize eccentric orbits, and then drive a gradual inspiral of roughly circular orbits to the innermost stable relativistic orbit (ISCO).The time scale for a full inspiral for an orbit initially with sma r and periapse distance r p is (Peters 1964;Hopman & Alexander 2005) where R S = 2GM • /c 2 is the Schwarzchild radius.This result applies when T GW (r, r p ) is longer than the star's orbital period, P (r).This condition is easily satisfied for any r < R h (Sari & Fargione 2019), and while it is only marginally satisfied at r ≈ R h , the general relativistic effects there are small enough to be ignored in relevant calculations.
As noted by Sari & Fargione (2019), there exists a region in the r − r p plane for which T GW (r, r p ) < T J 2B (r, r p ), depending, of course, on the steady-state stellar density profile, N (r).Orbits in this region will evolve toward an EMRI, rather than toward a TDE.Specifically, if the steady state follows the Bahcall & Wolf (1976) power law of N (r) ∼ r 5/4 (hereafter a "BW profile"), the r p (r) relation for which T GW (r, r p ) = T J 2B (r, r p ) has the form For stars with sma r and a periapse r p smaller than that of equation 5, the orbit will evolve into a series of quasi-circular orbits, and terminate as an EMRI.This analysis yields an estimate of the EMRI rate (Sari & Fargione 2019).The largest influx of stars into combinations of (r, r p ) at which GW dominates the orbital evolution is from eccentric orbits with an sma r 0 which satisfies equation 5 with r p = R T .This sma is r 0 = R h (R S /R T ) 2 , and in a similar fashion to equation 2 (again, assuming a BW profile of N (r) ∼ r 5/4 ) (6) For the Milky Way SMBH and sun-like stars, and assuming a BW density profile this results in R EM RI ≈ 10 −2 R T DE .

Destructive stellar collisions
Stars will physically collide if their centers come within a distance f R R ⋆ of each other.Henceforth we track collisions only in the inner part of the cluster, where the orbital velocities are greater than the escape velocity of the stars.Correspondingly, we ignore gravitational focusing, so that f R ≤ 2. In the following we consider f R = 0 (no collisions), f R = 1 and f R = 2.The latter is an upper limit and probably an overestimate, since grazing collisions, even at hypervelocities, may not be completely disruptive.
Stars on eccentric orbits sample different stellar densities along their path.In density profiles steeper than n(r) ∼ r −1 the optical depth for collisions is dominated by stars near the periapse (Sari & Fargione 2019) and the typical collision time for a star with orbital parameters (r, r p ) is of order Although they did not include collisions in their work, Sari & Fargione (2019) did point out that the collisional time scale (equation 7) is usually shorter than both the two-body and the GW timescales (for Milky Way parameters and assuming a BW cusp).Quantitatively, they assessed that collisions dominate the evolution for practically all orbits with r < 10 −2 R h .Furthermore, consider the rate of hypervelocity collisions in a power-law density profile, n(r) = Ar −α .Integrating over all radii in a range r min to r max , the DC rate is where v(r) ≈ (GM • /r) 1/2 is the typical velocity at radius r.The constant A is determined by the total amount of stars in the cluster (see also Amaro-Seoane ( 2023)).Equation 9applies for 5/4 < α < 3 so that N is dominated by stars in the vicinity of R h , but the collisions are dominated by the vicinity of r min .For a BW profile with α = 7/4, Milky Way and sun-like values, and setting r min = R T , we find a dramatic rate of R DC ≈ 10 −3 yr −1 (for f R = 1).This result is unphysical and clearly indicates that collisions between stars will tend to deplete the inner part of the cusp.Depletion will dominate out to a radius R col , where the orbital velocity is approximately equal to the escape velocity of the stars v(r) ≈ GM • /R col = v esc .For the Milky Way and sun-like parameters this is roughly also the radius where T 2B = T col (Sari & Fargione 2019), so we expect that that the entire density profile will include a fully relaxed cusp (with a BW profile) which extends from R h into R col , where DCs take over and the number of stars plummets.The collision rate is therefore determined by the BW density at where we used the BW relation For the Milky Way and sun-like values we repeatedly assume, this total estimate is R DC ≈ 10 −7 yr −1 .Notably, if DCs do indeed deplete the inner part of the cusp, stars are unlikely to evolve through GW emission in multiple inspiraling orbits all the way to R T .The EMRI rate must therefore be reduced by the possibility of DCs.The TDE rate may also be affected, since some stars on eccentric orbits and r ∼ R h could suffer collisions before finally being scattered into J < J LC , but since only several orbits are required prior to the final plunge below R T , this effect should be small.In any case, it is clear that collisions must be accounted for in any analysis of the rates of destructive processes in stellar cusps.

THE MONTE CARLO SIMULATION
Despite continuous progress in N -body simulations (see, e.g., Arnold et al. (2022)), tracking ∼ 10 6 stars over long time scales is still impractical.This is especially true in our case, in which we need to follow stars close to the SMBH accurately where the periods are very short, while allowing the entire system to evolve over a cluster relaxation time.
Instead, we opt for a simplified Monte Carlo algorithm.Each star is tracked individually in (E, J) space, but the cumulative effect of other stars is averaged out according to the instantaneous density profile.The basics of such an algorithm were originally suggested by Henon (1971), and here we follow the technical approach presented by Fragione & Sari (2018) and Sari & Fargione (2019), to which we add a consistent treatment of stellar collisions.

Setup
We consider a Milky Way-like cluster where M • = 4 × 10 6 M ⊙ , which is surrounded by a cluster of N = 4×10 6 identical stars with mass m = 1 M ⊙ .The cluster is assumed to extend up to a radius of influence of where σ is the measured velocity dispersion external to the radius of influence (Merritt & Ferrarese 2001).Stars occupy the range between the tidal radius R T = 1.5 × 10 13 cm (1 au) and R h .At the beginning of the simulation, all stars are assigned values of the sma, r, and specific angular momentum, J.The sma is drawn randomly from a weighted distribution that produces a BW profile of n(r) ∼ r −7/4 , and the specific angular momentum is set by uniformly sampling (J/J C (r)) 2 over the range 0 to J C (r) = √ GM r (the circular angular momentum given the star's sma, r).The latter is essentially a thermal distribution, which generates an isotropic velocity field.Once r and J are set, the specific energy, E, period, P , eccentricity, e, and periapse, r p , are determined for each star assuming Keplerian orbits (12)

Evolution by effective two-body scatterings
Even though r and r p are continuous and are stored individually for each star, it is numerically convenient to divide the (r, r p ) plane into evenly spaced logarithmic bins.We set r i+1 = 2r i , r p,j+1 = 2r p,j , so i and j serve as indices of the resulting two-dimensional grid.Following Sari & Fargione (2019), we track at each time step the number of stars in any particular bin, N (i, j), and the total number of stars with an sma enclosed in a specific bin, regardless of periapse For each bin (i, j), we calculate the momentary eccentricity-dependent time step, where is a universal parameter of the system.We set ln Λ = 10, and the control parameter γ = 1.5 that takes into account the bin size (Fragione & Sari 2018), and also f T = 0.1 as the fraction of the relaxation time allowed in a single evolutionary time step.Time in the simulation is advanced in discrete time steps, ∆t, which are determined at each time t by the fastest evolving grid point (i, j) Updating all stars at every time step is very expensive numerically.We use the scheme suggested by Sari & Fargione (2019) in which for each grid point (i, j) we store the last time it was updated, t L (i, j), and update it again only when the current time, t, satisfies If a grid point (i, j) is due to be updated, each individual star in that bin is allowed to evolve its orbital parameters.For each star, k, we identify the range of distances from the SMBH that it samples along its orbit, which correspond in the i-th component of the grid to the range [W 1 , W 2 ].For each bin W we assume that the inspected star follows a random walk by a series of scatterings with the stars in that bin.The cumulative effect has an average impact parameter where P * k is the period of the inspected star.The impact parameter is converted to typical changes in the specific energy and angular momentum due to the stars at grid cell i = W , with and where v(W ) = GM • /r W is the characteristic velocity at r W . Finally, the energy and the angular momentum of the star are updated with respect to their current values of E old and J old to with 0 ≤ χ < 2π and 0 ≤ Φ < 2π randomly drawn from a uniform distribution.The star's orbital properties, r,P , r p , and e, are then updated consistently.This process is repeated for the entire range [W 1 , W 2 ] for the inspected star, and for all stars in the (i, j) bin that is being advanced in time.
It is worthwhile to point out explicitly the codependency between the choice of grid separation (r i+1 /r i ) and the approximations used in Equations (18-20).The latter assumes that two-body interactions are dominated by multiple small scatterings at large distances.These assumptions break down for thin spatial bins (r i+1 /r i ≪ 2).Not only will the assumption N (r) ≫ 1 be affected, but also the identification implied in equation 18, that a star contributes to scattering other stars only in the bin that includes its sma.This identification is appropriate, of course, for roughly circular orbits, but with r i+1 /r i = 2 it is actually consistent even with eccentric orbits and an anisotropic velocity distribution: practically all stars do indeed spend most of their period in a single spatial bin, and can be approximated as such in terms of their contribution to the average gravitational scattering.This last point is especially relevant to stars which are captured following a binary disruption (see below).In summary, given the nature of the approximations, increasing the spatial resolution does not necessarily improve the accuracy of the results.These shortcomings are best addressed by direct N -body simulations of the problem, which we report separately (G.Yassur and S. Balberg 2023, in preparation).

Evolution by GW emission
In each time step, stars evolve their orbits due to GW emission.The changes in energy and angular momentum are calculated with standard formulae (Peters 1964;Hopman & Alexander 2005) where and f (e) = 1 + (73/24)e 2 + (37/96)e 4 (1 + e) 7/2 (27) For simplicity, we follow Sari & Fargione (2019) and limit GW evolution only to stars with r p ≤ r p,GR , where we set r p,GR = 0.1R h .

TDEs, EMRIs, and Ejections
When the new angular momentum of a star corresponds to an orbit with r p ≤ R T , it is considered to be tidally destroyed by the SMBH.If this occurs following an effective two-body step (equation 22) the event is registered as a TDE, and if following a GW step (equation 24) the event is registered as an EMRI.
Since our goal is to study the steady-state cusp profile and rates of destructive events, we set to maintain a constant number of stars.We assume that every destroyed star is replaced by a new one that has been scattered into the cusp from the r > R h cluster, which obviously is not accounted for in our simulation.When stars from disrupted binaries are not included (see below), replacement stars are randomly generated in the outer bin of r (weighted to maintain an n(r) ∼ r −7/4 profile in this bin), and again with a thermal distribution of eccentricities.
Another possible outcome of an effective two-body step is that a star scatters to energy E < GM • /2R h , indicating that that star was ejected from the r ≤ R h cusp.
Again, attempting to maintain a steady-state cusp, we return the star to its previous (E old , J old ) quantities.We essentially assume that every ejected star is replaced by a star with identical properties that has been scattered into the cusp from the r > R h cluster.

Destructive collisions
Our principal innovation in this work is the inclusion of destructive two-star collisions in the inner part of the cusp.For the purposes of the simulation, we set that DCs are possible at distances smaller than R col = 2.4 × 10 16 cm (∼ 8 mpc), where the typical kinetic energy of a single star is about 10 49 erg, exceeding the binding energy of a Sun-like star.We thus neglect destruction through multiple low-velocity collisions at larger distances (Rose et al. 2023).These will certainly affect the stellar population, but are not relevant as sources of an energetic optical flare in a single event.
Our procedure for tracking collisions is as follows.We divide the radial positions from the SMBH into bins, and for convenience we use the same grid spacing as for the smas.Physical space is then divided into spherical shells with volume V i = 4π(r 3 i+1 − r 3 i )/3.For each star k that has an orbit with r p < R col , we calculate the probability p * k (i) to find the star in the volumetric bin i.The simplest cases are, p * k (i) = 1, if the orbit is entirely contained in a single spherical shell, or p * k = 0, if the star does not enter the volumetric bin i.For an orbit that crosses in and out of a bin, the probability is calculated through the fraction of time that the star spends in the shell out of the entire orbit.Denoting the time a star spends inside the spherical shell i as ∆t * k (i), the probability of finding the star in the i spherical shell is The integral in Equation ( 29) includes v r * (x), the radial component of the star's velocity when at x, which for Keplerian motion with with sma r * k and eccentricity e is .
(30) The limits on the integral in Equation ( 29) can be x min * k (i) = r i and x max * k (i) = r i+1 if the orbit passes through both surfaces of the shell (i), or x min * k (i) = r p * k or x max * k (i) = r a * k if the star only partially penetrates the shell, where r p * k (r a * k ) is the star's periapse (apoapse).
The vector p * k is converted into an effective, timeaveraged, density of stars in each spherical shell This quantity is used to calculate an effective optical depth for collisions in the ith shell.We denote the total distance the k-th star travels though the shell in a single orbit as ∆x * k (i).It may vary from 2πr * k for circular orbits to 2(x max * k (i) − x min * k (i)) for very eccentric orbits.The optical depth the kth star experiences in the ith spherical shell is then We remark that in Equation ( 32) the weight of the specific star k is subtracted from the effective density, in order to remove the fictitious option of the star colliding with itself.Finally, we calculate the probability that the kth star will collide with another star during a single orbit is Note that the factor 1/2 in the left-hand part of Equation 33 is applied as an effective correction, appropriate for small τ , in order to avoid double counting each possible collision when the vector p col, * k is compiled over all stars.The probabilities p col, * k are typically minuscule, given the very small size of the stars.The values of p * k (i) and p col, * k are calculated at the beginning of each evolutionary time step.However, applying them in the calculation presents technical challenges, due to the very wide range of time scales involved.First, the periods of the candidate stars for collisions vary over more than four orders of magnitude.Second, in our evolutionary scheme the time step size, ∆t, is set by the stars with the shortest T J 2B (equation 14).These are typically the stars with sma r ≈ R h and the highest eccentricities.The evolutionary time step is therefore of order (R 3 h /(GM • )) 1/2 , which is several orders of magnitude greater than the orbital periods of the inner stars.Tracking collisions per orbit per star is computationally prohibitive.
We overcome this complexity by estimating the probability for a collision per star over multiple orbits.The probability that a star will not suffer a collision over the entire evolutionary time step ∆t is Even though ∆t/P * k is large, p col, * k is typically very small.The result, p nocol, * k (∆t), is almost always very close to unity, and is readily evaluated.
Equation 34 includes a nontrivial assumption, which is that each star randomizes the orientation of its orbit over a single period.Obviously this assumption is incorrect in Newtonian dynamics, but is more realistic when considering relativistic precession, which is relevant in the inner cusp close to the SMBH, where we attempt to estimate the collision rate.A single simulation time step is much longer than an orbital period (and both are much shorter timescales than T GW which evolves the orbit in terms of (r, r p )), so randomization appears to be a crude, but reasonable, assumption.We examine this point separately in specific N-body simulations (G.Yassur and S. Balberg, in preparation).Given the final vector p nocol, * k (∆t), we scan at the beginning of each evolutionary time step all the stars with r p ≤ R col .For each star we draw from a random number 0 < p 0 < 1, and in the rare cases in which p 0 < 1 − p nocol, * k (∆t) the star is decreed as having experienced a collision.If this happens, the simulation determines the spherical shell where the collision happened, by drawing from the weighted distribution of exp[−τ * k (i)] (Equation 32).The partner star is determined by drawing from the weighted distribution of p * k (i) (Equation 29) in the determined ith shell.The collision is registered, along with its location in terms of the index i of the spherical shell, and the two collided stars are removed from the simulation.Again, since we seek a steady state, the removed stars are replaced in a similar fashion to those lost in TDEs and EMRIs, as described in section 3.4.
We comment that in principle a new star can itself have a nonzero probability for collision.This is a very rare possibility, but in order to account for it we examine a secondary step, in which we calculate p col * k for the replacement stars, and check for the probability of them colliding with other relevant stars.This is done for a shorter time ∆t ′ < ∆t, drawn exponentially from 0 < ∆t ′ /∆t < 1 (since we do not compute the actual time at which the first collision took place).We find that this precaution is basically meaningless when stars from disrupted binaries are ignored, while when they are included (Section 5) it does generate a small correction to the results.
It is noteworthy that an initial setup according to the BW profile does create relatively large collision probabilities.For such a profile several collisions occur during a time of order P (R h ).Obviously, this is an artifact of the steep density gradient of the BW profile, which does not account for collisions.To remedy this conflict, we processes the initial stellar distribution setup (Section 3.1) with the collision algorithm described above once for time step ∆t 0 prior to actually starting a timedependent evolution.We set ∆t 0 = (R 3 h /GM • ) 1/2 and remove and replace stars that are determined to have collided (typically of order 10, depending on the choice of f R ).These collisions are not included in the total count of collisions which occur during the simulation.

RESULTS PART I: DESTRUCTIVE RATES WITHOUT STARS FROM DISRUPTED BINARIES
We begin by examining the impact of stellar collisions on the rates of destructive events and the cusp structure when stars from disrupted binaries are not considered.The significance of adding such stars is discussed below in Section 5.
In figure 1 we show the calculated TDE, EMRI, and DC rates as a function of time.We ran the simulations for approximately one cluster relaxation time, T h (equation 1 for r = R h ).Notably, the system roughly settles on a steady state after a few tenths of T h (since the simulation is initialized with a BW n(r) ∝ r −7/4 cusp that does not account for GW emission or stellar destruction).
When collisions are not accounted for, we generally reproduce the qualitative and quantitative results found by Sari & Fargione (2019).Specifically, we find that the steady-state TDE rate is about 4 × 10 −6 yr −1 , in good agreement with the order-of-magnitude estimate presented above in Section 2.1.The steady-state EMRI rate settles on about 3.5 × 10 −8 yr −1 , or about 1% of the TDE rate, again in good agreement with the estimate in Section 2.2.
Figure 1 also includes the rates when DCs are taken into account.Shown are the results for both the conservative estimate f R = 1 and the upper limit f R = 2.We see that the TDE rate is almost unaffected by the collisions.This is to be expected, since DCs occur close to the SMBH, while the dominant source of TDEs are stars in the vicinity of R h which plunge into the SMBH after being scattered into the loss cone.These are scattered into an orbit with J ≤ J LC over a time scale of a few P (R h ), so most TDE candidate stars sample the collision-prone region of the cusp only a few times before being disrupted.This insight is consistent with the DC rates presented in the bottom panel in Figure 1.These rates tend to settle on a steady-state value of slightly over 10 −7 yr −1 , in very good agreement with the estimate presented in Equation 10; DCs are mostly due to stars which diffuse in the r ≤ R col region gradually and spend most of their orbits in the inner cusp.The fraction of TDE candidates that were lost to collisions is small.
The effect on the EMRI rate, on the other hand, is dramatic.Collisions essentially deplete the inner part of the  cusp, and since stars can collide prior to any orbital evolution, EMRIs are suppressed almost completely.The EMRI rate drops to few times 10 −9 yr −1 , 10% or less of the estimated rate without collisions.
Note that the steady-state value of R DC is essentially insensitive to the value of f R (as long as it is nonzero).Collisions dominate destruction in the inner cusp, and so simply balance the rate at which stars are supplied to the region R ≤ R col , which is a general property of the cluster.Correspondingly, the density profile in the inner cusp adjusts so that integrating over (n(r)f R ) 2 at every distance r from the SMBH yields the appropriate collision rate.Setting f R = 2 mostly reduces the steadystate densities by about a factor of 2 compared to the f R = 1 case.The enhanced depletion of n(r) for f R = 2 does have a nonnegligible effect on the EMRI rate, since it is proportional to the number of stars in the inner region of the cusp.But this is of secondary importance, given that R EM RI is very small.
Figure 2 shows the fraction of collisions that occurred in the different volumetric bins.In order to asses the steady state distribution of collisions, these fractions are averaged only over the later half (time wise) of the simulation.Quantitatively, we find that about 99% of collisions occur at distances ≳ 3 mpc from the SMBH.This is, of course, consistent with the total DC rate found above, recalling that most stars that are supplied into the r ≤ R col region of the cusp evolve first into an r ≲ R col orbit.The total energy per collision at R col is by construction of order 10 49 erg.The implications on the potential observability of such collisions are dire: the thermal energy in the collisions will be small and the optical flare will be very dim.Some of the ejected gas may eventually accrete onto the SMBH, but given the large initial distance, the accretion rate will be low, leading to a prolonged, low-luminosity event, if at all.This last result changes dramatically when stars from disrupted binaries are included in the simulations, as we demonstrate below.

RESULTS PART II: DESTRUCTIVE RATES INCLUDING STARS FROM DISRUPTED BINARIES
As mentioned above, disrupted binaries add a unique channel for transporting stars to inner regions of the stellar cusp.Generally, a binary with separation distance a b is disrupted when it approaches the SMBH closer than its tidal disruption radius, As is the case for single stars, binaries are far more likely to arrive at their tidal disruption radius after being scattered into an elongated orbit around the SMBH rather than by diffusing inward through a series of quasicircular orbits.As a result, a disrupted binary typically has an original total energy close to zero, and once disrupted the three-body process with the SMBH often results in one star being ejected with positive energy (Bromley et al. 2006;Perets & Šubr 2012), with the other being captured in a tight orbit (Rossi et al. 2014).Indeed, ejection through binary disruption is considered the favored production mechanism for the hyper-velocity-Stars (HVS) which are observed leaving the plane of the Milky Way with velocities exceeding 1000 km s −1 (Brown 2015).The bound companion is left in a highly eccentric orbit around the SMBH.Its periapse is approximately equal to the binary tidal disruption radius, r p ≈ R T b (a b ), and energy conservation requires an sma of approximately For the Milky Way SMBH and Sun-like stars, this implies eccentricities of e ≈ 0.99.
Hereafter we refer to stars that were disrupted from binaries and obtain bound, eccentric orbits as "injected stars".Their contributions to the TDE, EMRI, and DC rates, as well as their impact on the stellar density profile, obviously depend on two parameters: the total rate, η b , at which they are injected, and the distribution of injected stars over the (r, r p ) space.The former depends on the fraction of stars which exist in binaries in the vicinity of R h , while the latter reflects the original distribution function of binary separation distances.Some nontrivial assumptions are required to set the distribution of a b , since it is a combined result of the inherent characteristics of binaries and how they can survive being disrupted by other stars (i.e."ionized", Heggie & Hut 1993).Tighter (smaller a b ) binaries should clearly be more resilient to disruption by multibody forces.According to Hopman (2009), the probability of a binary with a given value of a b to survive disruption rises rapidly the further it lies from the SMBH.For Milky Way parameters and Sun-like stars, binaries with a b ≤ 0.2 au should survive disruption at R h .On the other hand, tighter binaries have smaller tidal disruption radii, and therefore smaller loss cones in angular momentum space and a smaller probability of disruption per binary.Conversely, wider (larger a b ) binaries which can survive somewhat outside of R h also have larger tidal disruption radii, and so they too can be scattered into disruptable eccentric orbits close to the SMBH with comparable probabilities.
A full analysis of the proper distribution of binary separation distances is beyond the scope of this work, and hereafter we limit ourselves to two simplified options.As a nominal model, we assume that the maximal value is a max = 1 au (significantly wider binaries are both unlikely to survive in or close to the cusp, and will not contribute to the density in the r ≤ R col region).We also consider the case a max = 0.1 au, which essentially assumes that only very tight binaries exist at any region that serves as a source for injected stars.In any case, we set the lower limit of a min = 0.01 au which is the minimal separation distance for two sun-like stars.

Including stars from disrupted binaries in the simulations
Our algorithm for including stars from disrupted binaries is as follows.When a star is destroyed, a random number 0 ≤ ξ < 1 is uniformly drawn and compared to a predetermined parameter p b .If ξ > p b the star is replaced by a new star in the outer sma bin (as discussed in 3.4); setting p b = 0 means no injection of binaries.If ξ ≤ p b a star is injected with the appropriate parameters from a disrupted binary: r p = r T b (a) and r = 100r p .Note that this method conserves the total number of stars, but the steady-state rate of injected stars is in turn consequential, rather than a predetermined parameter (our method is thus differs from Sari & Fargione (2019), who injected stars at a constant rate).For a given set of external parameters, we can fine tune the parameter p b in order to generate any specific steady-state rate η b .
The original binary sma a b is drawn according to a log-normal distribution which conforms with the observed population of binaries (Duchêne & Kraus 2013).This distribution is assumed to lie between two values, a min and a max .As mentioned above, we set a min = 0.01 au, and consider two values for a max : 1 au and 0.1 au.

The effects of binaries without collisions
Figure 3 presents the results of the simulations when collisions are ignored, in terms of the steady-state values of R T DE and R EM RI , as a function of the steady-state injection rate, η b .Again, we clarify that η b is itself a result of the simulation, which we fine tune (by trial and error) by setting the parameter p b .
We find that inclusion of injected stars from disrupted binaries increases both the TDE rate and the EMRI rate.The results are consistent with those of Sari & Fargione (2019).Given the small initial periapses and high eccentricities, injected stars have a much larger probability (compared to stars that enter the cusp at R h ) to be scattered into orbits that terminate in destruction at R T (see Bromley et al. (2012)).Numerically, we find that about 70% of the stars which were injected over the second half of the simulation were destroyed, in very good agreement with the analytical estimate of Fragione & Sari (2018).The majority of injected stars are destroyed as TDEs, but a significant minority inspiral to R T through the EMRI channel, either because they are initially injected into a combination of (r, r p ) for which GW emission leads to an EMRI, or because they are quickly scattered into such a combination.Hence, both R T DE and R EM RI are enhanced by an addition of order η b , which for the latter can result in an increase of over an order of magnitude in comparison to the case of no binary disruptions.In particular, setting p b = 0.81 yields a corresponding injection rate of η b ≈ 10 −5 yr −1 , and results in steady state-rates of R T DE ≈ 9 × 10 −6 yr −1 and R EM RI ≈ 3 × 10 −6 yr −1 .
An interesting quantity of the destroyed injected stars is the lifetime they experience in the system.Figure 4 presents the distribution function of the lifetimes of injected stars that were destroyed during the later half of the simulation.We emphasize that these are only injected stars, not stars from the original stellar profile or stars added to the cusp at R h .The total function, as well as the breakdown between the TDE channel and the EMRI channel are displayed.We comment that a few percent of the stars are destroyed within one time step of their injection (typically if injected very close to the SMBH), for which the simulation does not have the appropriate temporal resolution.These stars are not depicted in the figure.
We see that the lifetimes are spread over the entire range of a few times 10 5 to 10 10 yr (essentially the half-length of the simulation).The long-lived stars are generally scattered into more circular orbits and then merge with the preexisting population, for which there is a finite probability of being rescattered toward destruction.On the other hand, there is a maximum likelihood for lifetimes of a few 10 6 yr, both in TDEs and in EM-RIs.This reflects the fact that the surplus of stars over the range of injection radii (compared to a BW profile) shortens the two-body scattering time scale (equation 3).As a result, the injected stars, which inherently have large eccentricities, evolve faster than expected in a BW profile, so that a large fraction reach r p = R T either by two-body scatterings or by GW emission.

Including binaries and collisions
Our results change quantitatively when collisions are taken into account.In Figure 5 we show the calculated TDE, EMRI, and DC rates as a function of time, again over approximately one cluster relaxation time, T h .All simulations have the nominal values of f R = 1 for collisions and a max = 1 au for the distribution of the (r, r p ) combinations for the injected stars.In the different simulations we varied the constant p b , generating a different steady-state injection rate, η b .Similar to Section 4, the system roughly settles on a steady state after a time of a few tenths of T h .
As is to be expected, increasing the injection rate means more stars are supplied to the r ≤ R col region, and so significantly more DCs and EMRIs, but only a minor effect on the TDE rate is seen.In this case as well DCs dominate over EMRIs by more than an order of magnitude.This result is directly depicted in Figure 6, which presents the results of the simulations in terms of the steady-state values of R DC , R T DE and R EM RI as a function of the steady-state injection rate, η b .
We find that the inclusion of stars from disrupted binaries causes a dramatic increase of the collision rate, when compared to the same rate when binaries are not included.For low injection rates the new stars combine with those supplied by the cusp, creating larger densities and a significant growth of R DC .For higher values of η b collisions become the dominant destruction mechanism in the inner cusp; they self-regulate the density profile, and the slope of the R DC (η b ) curve is moderated to a roughly linear one.Nonetheless, for the higher range of η b ≳ 5 × 10 −6 yr −1 , R DC is enhanced by more than an order of magnitude (compared to the value at η b = 0), and becomes comparable to R T DE .The latter is now essentially a constant, limited almost entirely to stars plunging from R h on eccentric orbits.
Note that as in Section 5.2, R EM RI is also enhanced by a factor of ∼ 10 for the higher end of η b values.However, since the evolutionary paths to EMRIs are still significantly suppressed by collisions, their total rate remains small.The enhancement of R EM RI due to injected stars is of order 10 −2 η b , rather than by ∼ η b , the result when collisions are ignored.
Another significant consequence of a high injection rate is on the distribution of distances from the SMBH   where collisions occur.This is demonstrated in Figure 7, which shows the resulting distributions (averaged over the second half of the simulation) for the different injection rates presented in Figure 6.
Higher injection rates naturally lead to a higher fraction of collisions which occur at smaller distances from the SMBH.Note that if at the injection radius the collision time scale is shorter than any evolutionary time scale (relaxation or GW emission), the distribution of collision distances essentially reflects the assumed distribution of the binaries' original separation distances.The log-normal distribution of a b between a min = 0.01au and a max = 1au leads to collisions of injected stars being spread roughly in a log-normal fashion in the range 10 − 1000 au from the SMBH.
This result certainly reflects our simplifying assumption that disruption always occurs exactly so that r p = R T b (a b ).In reality there is a finite probability for disruption over a range of distances from the SMBH, which depends also on the distance of closest approach and orbit orientation (Hills 1988;Bromley et al. 2006;Sari & Rossi 2010;Bromley et al. 2012).A full analysis of this effect is beyond the scope of this work, but we do expect our general results to be valid as long as the typical variability is much smaller than the total range of values for a b .For example, Bromley et al. (2006) demonstrated one case where for a given value of a b , the resulting distribution of energies for the captured star is a Gaussian with a standard deviation of about 40% (20% in the velocity of the ejected star).Such a variability, which is both symmetric and relatively small (compared to the range of binary separation distances), should not significantly affect the spatial distribution of collision radii, especially since stars are always captured in very eccentric orbits.In such orbits the collision time is shorter than the relaxation time (even for stars whose apoapse lies slightly above R col ), and collisions mostly occur close to the periapse.It is noteworthy that allowing for a probabilistic outcome in binary disruption also implies a finite probability for tighter binaries not to be disrupted at all.Some of the tightest binaries will be destroyed as a combined TDE (Mandel & Levin 2015), modifying the balance between TDEs and DCs.
Our results suggest that collisions become dominant in the total collision rate already around η b ≈ 10 −6 yr −1 .We emphasize that for η b ≈ 10 −5 yr −1 collisions at and below R ∼ 100 au make up more than 60% of the total collisions.This radius roughly corresponds to total kinetic energies of 10 51 erg, comparable to supernovae in terms of the driving energy source.Our result concurs with the analytic estimate by Balberg et al. (2013), and implies that for steady-state binary disruption rates of order 10 −5 yr −1 , a Milky Way-like galaxy should exhibit "collisional supernovae" with similar rates as TDES.
Applying the parameter survey (f R = 2 or a max = 0.1 au) does indicate an important sensitivity in terms of the distribution of distances from the SMBH where DCs occur.Figure 8 shows the relative fraction of collisions which occur at different distances from the SMBH for simulations with the two variations, in both While the change to f R = 2 only slightly modifies the distribution of distances (again, it simply rescales the local densities), setting a max = 0.1 au funnels the supply of injected stars to smaller distances from the SMBH, thus enhancing the DC rate at such radii.This implies that while the total DC rate is determined by η b , the energy distribution function associated with collisions becomes highly localized.In particular, about one half of the collisions correspond to energies of order 10 52 erg.This result, if realistic, has significant implications concerning the prospects of observing the early optical flare that will follow a collision.Some caution is called for with regards to this conclusion, since if only the tightest binaries are candidates for generating injected stars, their fraction among the overall population in the cusp may be small, resulting in lower injection rates.
The role of the maximum binary separation distance is also reflected in the lifetime of injected stars.The distribution function of lifetimes for the two cases (f R = 1, a max = 1.0 au) and (f R = 1, a max = 0.1 au) over the second half of the simulation is shown in Figure 9. Since more than 90% of the injected stars are destroyed by collisions in both cases, we only show the total distribution function (as opposed to the breakdown of destruction channels shown in Figure 4).
We find that in both cases the bulk of the lifetime distribution function centers around times of order 10 7 − 10 8 years.Since collisions deplete the density of stars in the inner part of the cusp, the two-body gravitational timescales the injected stars experience are longer than in a BW profile.The overabundance of shorter (∼ 10 5 − 10 6 yr) lifetimes, seen in Figure 4, is therefore absent here, consistent with the suppression of the TDE and EMRI channels for injected stars.We note that setting a max = 0.1 au confines the injected stars to a tight range of radii, and so decreases the spread of lifetimes.Quantitatively we find that the average lifetime of an injected star is about 6 × 10 8 yr for a max = 1 au, and reduces to about 10 8 yr for a max = 0.1 au.Since these steady-state lifetimes correspond to an injection rate η b = 10 −5 yr −1 , these results imply that the momentary number of injected stars in the profile is about 6000 for a max = 1 au, and 1000 for a max = 0.1 au.
To complete the analysis, we present in Figure 10 the steady-state profile of the cusp in terms of N (r) for four representative cases.This profile is calculated by averaging the vector N (i) over the last 10% of the simulation time.The four cases are (i) no disrupted binaries and no collisions, (ii) stars injected due to disrupted binaries with η b = 10 −5 yr −1 but no collisions, (iii) no disrupted binaries and including collisions, and (iv) injected stars at η b = 10 −5 yr −1 and including collisions.We see that without disrupted binaries and without collisions the profile tends to maintain the original BW N (≤ r) ∝ r 5/4 result.The innermost part of the cusp is slightly underpopulated, as is to be expected due to the dilution caused by the accelerated losses through and again with collisions (red and blue symbols, respectively).Shown is N (≤ r), the total number of stars with sma ≤ ri, averaged over the last 10% of the simulations.Also shown (dashed line) is the theoretical profile of Bahcall & Wolf (1976), N (≤ r) = N × (r/R h ) 5/4 , which does not account for disrupted binaries, collisions or relativistic effects on the orbits.Stars are injected with amax = 1 au, and collisions are estimated with fR = 1 in the collision cross section.
GW emission.When the injection rate is significant (but still ignoring DCs), there is an obvious excess at small radii with respect to the BW profile, since stars are injected specifically into the inner cusp, at the expense of the region near R h .These results are naturally consistent with the original analysis by Sari & Fargione (2019).When collisions are taken into account the interior part of the cusp settles into a different steady state, while the exterior (near R h ) steady-state profile is unchanged (and depends only on the mechanisms which supply stars to the outer cusp).Collisions clearly create a significant depletion mechanism for the inner part of the cusp.Without injected stars, the main source of stars for this region is through the diffusion of stars from larger radii, so the steady state is achieved through collisions close to the (fiducial) value of R col ≈ 8 mpc set in the simulation.Since most stars arrive at r ∼ R col along quasi-circular orbits, N (r) drops over a very short range from the BW slope to practically zero (see Section 2.3).When the injection rate is large, it establishes a steady population which balances DCs in the inner cusp and N (r) there is significant.Moreover, recall that injected stars have highly eccentric orbits, so most collisions occur closer to the periapse, rather than at the sma (compare the distribution of collision distances in Figure 8 to the N (r) profile).

SUMMARY AND CONCLUSIONS
Recent and upcoming advances in time-domain astronomy offer new opportunities for studying energetic transient sources.Dedicated surveys, such as ZTF1 (Bellm et al. 2019) and ASAS-SN2 (Kochanek et al. 2017) have dramatically increased detection rates for energetic high-cadence events, and these will be enhanced further in terms of discovery rates and observed wavelength (near UV) by the upcoming ULTRASAT3 mission (Shvartsvald et al. 2023).The prospects of detecting luminous and rare transients in general are on the rise, and these include sources which are unique to galactic nuclei where an SMBH reigns over the dynamics and the fate of the stars surrounding it.
In galactic nuclei harboring an SMBH, the stellar orbits are dominated by the deep potential well of the SMBH, while two-body interactions among the stars modify their energy and angular momentum, and shape their distribution over longer timescales.Stars may be scattered into extremely eccentric orbits, which will occasionally result in a star being tidally disrupted if it ventures too close to the SMBH.Moreover, stars sufficiently close to the SMBH also evolve their orbits through energy loss due to GW radiation emission.If this energy loss is rapid enough, these stars will inspiral individually toward the SMBH on a general relativistic timescale, and disrupt "gradually" as they approach and cross the tidal disruption radius.In this case, a possible electromagnetic flare will be accompanied with a GW event with emission in the millihertz frequency band (Alexander 2017).
In this work, we examined the rates of high-velocity destructive stellar collisions which occur close to the SMBH.Such collisions compete with the other destructive channels.We used a simplified numerical Monte Carlo analysis of stellar dynamics in galactic stellar nuclei.We applied the specific version by Sari & Fargione (2019) for simulating two-body relaxation and GW emission effects on stars orbiting an SMBH, to which we added self-consistent tracking of DCs.Our main results concern the relative rates of the different channels of stellar destruction near the SMBH, which are TDEs, EMRIs, and DCs.
We find that for a gravitationally relaxed Milky Waylike stellar cluster and SMBH, the DC rate is of order a few percent of the TDE rate.The latter is dominated by the scattering of stars into eccentric orbits in the vicin-ity of the radius of influence, and the former is the main destruction mechanism of stars which diffuse into the inner part of the cusp.When compared to previous works which did not account for collisions, we find that DCs replace EMRIs as the secondary destruction channel in the cusp.Once collisions are included, main-sequencestar EMRIs become practically nonexistent (of course EMRIs of compact objects should not be affected).
We also examined the role of binary tidal disruptions in the destructive processes in galactic nuclei.For stars injected on highly eccentric orbits close to the SMBH due to Hills binary disruption, we find that they are almost exclusively destroyed through high-velocity collisions.In other words, the DC rate essentially settles on a steady state with the injection rate.The TDE rate is unaffected (since it is dominated by stars scattered into eccentric orbits far from the SMBH), and the EMRI rate remains very small (although binary disruption does increase this rate due to the finite probability of injected stars to avoid collisions and eventually inspiral toward the SMBH).For a Milky Way-like system with Sun-like stars, we find that injection rates of order 10 5 yr −1 will generate a DC rate which is comparable to the TDE rate.
Again, it is noteworthy that these results significantly modify the expected EMRI rates for high injection rates when collisions are not considered.In such simulations (see also Sari & Fargione (2019)), injected stars have sizable probabilities to either be disrupted by the SMBH on plunging orbits or on gradually circularized orbits due to GW emission, so high injection rates predict comparable EMRI and TDE rates.Allowing for collisions is therefore crucial when considering the prospects of detecting luminous GW emission due to main-sequence-star EM-RIs, or for interacting EMRIs to create quasiperiodic emissions (Metzger et al. 2022).
An important aspect of our results is the spatial distribution of high-velocity stellar collisions.In our simplified model we assumed that collisions are "destructive" when the typical velocity, estimated according to the SMBH potential, is greater than the stellar escape velocity.For Sun-like stars DCs correspond to energies of order ∼ 10 49 erg and higher.In terms of the distance from the SMBH, for Milky-Way parameters this corresponds to a distance of order R col ≤ 8 mpc.In a steady state, gravitationally relaxed system, stars are supplied to the inner part of the cusp through diffusion in energy and angular momentum space; in other words, from the outer cusp to its interior.As a result, almost all collisions occur close to R col , which implies that the optical flare following the collision will be very dim, and the debris will, at best, accrete very slowly onto the SMBH.The observational potential of such collisions is poor.However, stars from disrupted binaries are injected much closer to the SMBH than R col , and their collisions will respectively be more energetic.We reproduce the rough analytic estimate of Balberg et al. (2013), that most DCs of injected stars will be at energies of order 10 51 − 10 52 erg, implying a short but bright optical flare.These are "collisional supernovae", whose rate might be comparable to the TDE rate, if the binary disruption rate is significant.Some simplifications in our model should be highlighted as points for further consideration.The Monte Carlo two-body averaged scheme cannot take into account coherent torques between slowly precessing orbits, namely resonant relaxation processes (Kocsis & Tremaine 2011, 2015).On average resonant relaxation is relevant at radii beyond the DC regime (r ≳ 1000 au), since at smaller distances relativistic precession decouples the GW inspiral from the residual torques of the background stars.This result should be reevaluated when considering a stellar profile which is modified by injected stars.We plan to examine this issue separately with designated N -body simulations (G.Yassur and S. Balberg, 2023, in preparation).
In our model we also ignored the role of low-velocity collisions in the outer parts of the stellar cusp.These are generally not destructive, and more likely result in mass transfer, mass loss, and mergers (Freitag & Benz 2005).Multiple collisions may be however be destructive, and in any case, changes in the number and masses of single stars may accumulate to quantitative effects on the dynamics and steady-state profile in the outer cusp (Sills et al. 2005;Dale & Davies 2006;Rose et al. 2023).In particular, they might enhance the formation of stellar mass black holes (Rose et al. 2022), which will eventually affect the cusp dynamics.
Finally, we clearly oversimplified our analysis by considering a population of single-mass stars.Realistically, there should be a mass function both for the stars in the cusp and the injected stars, which implies a complex relation between the mass and number densities.Moreover, a full mass function will result in segregation due to dynamical friction, and objects with larger masses will generally sink toward the center thus adjusting the cusp density profile.Most notably, there is likely to be segregation of compact objects towards the inner part of the cusp (Keshet et al. 2009;Linial & Sari 2022), with stellar black holes leading with the steepest density profile.Compact objects will affect main-sequence stars dynamically by shortening the two-body relaxation timescales, without participating in collisions (although a compact object may still disrupt a main-sequence star, depending on the impact velocity; see the appendix in Metzger & Stone (2017)).Thus a large population of compact remnants in the inner cusp may change the relative mix of DCs and main-sequence EMRIs.The dynamics and destruction rates can be complicated even further if one or more intermediate-mass (∼ 10 3 M ⊙ ) black hole is present in the inner cusp (Frgaione et al. 2018a,b).We address some these complications with N -body simulations of the inner cusp in our separate work (G.Yassur and S. Balberg, 2023, in preparation).

Figure 1 .
Figure 1.Top to bottom: the TDE, EMRI, and DC rates, respectively, as a function of time with and without allowing for destructive stellar collisions.The collision cross section is taken as π(fRR⋆) 2 ; setting fR = 0 corresponds to no collisions.Notably, including collisions hardly affects the TDE rate, but suppresses EMRIs almost completely.

Figure 2 .
Figure 2. The fraction of collisions as a function of distance from the SMBH out of the total number of collisions, averaged over the later half of the simulations.

Figure 3 .
Figure 3. Steady-state TDE and EMRI rates as a function of the injection rate of stars from disrupted binaries.Collisions are ignored.The rate η b is generated by varying the probability p b , where stars destroyed in a TDE or an EMRI are replaced with an injected star with probability p b , and by a star at R h with probability 1−p b .Shown are the results for the nominal assumption, amax = 1 au and for the tight binary assumption, amax = 0.1 au.

Figure 4 .
Figure 4.The distribution of lifetimes experienced by the injected stars, for simulations with η b = 10 −5 yr −1 , and for the nominal value amax = 1 au.Shown are the total fraction and separately the fractions for stars destroyed as a TDE and stars destroyed as an EMRI.Collisions are ignored.

Figure 5 .
Figure 5. Top to bottom: the TDE, EMRI, and DC rates, respectively, as a function of time, for various values of the injection rate, η b .Different curves correspond to different injection rates, noted in the legend.The collision cross section is taken as π(fRR⋆) 2 , with fR = 1, and the injected stars are drawn from a log-normal distribution with amax = 1 au (see Section 5.1).

Figure 6 .
Figure 6.TDE (black), EMRI (blue), and DC (red) rates as a function of the injection rate of stars from disrupted binaries, η b .Shown are the results for the nominal combination fR = 1 and amax = 1 au, and also for changing (separately) to fR = 2 and amax = 0.1 au.Note that these changes have a noticeable effect only on REMRI .

Figure 7 .
Figure 7.The fraction of collisions as a function of distance from the SMBH out of the total number of collisions, averaged over the later half of the simulations.Different curves correspond to the different injection rates shown in Figure 5, as noted in the legend.

Figure 8 .
Figure 8.The fraction of collisions as a function of distance from the SMBH out of the total number of collisions, averaged over the later half of the simulations, for an injection rate of η b = 10 −5 yr −1 .Shown are the results for the nominal values of fR = 1, amax = 1 au, and also for the variations fR = 2 or fR = 1 amax = 0.1 au.

Figure 9 .
Figure 9.The fraction of lifetimes experienced by the injected stars, for simulations with collisions and η b = 10 −5 yr −1 .For this injection rate over 90% of the injected stars are destroyed through collisions.Shown are the results for the nominal values of fR = 1, while varying the maximum initial binary separation distance by amax = 1 au (black) amax = 0.1 au (red).

Figure 10 .
Figure10.The cusp profile without (η b = 0) and with stars form disrupted binaries (injected at a rate η b = 10 −5 yr −1 ), without collisions (black and green symbols, respectively) and again with collisions (red and blue symbols, respectively).Shown is N (≤ r), the total number of stars with sma ≤ ri, averaged over the last 10% of the simulations.Also shown (dashed line) is the theoretical profile ofBahcall & Wolf (1976), N (≤ r) = N × (r/R h ) 5/4 , which does not account for disrupted binaries, collisions or relativistic effects on the orbits.Stars are injected with amax = 1 au, and collisions are estimated with fR = 1 in the collision cross section.