Segregation and Collisions in Galactic Nuclei: Rates of Destructive Events Near a Supermassive Black Hole

The centers of galaxies host a supermassive black hole (SMBH) surrounded by a dense stellar cluster. The cluster is expected to develop mass segregation, in which gravitational scatterings among the stars cause heavier objects to sink closer to the central black hole, while lighter objects will tend to be overconcentrated in the outer regions. This work focuses on the implications of mass segregation on the different channels for violent destruction of stars in the cluster: tidal disruptions, gravitational-wave-driven inspirals and high-velocity destructive collisions between stars. All such events occur close to the central black hole, where the heavier objects congregate. The analysis is based on a simplified Monte Carlo simulation, which evolves a two-mass population in a cluster surrounding a Milky Way–like SMBH. The simulation is based on the single-mass scheme used by Sari & Fragione and Balberg & Yassur, which has been extended to allow for the dynamical friction effects typical of unequal-mass populations. The effects of mass segregation on the rates of the different destruction channels are analyzed self-consistently in the overall evolution of the cluster. Also considered are stars which are injected into the cluster after being disrupted from a binary system by the SMBH. Such stars are captured in the inner regions of the cluster, and so their orbital evolution, as well as their destruction rate, are therefore influenced by heavy objects that might be abundant in the vicinity of the SMBH.


INTRODUCTION
Supermassive black holes (SMBHs) in the centers of galaxies are typically surrounded by a dense cluster of stars.The distribution and dynamics of the stars are mostly determined by a combination of the SMBH's dominant potential and gravitational interactions between the stars (Alexander 2017).These dynamics give rise to multiple channels for the violent destruction of stars in the cluster, specifically tidal disruption of stars (Rees 1988;Gezari 2021), gravitational-wave emission as objects inspiral into the SMBH (Barack & Cutler 2004;Amaro-Seoane & Gair 2007), and high-velocity stellar Corresponding author: Shmuel Balberg shmuel.balberg@mail.huji.ac.il collisions (Balberg et al. 2013;Amaro-Seoane 2023).While relatively rare, such events are likely to produce unique observational signatures, and reliable estimates of their rates are certainly required.
Several estimates of the rates for tidal disruption events (TDEs) and extreme-mass-ratio inspirals (EM-RIs) have been carried out in the past (see, e.g., Syer & Ulmer (1999); Magorrian & Tremaine (1999); Wang & Merritt (2004); Brockamp et al. (2011); Aharon & Perets (2016); Panamarev et al. (2019); Stone et al. (2020); Broggi et al. (2022)).Less attention has been given to destructive collisions (DCs), which in fact cannot be neglected among the processes which determine the steady-state profile of the cluster.In a recent work, Balberg & Yassur (2023) expanded the approximate Monte Carlo-based Henon (1971) method, developed by Fra-Balberg gione & Sari (2018) and Sari & Fragione (2019) for a single-mass population of stars, by including a model for DCs.In particular, they found that, once accounted for, DCs replace EMRIs as the secondary process which depletes the innermost part of the stellar "cusp."This is especially true when considering stars that are captured near the SMBH following binary tidal disruption -the Hills (1988) mechanism.If this capture rate is significant, DCs between captured stars might occur at a rate comparable to that of TDEs.
The single-mass population approximation mentioned above is generally considered a reasonable representation of a stellar cluster and has indeed been used in many previous studies.Specifically, the general setting can be assessed based on the groundbreaking analysis of Bahcall & Wolf (1976), who developed a fundamental solution for a single-mass population of stars assumed to have an isotropic velocity field.By approximating the gravitational potential of the SMBH as Newtonian, and treating all interactions between stars as random two-body gravitational scatterings of point objects, they found that the steady state number density profile follows a power law of n(r) ∼ r −α with α = 7/4 (hereafter, the BW profile).This result, which applies when the SMBH dominates the gravitational potential, serves as a point of reference for analyses of galactic nuclei.
In reality, of course, stars in the cluster, are distributed in a continuous mass function.It is inevitable that such a stellar cluster will settle on a steady state near an SMBH that includes "mass segregation," in which heavier objects are overconcentrated at smaller distances from the SMBH and lighter objects are overconcentrated further out.This is due to the tendency toward energy equipartition in multiple gravitational encounters, which will cause the heavier objects to sink in the potential well of the SMBH.In essence, this is another manifestation of dynamical friction (Chandrasekhar 1943).The point of reference for mass segregation in galactic nuclei was again set by Bahcall & Wolf (1977), who generalized their original work by including multiple mass groups.By limiting their analysis to the case where the massive stars are also most abundant, they found a simple analytic solution, in which the most massive stars of mass m = m max settle on the single-mass BW profile, while lighter stars with m < m max have a somewhat shallower steady-state profile of n m (r) ∝ r −(3/2+m/4mmax) .
A sufficiently top-heavy stellar mass function is not necessarily realistic.Alexander & Hopman (2009) coined the term "weak segregation" to describe this case when heavier objects dominate their own scatterings all through the cluster.They demonstrated that in the opposite scenario, where heavier objects are rare and their gravitational scatterings are dominated by the lighter objects all the way to the inner boundary, there is "strong segregation" (see also Preto & Amaro-Seoane (2010)).In this case the "cusp" of heavier objects settles on a density profile which is even steeper than the BW solution, with α = 1.75 − 2.25, depending on the details.Results for weak and strong segregation were since investigated in additional studies (Keshet et al. 2009;Merritt 2010;Amaro-Seoane & Preto 2011;Antonini 2012;Aharon & Perets 2016;Baumgardt et al. 2018;Zhang & Amaro-Seoane 2023), and generally reproduced the analysis by Alexander & Hopman (2009) (see, however, the recent critique by Linial & Sari (2022) and below in the results of §4).
Since they occur close to the SMBH, rare destructive events may be particularly sensitive to the existence of a component of heavy objects.While in some cases general relativistic effects and the impact of segregation on the EMRI rate were explicitly considered (Hopman & Alexander 2006;Aharon & Perets 2016;Bar-Or & Alexander 2016), the possibility of stellar collisions in conjunction with mass segregation has not yet been addressed.Most notably, there is the possible impact on stars from disrupted binaries, which are generally "injected" close enough to the SMBH so that they constitute a major component of the stellar destruction rate (Sari & Fragione 2019).Moreover, their collision time scale is significantly shorter than the two-body relaxation time scale, and they dramatically boost the DC rate (Balberg & Yassur 2023).
The goal of the present work is to examine the potential effect of mass segregation on the expected rates of the violent destruction events of stars.Heavier objects that have segregated to the interior regions of the cluster will reduce the two-body relaxation time, and also cause the lighter stars to gain energy on average through dynamical friction.Both effects, if significant, may indeed modify the rates of destructive events.The evaluation of these rates is done by generalizing the code described in Balberg & Yassur (2023) to a two-mass population.The code tracks two-body scatterings, GW losses and collisions self-consistently, and has been extended to include approximate treatment of the transfer of energy and angular momentum between the two components.The analysis includes the case when captured stars from disrupted binaries contribute a significant source to the main-sequence (MS) stellar density close to the SMBH.
This paper is structured as follows.Section 2 reviews the principal equations that describe two-body relaxation and dynamical friction in a two-mass system, and presents the approximations used later in this work.For completeness, the main features of the numerical code, presented in Balberg & Yassur (2023), are summarized in Section 3, noting the extensions introduced to describe dynamical friction.The steady-state density profiles of test simulations for a case with no relativistic effects, no collisions, and no captured stars from disrupted binaries are shown in Section 4, and compared to standard segregation theory as mentioned above.Section 5 presents the main results of this work, when including GW losses, DCs, and possibly a significant source of stars from binary disruptions.The presentation focuses on the corresponding rates of destructive events.Conclusions and a summary are offered in Section 6.

RELAXATION AND DYNAMICAL FRICTION IN A TWO-MASS SYSTEM
The principal features of the stellar cluster surrounding the SMBH can be assessed based on relatively simple models (Binney & Tremaine 2008;Merritt 2013).If two-body gravitational scatterings are completely random, they cause stars in the stellar cluster to diffuse in angular momentum and energy space.The typical time scale for a star orbiting the SMBH on a circular orbit with a semimajor-axis (sma) r to change its specific angular momentum, J c (r), by order unity is (Lightman & Shapiro 1977) . (1) In the first equation n(r) and v(r) are the stellar number density and typical velocity at r, respectively, and ln Λ is the appropriate Coulomb logarithm.The quantity ⟨m(r)⟩ is the average mass of the stars in the vicinity of r.In the second equation it is assumed that the stellar cluster is a two-mass population with masses m L and m H for the lighter and heavier stars, respectively (and using the approximate relation n(r) ∼ (N L (r) + N H (r))/r 3 ).Other aspects of the second equation reflect a simplifying assumption that the SMBH mass, M • , completely dominates the gravitational potential.This should apply inside of the radius of influence, R h , defined as the radial distance at which the enclosed stellar mass is equal to M • .Substituting appropriately that v(r) ∼ (GM • /r) 1/2 , and P (r) ∼ (r 3 /GM • ) 1/2 is the orbital period of the star yields the final result of the right-hand side.In the case of a single-mass population equation 1 reduces to a simple expression with N (r)m 2 in the denominator.Note that N L (r) and N H (r) represent the number stars of each type in the vicinity of r, which can effectively scatter a test star at r.These val-ues roughly correspond to the total number of stars of each type enclosed up to a radius r (as long the number density profile is not n(r) ∼ r −4 or steeper, which is a reasonable assumption in all that follows).The energy relaxation time is similar (slightly longer), and changes modestly even for more eccentric orbits with the same sma r.However, the angular momentum relaxation time depends on the eccentricity of the orbit.The diffusive (quadratic) time evolution of angular momentum with time implies that the time scale for an order-of-unity change in the specific angular momentum of an eccentric orbit is smaller by a factor of (J(r, r p )/J c (r)) 2 = r p /r, where r p is the periapse distance.The relevant two-body time scale eccentric orbits is therefore (Binney & Tremaine 2008), The relaxation times in equations 1 and 2 are independent of the mass of the test star.This simply reflects the fact that gravitational acceleration depends only on the scattering stars, which are represented by the second moment of the mass distribution.Equation 2 thus provides a first indication regarding the possible influence of segregation on the DC rate.When stars captured in eccentric orbits following binary disruptions are the dominating source for collisions, both the rate and the spatial distribution are influenced by the fact that these stars are captured on very eccentric orbits.Specifically, Balberg & Yassur (2023) found that for the Milky Way SMBH, if stars are captured near the SMBH at a rate of 10 −5 yr −1 , equilibrium created by collisions is characterized by a steady-state number of a few thousand captured stars.Given the quadratic dependence on mass in the relaxation time, even a small number of heavy objects which have segregated to the inner part of the cluster can shorten the relaxation time, and possibly modify the DC rate and distribution.The two-body relaxation time is derived by assuming that the effect of gravitational scatterings is limited to a random-walk diffusion in angular momentum space.Any specific test star will also experience a net drift in angular momentum and kinetic energy due to dynamical friction created by the other stars.In general, the typical rates of drift in the velocity in the direction of motion, v ∥ , and the specific kinetic energy, E, of a single test star of mass m t in a field of stars with mass m f are (Binney & Tremaine 2008;Merritt 2013) (numerical prefactors of 16π 2 have been dropped).In equations 3 and 4, v t is the velocity of the test star, and f (v f ) is the velocity distribution function of the field stars, assumed to be uniform and isotropic.Unlike random scatterings, these dynamical friction formulae are linear and also depend on the mass of the test star, as it creates a wake in the field stars and friction is the result of the back reaction on the test star itself.The sign in equation 4 conforms with the standard of setting the energy of the bound stars as positive.
Equations 3 and 4 explicitly reflect the standard approximation that the net specific contribution to dynamical friction of a field star depends on whether its velocity is lower or higher than v t .Specifically, in equation 4 it is assumed that field stars with velocity v f < v t combine to affect a cooling term, while stars with v f > v t induce a heating term.Henceforth, it is assumed that at any radial distance from the SMBH r x , all stars have a velocity distribution function which is tight, symmetric, and centered around v(r x ) = (GM/r x ) 1/2 .The rates of change of the parallel component of the velocity and the specific energy in a circular orbit are correspondingly simplified to and where n f (r) is the number density of field stars at r.The factor 1/2 reflects the assumption that approximately half of the field stars have velocities v f < v t .
In general, dynamical friction will occur in any combinations of masses, including a single-mass case (m t = m f ).However, given that the focus of this work is on mass segregation which arises in a two-mass system (m t ̸ = m f ), only net transfer of angular momentum and energy between the two populations is considered.It is also naturally accessible within the approximations used in the Monte Carlo code described below (as opposed to dynamical friction within a single-mass population, which requires an actual velocity distribution).This constraint arises directly for the rate of change of the energy (equation 6 which is identically zero for m t = m f ), but the approximate form for transfer of angular momentum (equation 5) is also only applicable for the two-mass case.Correspondingly, the signs of changes in equations 5 and 6 are set so the top (bottom) option applies to the lighter (heavier) star.This choice leads to dynamical friction creating the required bias in the diffusive nature of gravitational scatterings.It effects a deterministic gain in angular momentum and loss in energy for the lighter stars, and conversely for the heavier stars, again -maintaining the convention that bound orbits have positive energies.
The approximate forms of equations 5 and 6 obviously do not contain a full description of dynamical friction processes.Specifically, Antonini & Merritt (2012) demonstrated that this approximation is not truly appropriate for stars orbiting an SMBH (as opposed to a Maxwellian distribution, for which Chandrasekhar (1943) introduced the approximation).In fact, even light stars with v f > v t will on average generate some dynamical friction on the heavier star's velocity.Nonetheless, given the level of the approximation used in the present work, equations 5 and 6 will be used as an approximate description of dynamical friction among the stars, which ultimately leads to mass segregation.
Dynamical friction sets another gravitational time scale in the system.The time scale for a change of order unity in the angular momentum of a star on a roughly circular orbit through the drift effect can be estimated from equation 5 to be Given that the rate of change of specific angular momentum is J −1 dJ/t = v −1 dv ∥ /dt regardless of eccentricity, equation 7 actually applies to all orbits, independent of the ratio r p /r.Thus, equation 7 shows the general dynamical friction time, T J DF (r, r p ; L/H) of a star with sma r and periapse r p .Note the distinction between the lighter (L) and heavier (H) stars, explicitly introducing a dependence on the mass of the scattered star.Generally, we will find that T DF (r, r p ; L) ̸ = T DF (r, r p ; H).Moreover, the ratio T J DF (r, r p ; L/H)/T J 2B (r, r p ) indicates which of the two processes is more efficient in terms of driving the evolution of a particular component of the population for a given combination of (r, r p ).For the heavier stars, The dynamical friction time scale can be either longer or shorter than the diffusive two-body relaxation time scale, depending on the ratio m L /m H and on the local ratio of N L (r)/N H (r). In fact, the dominant mechanism (swing or drift) may be different for the lighter and heavier stars.It is noteworthy that drift is most relevant for the heavier stars in the strong segregation limit where N H (r)m 2 H << N L (r)m 2 L , so that the dynamical friction time scale for the heavier stars on circular orbits is shorter by a factor of m L /m H .In the opposite weak segregation (Bahcall & Wolf 1977) case in which the heavier stars dominate, N H (r)m H > N L (r)m L , the dynamical friction time scale may also be the shortest for lighter stars, but only by a moderate factor of (1 + m L /m H ) −1 .

MONTE CARLO SIMULATION WITH TWO MASSES AND DYNAMICAL FRICTION
Stellar clusters surrounding an SMBH are often analyzed by Fokker-Planck solutions of diffusion in energy and angular momentum, specifically in a reduced onedimensional formulation by assuming an isotropic velocity field (Binney & Tremaine 2008;Merritt 2013).However, studying specific complex phenomena related to the stellar cluster requires more comprehensive solutions for the properties of the cluster.Exact solutions are currently unavailable: First and foremost, because of the need to track a very large number of objects over an extreme range of time scales.Correspondingly, many approximate methods have been suggested, which vary in content and technique depending on the particular applications.
The cluster evolution numerical code used in this work is an extension of the single-mass code presented in Balberg & Yassur (2023).In the following, the general key features and approximations used are summarized, and the main focus of the section is on the extension to a two-mass population of stars: The implementation in two-body random scatterings and the introduction of dynamical friction ( §3.3).Other features of the code are described briefly in §3.4 and the reader is referred to Balberg & Yassur (2023) for further detail.

Monte Carlo Code Setup
The simulations were carried out with a simplified Monte Carlo algorithm.It is based on the method originally suggested by Henon (1971), and applies the technical approach presented by Fragione & Sari (2018) and Sari & Fragione (2019).The code evolves a cluster of stars surrounding an SMBH, including an effective treatment of random two-body scatterings, GW losses and a procedure for allowing and tracking destructive highvelocity collisions.
Each star is tracked individually by its specific energy, E, and specific angular momentum, J.These are translated to an sma r, periapse r p , eccentricity e and period P , assuming an instantaneous Keplerian orbit in the potential of the SMBH.The specific energy is therefore E = GM • /2r.The cumulative effect of gravitational scatterings off other stars is averaged over the instantaneous stellar profile, assuming spherical symmetry.
While 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.The intervals are set as 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.The structure of the cluster at each time step is then also represented by the number of stars of each type in every particular bin, N L (i, j) and N H (i, j).The time-dependent profile of the cluster is described by the instantaneous total number of stars with an sma enclosed in the radial bin i, regardless of periapse,

Time step
It is computationally prohibitive to evolve all of the stars in the cluster in every time step.Grouping the stars according to bins in the (r, r p ) plane allows to selectively evolve only a subset of the stars -specifically, those for which a sufficient time has passed since their last evolutionary update (Fragione & Sari 2018;Sari & Fragione 2019).In order to accomplish this, at the beginning of each time step the simulation scans the characteristic time scales for every bin in the (r, r p ) plane.These are the two-body relaxation time scale, and the light-mass dynamical friction and heavy-mass dynamical friction time scales, respectively: where P (R h ) is the orbital period at the radius of influence, γ = 1.5 is a control parameter used to account for the finite bin size (Fragione & Sari 2018), and f T = 0.1 is the fraction of the relevant time scale allowed in a single evolutionary time step.The Coulumb logarithm is set as ln Λ = 10 (Merritt 2013).The Balberg above equations are derived using the approximate relations i (for each species separately).At the beginning of each time step the code finds the minimum of the three time scales for each (r, r p ) bin, ∆t min (i, j), and then sets the global time step, ∆t, according to the minimum over the entire grid: ∆t = min i,j t min (i, j) . (13) Finally, the evolution due to gravitational scatterings is tracked with the scheme suggested by Sari & Fragione (2019).For each grid cell (i, j) the simulations keep the last time it was updated, t old (i, j).Stars in the grid cell are updated with gravitational scattering effects only when the current time, t, satisfies The time steps are typically a few tenths of P (R h ), while the cluster relaxation time, T 2B , is of order 10 5 P (R h ).
It is this scheme that allows to evolve what is a usually a small fraction of the stars in each time step, and still track ∼ 10 6 stars over the entire length of the simulation.Equations 9-11 hold an implicit approximation regarding the number of stars which effectively scatter a given test star at its radial position r ⋆ .The equations substitute this number with the number of stars with an sma that corresponds to the i-th radial bin which includes r ⋆ .This is a reasonable assumption for large N (r ⋆ ), which is generally the case in the outer part of the simulation zone.On the other hand, close to the SMBH the numbers of stars are smaller, but the relative effect of two-body scattering is also small to begin with.Hence, this approximation is viable and is applied in the calculations reported below.For further validation, several calculations were tested with a different approximation in which the code calculates in every time step the actual effective number density of stars, n(r), in each radial shell by finding the fraction of the orbit each star spends in every radial shell (easily found for Keplerian orbits around the SMBH, as is done for tracking DCs : see equations 28-30 in Balberg & Yassur (2023)).The results were generally similar, but this method is inferior since it includes an inevitable inconsistency.At some radius the scheme for counting stars must be changed back from n(r i )r 3 i to N (r i ), since in the region r > 0.1R h , where most of the stars reside, some stars spend parts of their orbits beyond R h .In reality, other stars penetrate from r > R h to complete n(r), but these are not included in the simulation.The qualitative and quantitative results were found to be very similar, and therefore the approximation that the effec-tive number of scattering stars at r ⋆ is N (r ⋆ ) is used in all simulations reported here.

Evolving the Specific Energy and Angular Momentum by Gravitational Scatterings
In every time step the (r, r p ) bins which need to be updated are determined according to equation 14.The specific energy and angular momentum of stars which correspond to these bins are then evolved as follows.Each star is processed with a time ∆t = ∆t min (i * , j * ), where i * and j * are the indices which correspond to its location in the (r, r p ) grid.The inspected star travels along its orbit through a range of distances from the SMBH lying between the periapse r p * and the apoapse r a * .These correspond to a range of [W 1 , W 2 ] values in the radial grid.During a time step ∆t the inspected star crosses through the W -th bin ∆t/P * times, for a total of N (r W )∆t/P * interactions.
Following Sari & Fragione (2019), the cumulative effect of random scatterings by stars in a given bin is estimated as follows.The typical change in the specific energy of the inspected star by a single field star with mass m f in the W -th bin is of order ∆E 1 ≈ Gm f /r W .Given the random-walk nature of the diffusion in energy space, the cumulative effect of random scatterings by all N L (r W ) + N H (r W ) stars in the W -th radial bin during a time step ∆t min (i * , j * ) is therefore a characteristic change in specific energy of .
(15) Here P * is the orbital period of the inspected star (typically shorter than the time interval ∆t min (i * , j * ), except for stars with R ⋆ ≈ R h ).The square root reflects the diffusive (random-walk) nature of repeated scatterings.
The characteristic change in angular momentum due to random scatterings with the stars in the W −th bin during the given time step is where v(W ) = GM • /r W is the characteristic velocity at r W .The effects of dynamical friction must be calculated separately for the lighter and heavier stars.The typical changes in the specific energy due to dynamical friction between the two populations are calculated per bin W by applying equation ( 6) above and substituting for each species of stars in the partial number densities, and maintaining consistency with the approximate forms of the typical changes in energy and an-gular momentum from random scatterings at r W (equations 15 and 16).Approximating that for each species n and for the changes in energy of the lighter and heavier stars, respectively.Similarly, by applying equation 5, the changes in specific angular momentum through dynamical friction are calculated as and These choices conserve total energy and angular momentum in terms of the transfer between the two populations.Note that the changes induced by dynamical friction are linear in time, representing the deterministic drift caused by this effect.Finally, the changes in specific energy and angular momentum of the inspected star due its passage in the Wth bin during the ∆t(i * , j * ) time interval are updated by adding the effect of random two-body scatterings (Fragione & Sari 2018; Sari & Fragione 2019; Balberg & Yassur 2023) and the deterministic drift of dynamical friction.Denoting the values at the beginning of the time step as E old and J old , the new specific energy and angular momentum are with 0 ≤ χ < 2π and 0 ≤ Φ < 2π drawn randomly from a uniform distribution.The star's orbital properties, r * ,P * , r p * , and e * , are then updated consistently.The process is repeated for the entire range [W 1 , W 2 ] relevant for the inspected star, and for all stars in the (i, j) bins that are being advanced in the global time step.

Additional Features of the Monte Carlo Code
Other physical processes are tracked in the Monte Carlo code by the methods described in Balberg & Yassur (2023).
1. GW losses due to the gravitational field of the SMBH are calculated separately for each star in each time step with standard formulae (Peters 1964;Hopman & Alexander 2005;Sari & Fragione 2019).
2. A MS star, is considered lost when its periapse drops below its tidal disruption radius, R T , which depends on the mass of the star m ⋆ and its radius, R ⋆ .The standard value is (Stone et al. 2020) If the star drops below R T following a gravitational scattering step, the star is considered as having undergone a TDE, while if this occurs after a GW loss step, the star is considered as having been disrupted gradually in a MS-EMRI.
3. A compact object is considered lost when its periapse drops below 4R S (Merritt 2013), where is the Schwarzschild radius (and c is the speed of light).In principle, compact objects are also lost in two types of evolutionary scenarios, either being scattered into eccentric orbits with r p < 4R S or by gradually losing energy in GW emission over multiple quasi-circular orbits.The former, often referred to as "direct plunges" (Hopman & Alexander 2005; Bar-Or & Alexander 2016), are the equivalent of TDEs for the MS stars, while the latter produce "clean" (no disruption) general relativistic EMRIs.While it would be interesting to quantify the details of the evolutionary tracks of individual compact objects (and estimate the resulting GW signal), accurate modeling of the strong relativistic effects close to R S are beyond the scope of this work, and will be addressed separately.
4. DCs are assumed to occur only between MS stars and only below a radius R col , at which the typical gravitational energy in the orbit, GM • m ⋆ /R col , equals the gravitational binding energy of the star, Gm 2 ⋆ /R ⋆ .For the Milky Way SMBH and Sun-like stars, R col ≈ 4 × 10 −3 R h .The probability a star Balberg that has a periapse r p ≤ R col will collide with another MS star during the time step is estimated by finding the instantaneous optical depth for the star along its orbit.The optical depth is calculated per radial bin, i, using the weighted probability of all other stars to be found in the i-th bin (derived from each star's orbital parameters).This allows to assess the overall number density of other MS stars in the bin, n L (r i ).The number density is coupled to a cross section of π(f R R ⋆ ) 2 , with f R being a free parameter, taken here as f R = 1.The probability for a collision of a particular star during the current time step is then compared to a randomly drawn number, and if the star is determined to have undergone a DC, the details are completed in two additional stages.First, the radial bin of the collision is determined by drawing a second random number and comparing it to the partial probabilities for a collision in each bin.Given the radial bin where the collision occurred, the partner MS star is found by comparing a third random number to the weighted probabilities of finding each of the other stars in that radial bin. 5. Stars which were lost due to a TDE, an EMRI, or a DC are replaced with new stars, maintaining a constant number of stars in the simulation.The new star is placed in the simulation through one of two options: (i) representing a star that was scattered into the cluster from outside of the SMBH radius of influence and thus placed in the outer radial bin (between R h /2 and R h ) or (ii) representing a star that was captured close to the SMBH following binary tidal disruption.The relative probabilities for each option are a control parameter of the simulation, and can be adjusted in order to produce a desirable value for the injection rate of stars from disrupted binaries.The properties of injected stars are discussed below in Section 5.2.

RESULTS FOR FUNDAMENTAL CASES OF MASS SEGREGATION
All the simulations described below were carried out assuming a Milky Way-like system with an SMBH of mass M • = 4 × 10 6 M ⊙ , surrounded by a cluster of stars which extends up to a radius of influence of where σ is the measured velocity dispersion external to the radius of influence (Merritt & Ferrarese 2001).The two-mass population is divided between stars with mass m L = 1 M ⊙ and mass m H = 10 M ⊙ , set up so that the total mass is m = M • .The lighter stars are assumed to be MS, Sun-like stars with a radius R ⋆ = R ⊙ , while the heavier stars are essentially stellar-mass black holes and approximated as point objects.The MS stars occupy the range between the tidal radius R T ≈ 1.1 × 10 13 cm, and R h , while the stellar black holes can sink as low as 4R S .This setup of a twomass population serves as a reasonable approximation for studying the general physics of mass segregation of a multi-mass system (Alexander & Hopman 2009).
Simulations are initiated by predetermining the total abundances of each species, N L and N H .At the beginning of the simulation, all stars are assigned values of their sma, r, and specific angular momentum, J.The sma is drawn randomly from a weighted distribution that produces a power-law profile n(r) ∝ r −α , or N (r) ∼ r 3−α .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).This choice generates a thermal distribution, equivalent to an isotropic velocity field.The values of r and J are then used to find the specific energy, E, period, P , eccentricity, e, and periapse, r p , for each star assuming Keplerian orbits around the SMBH.The initial profile is set as α L = α H = 7/4 for both species of stars.This choice generally allows for convergence over time scales of a few tenths of T 2B (R h ) (although in cases with N H ≤ 4 × 10 3 convergence is faster when starting with α H > 7/4, as is predicted by the theoretical results for strong segregation; Alexander & Hopman (2009); Keshet et al. (2009)) .
Figures 1 and 2 present the results of two simplified simulations which serve as tests of the code, as well as a basic demonstration of mass segregation.The values of N (r) are averaged over the last 10% of the simulation time, which is one cluster relaxation time, T 2B (R h ).In these simulations (and all that follow) the steady-state is typically reached after several tenths of T 2B (R h ), so the assumption of steady state appears to be robust.The simulations were carried out without relativistic corrections (the GW loss terms) and ignoring the possibility of collisions.Stars from disrupted binaries are also not taken into account.The figures show, respectively, the steady state stellar profiles found for the two combinations N L = 3.5 × 10 6 ; N H = 5 × 10 4 and N L = 3.96 × 10 6 ; N H = 4 × 10 3 .The profiles are displayed in terms of the number of stars in the i-th radial shell, with r i as the outer radius of that shell (N (r i ) is the number of stars between r i /2 and r i ).
In both cases the code reproduces the general trends expected from theory, and also the standard results  and N -body simulations with a smaller numbers of stars (Freitag et al. 2006;Amaro-Seoane & Preto 2011;Antonini 2012;Baumgardt et al. 2018;Panamarev et al. 2019).For N H = 5 × 10 4 the stellar-mass black holes dominate the two-body scatterings over almost the entire range of the simulation, constituting a weak segregation scenario.Interior to about 0.2R h the steady-state radial profile of the MS stars is numerically consistent with N L (r) ∝ r 1.475 , which is the theoretical estimate of α L = 3/2 + m L /4m H (Bahcall & Wolf 1977) for weak segregation when the heavier objects dominate.It is only in the outer 0.2R h < r region that the lighter stars dominate their own scatterings and their profile turns over to the value typical of strong segregation, with α L ≈ 7/4.These trends coincide with the results for the black holes, which interior to about 0.05R h develop a profile that is consistent with the standard BW power law of N H (r) ∝ r 5/4 (α H = 7/4), while beyond this radius the density profile becomes steeper, as is appropriate for a transition to the strong segregation regime.
The general trends in the case N H = 4 × 10 3 are similar, except that the transition between strong and weak segregation occurs around 0.005 − 0.01R h .Again, in the outer part of the cluster the MS stars dominate the scatterings, and so follow a radial profile that is consistent with the BW value of α L ≈ 1.75, while the black holes follow a steeper profile, roughly consistent with α H = 2.The inner part of the cluster is again consistent with the weak segregation theoretical values of α L = 1.475 and α H = 1.75.
While this transition from weak to strong segregation at increasing radii is not the main topic of the current work, it is noteworthy that that it is an inevitable consequence of segregation.As pointed out recently by Linial & Sari (2022), even for N H << N L , it is always expected that at sufficiently small radii the heavier stars become abundant enough to dominate gravitational scattering, so that in the inner region of the cusp segregation is weak, rather than strong.
This transition from strong to weak segregation in any given profile leads to some further insight.Following Alexander & Hopman (2009), it is customary to evaluate the regime of segregation through the parameter ∆ AH Values of ∆ AH > 1 correspond to weak segregation, while ∆ AH < 1 implies strong segregation.
Most studies use the global values of N L and N H in the cluster to describe the entire system with a single value of ∆ AH .Numerical results are then verified by extracting some characteristic value of the powers α L and α H (typically at some fiducial radius or by averag-Balberg ing over some finite range of radii).However, by nature segregation causes the value of ∆ AH to change across the cluster, increasing inward.This is demonstrated in the lower panels in figures 1 and 2, which show the local values of ∆ AH (r) as a function of radius.Note that in both figures the transition between strong and weak segregation occurs where ∆ AH has a value of a few tenths, as is to be expected based on the relation between ∆ AH and the relevant gravitational times scales (Alexander & Hopman 2009;Preto & Amaro-Seoane 2010).
Finally, another noteworthy result is that in both simulations N H (r) tends to turn over and actually decrease in the very outer radial shell (between R h /2 and R h ), even though this shell is the most voluminous.This result is qualitatively consistent with another prediction by Linial & Sari (2022), that heavier stars should be exponentially rare at large radii.However, since the Monte Carlo code described here focuses on occurrences near the SMBH, it is ill-equipped in terms of resolution to quantitatively reproduce this exponential decline.

TDES, DCS AND EMRIS IN A SEGREGATED CLUSTER
We may now turn to consider the main topic of the current work, which is the possible impact of segregation on destructive events near an SMBH.By selectively adding GW losses and high-velocity stellar collisions, the rates of these events can be estimated as a function of the stellar composition of the cluster, namely as a function of the combination {N L , N H }.

Stellar Profiles
Further simulations with N H = 5 × 10 4 and N H = 4 × 10 3 , were carried out while including (i) GW losses, and (ii) GW losses and destruction of stars in DCs.Figures 3 and 4 show the steady-state radial profiles of the two combinations of stars in the cluster, respectively.Again shown are the profiles found after about one cluster relaxation time, T 2B (R h ), well after a steady state has been reached.The results for the test case of no GW losses, no collisions, and no injected stars (figures 1 and 2) are also shown (black crosses).The top panel in both figures shows the profiles of the MS stars, and the bottom panels that of the black holes.
The inclusion of GW losses alone (no collisions) has a minimal effect on the radial profile.As is to be expected, GW losses play substantial role only close to the SMBH, where the general relativistic timescale is shorter than the gravitational timescales.This causes a minor depletion of stars of both types close the SMBH, so that the profiles deviate from the standard power laws expected from gravitational scatterings alone.As a result, the profiles are very similar to the test cases shown in §4.In fact, even those results also had a slight depletion of stars at inner radii due to the finite inner radius stars can sample (R T or 4R S ), which causes a deviation from a truly isotropic velocity distribution.
Allowing for DCs induces a dramatic effect on the profile of MS stars in the region interior to R col .The collisional time scale for a MS star with orbital parameters {r, r p } can be estimated by (Sari & Fragione 2019; Balberg & Yassur 2023) which for r < R col is shorter than all gravitational time scales in the steady-state profiles derived above.Consequently, regardless of segregation, the region interior to R col becomes practically devoid of MS stars.Scattering from r > R col is too slow and cannot resupply new stars at a sufficient rate to reproduce a profile typical of the collision-free case.Note, however, that there exists a quantitative difference between the profiles of the two cases N H = 4 × 10 3 and N H = 5 × 10 4 .The former is practically devoid of MS stars, and so is very similar to the single-mass (no black holes) profile.On the other hand, a sufficiently large number of black holes drives the MS stars to exhibit a shallower slope, implying that they spread further inward in the cluster.The qualitative reason is that the abundant heavier black holes decrease the relaxation timescale (equation 1) of the MS stars at r ≳ R col .In spite of this diffusion being biased (the net drift of dynamical friction drives the MS stars outward), the enhanced diffusion coefficient allows some MS stars to diffuse inward prior to experiencing a collision.A large number of heavier stars thus leads to a finite steadystate penetration of MS stars to smaller radii.
Finally, note that allowing for collisions has practically no effect on the distribution of the black holes.This is to be expected, since at radii interior to R col the black holes dominate the scatterings in the weak segregation regime.Further depletion of the MS stars through collisions is inconsequential.

The Effects of Disrupted Binaries
Injection of MS stars represents the result of binary disruptions near the SMBH (it is assumed that only MS-MS binaries are relevant for disruptions).The procedure for injecting stars is as follows.For every MS star that is destroyed in any of the channels (TDE, EMRI, or DC), a random number, ξ is drawn and compared to a predetermined probability 0 ≤ p b < 1.If ξ > p b a new star is placed in the outer radial shell (R h /2 < r < R h ) with its orbital parameters drawn randomly assuming a BW profile and an isotropic velocity distribution.For ξ < p b the star is replaced by an injected star, whose properties are drawn based on additional assumptions regarding the original binary population.In the results shown here it is assumed that the binary separation distance, a b , follows 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 = 0.01 au, and a max =1 au.The choice of a max is somewhat arbitrary (although it has to be a relatively small value, since only tight binaries can survive in the dense stellar cluster).
Further detail and sensitivity checks can be found in Balberg & Yassur (2023).
Once the original value of a b has been determined, the injected star is assumed to be captured in an orbit with a periapse r p (a b ) and an sma r(a b ) of where R T (a b ) is the distance from the SMBH where binary tidal disruption is expected.The injected orbit is thus very eccentric, with e ≈ 0.99.These estimates are consistent with the inferred physics of binary disruption (see, e.g.Bromley et al. (2006); Perets & Šubr (2012); Rossi et al. (2014)), in which the three-body process with the SMBH often results in the partner being ejected with a large energy.This inference is supported by the existence of hypervelocity stars (HVS) which are observed leaving the plane of the Milky Way with velocities exceeding 1000 km s −1 (Brown 2015); binary disruption is the favored production mechanism of these stars.
Figures 3 and 4 also show the radial profiles of the stars when the parameter p b is adjusted in each case to produce an assumed injection rate of η b = 10 −5 yr −1 .GW losses and collisions are also included.
As is to be expected, binary disruption creates a significant source of MS stars for the inner region of the cluster.Correspondingly, the abundance of MS stars is obviously greater than in the corresponding case with no injected stars (η b = 0).Since the shortest distances dominate the collision rate, the highly eccentric orbits of the captured stars lead to much shorter collision time scales (equation 27) than stars which diffuse inward at R col .As a result, most injected stars undergo collisions, as is the case for a single-mass population (Balberg & Yassur 2023).However, some minor differences do arise between the two cases, N H = 4 × 10 3 and N H = 5 × 10 4 .While the collisions dominate, the two-body relaxation and dynamical friction time scales are not infinitely longer, so to some extent this effect modifies the distribution and orbits of the injected stars.As discussed above, this effect is more pronounced for a larger black hole population, causing the injected stars to sustain a shallower density profile in the inner part of the cluster.These effects carry on into the rates of destructive events that occur near the SMBH, which is the topic of the following subsection.

Destructive Events near the SMBH and the Effects of Segregation
All three destructive channels of MS stars -TDEs, EMRIs, and DCs -affect the steady-state profile of the stellar cluster.Since all three processes occur close to the SMBH, there is some codependence between them and gravitational scatterings, and specifically consequences of mass segregation.
Figure 5 shows the results of the simulations in terms of the TDE rates in the steady state, as a function of the total number of stellar-mass black holes, N H (when in all cases m = M • ).The curves in the figure correspond to the physical scenarios considered above: no collisions, allowing for collisions, and allowing for collisions when stars are injected (representing disrupted binaries) at a rate of η b = 10 −5 yr −1 .GW losses are included in all simulations.
Recall that TDEs are the result of stars being gravitationally scattered into eccentric orbits which end in them plunging to and below R T .In a gravitationally relaxed cluster, this rate is dominated by stars with an sma close to R h , which are the vast majority of the MS stars.The TDE rate can therefore be roughly approximated by Sari & Fragione (2019) where J C (R h ) is the angular momentum of a circular orbit at R h , and J LC = √ 2GM • R T is the "loss cone" value in angular momentum space, which corresponds to an eccentric orbit with r p = R T .For the Milky Way SMBH and a cluster of Sun-like stars, this value is R T DE ≈ 7 × 10 −6 yr −1 , which is consistent with the result here, particularly for N H = 0. Consequently, mass segregation has a relatively minor effect on the TDE rate, since the black holes are subdominant close to R h .This is indeed reflected in the fact that the TDE rate is almost independent of the existence and number of black holes.The very minor dependence that does exist can be assessed by the rough approximation in equation 30: As N H is increased in the simulations, N L decreases linearly (since the total mass in the cluster is kept fixed), but the two-body relaxation time scale becomes shorter with a quadratic dependence on N H .The latter effect is weak because of the exponential-like decline in the number of heavier stars close to R h , but it still dominates over the former, linear scaling of N L .The TDE rate thus slightly increases for larger values of N H .The opposite is true only for very large values of N H , when N L ≤ 2 × 10 6 .
The effect of mass segregation is far more pronounced when stars from disrupted binaries are injected at the relatively high rate of η b = 10 −5 yr −1 .In the case of a single-mass population (N H = 0), the relaxation time experienced by the injected stars is significantly longer than the collision time and so practically all injected stars are destroyed in DCs (Balberg & Yassur 2023).Hence, they hardly contribute to the TDE rate, which remains unchanged (essentially driven only by stars at R h , as mentioned above).However, in a two-mass system the segregated heavier stars are overconcentrated in the inner part of the cluster (where the stars are injected), and the two-body relaxation time there is clearly shortened, and some injected stars evolve their orbits.Although on average the injected stars gain angular momentum (due to dynamical friction), the stochastic nature of two-body scatterings allows some stars to lose angular momentum so that their periapse drops below R T and they are destroyed as TDEs.The net result is a sizable increase in the total TDE rate, and for large values of N H it reaches a factor of 2 with respect to the single-mass case.
Figure 6 shows the DC rates in the steady state in the simulation as a function of N H .The results are for either collisions but no injected stars, or collisions and stars injected at a rate of η b = 10 −5 yr −1 .
In the absence of segregation (no stellar-mass black holes) the collision time of the MS stars at R col is much shorter than the two-body relaxation time.Since the region interior to R col is essentially depleted (figures 3 and 4), the point of reference in terms of the steady-state DC rate in a cluster without injected stars is approximately where v(R col ) ∼ GM • /R col is the typical velocity at R col .For the Milky Way parameters and N L (N H = 0) = 4 × 10 6 , this corresponds to R DC ∼ 10 −7 yr −1 , which is indeed the result found in the simulations.As already demonstrated above, when heavier black holes are included in the cluster their primary effect is to redistribute the lighter stars.Since in the steady-state profile the MS stars tend to larger radii, their density at R col is smaller than in a single-mass profile, and the total DC rate decreases appropriately.This effect is marginal for N H ≲ 5 × 10 3 , but becomes significant at higher values.Specifically, note that the suppression of the DC rate exceeds the simple effect of N L being smaller than 4 × 10 6 in the single-mass case.For example, for N H = 2.5 × 10 4 , the collision rate has dropped to about 0.6R DC (N H = 0), whereas the total number of MS stars in the simulation is reduced by a factor of only ∼ 0.875 (and N 2 L is smaller by a factor of ∼ 0.765).For very high values of N H (which are probably not realistic), DCs are almost nonexistent.
The DC statistics are substantially different when stars from disrupted binaries are taken into account.Firstly, the reference DC rate is high, with R DC (N H = 0) ≲ η b /2 = 5 × 10 −6 yr −1 .Again, this is simply because the injected stars are almost exclusively destroyed by collisions, with one collision for every two injected stars.Secondly, while the impact of segregation on the total collision rate is also quite significant (for high values of N H ), the origin of this effect is different than in the case of the relaxed steady-state profile discussed above.Segregation does not affect the disrupted binaries as the source of candidate stars for collisions (with η b being an external parameter).The main effect here is that a shorter relaxation time allows injected stars to evolve their orbits prior to undergoing a collision.This is a mirror image of the effect that we have already seen in figure 5: A much greater fraction of the injected stars is destroyed in TDEs, at the expense of DCs.While the total DC rate is still substantial enough to be potentially observable, the suppression induced by mass segregation with respect to R DC (N H = 0) can be important for N H ≳ 5 × 10 3 , and reaches ∼ 50% at the highest values.
The MS-EMRI channel also terminates with stars being tidally disrupted at R T , but after a "gradual" descent, driven mostly by GW losses.Figure 7 shows the MS-EMRI rates calculated in the simulations as a function of N H , again for the cases of no collisions, collisions but no injected stars, and collisions along with stars injected at a rate of η b = 10 −5 yr −1 .
By construction, MS-EMRIs are a minor destruction channel in a Milky-Way like cluster.A star is likely to undergo an EMRI only when its GW loss time (Peters 1964;Hopman & Alexander 2005) is the shortest of all evolutionary time scales.As shown by Sari & Fragione (2019), in the absence of collisions, in a single-mass cluster T GW is shorter than T J 2B only in a relatively narrow region in (r, r p ) space, leading to a typical MS-EMRI rate of or about 1% of the TDE rate for Milky Way parameters.This is consistent with result R EM RI (N H = 0) ≈ 4 × 10 −8 yr −1 found for this case in the simulations.
It is evident that in the simulations without collisions mass segregation has an inhibiting effect on MS-EMRIs.This is a direct result of the dynamical friction aspect of segregation.The presence of black holes close to the SMBH opposes the GW loss effect on the orbital evolution of the lighter MS stars.A similar trend was naturally found by Aharon & Perets (2016), that heavier black holes reduce the EMRI rate of other, lower-mass compact objects.Unlike TDEs, which occur through a small number random scatterings, an MS-EMRI progresses through multiple, mostly quasi-spherical orbits.This becomes very unlikely when the drift drives a net energy gain through dynamical friction which counters the energy loss due to GW emission.Systematic inspiral toward R T is thus suppressed, by an order of magnitude for N H ≳ 10 4 and increasingly so (with some numerical scatter) at larger values.
As demonstrated in Balberg & Yassur (2023), collisions also suppress MS-EMRIs, since the probability of a star systematically inspiraling all the way from R col to R T before undergoing a collision is very small.Combined with the aforementioned effect of mass segregation, including collisions in simulations of a relaxed cluster drives the MS-EMRI rate to the level of the computational noise, and is essentially zero.
The MS-EMRI rate is resurrected if there exists a significant source of stars from disrupted binaries.These stars are injected relatively close to the SMBH, where T GW is comparable with the other timescales.If DCs are ignored, the MS-EMRI rate in such a scenario can rise to as much as a few 10 −6 yr −1 , on the same order of magnitude as TDEs (Sari & Fragione 2019).Once DCs are taken into account, this is no longer the case (since T col < T GW , but in the absence of the heavier black holes the MS-EMRI rate can still be as high as a few 10 −7 yr −1 (Balberg & Yassur 2023)).It is clear in figure 7 that a sufficiently large population of heavier black holes reduces this value significantly, again due to the opposing drift in energy induced by dynamical friction, which counters GW losses.It appears that even if the injection rate is high, the combined effect of DCs and mass segregation (assuming a nontrivial component of black holes) will render MS-EMRIs to be very rare.

The Spatial Distribution of DCs
The observational potential of DCs strongly depends not only on the total rate, but also on their spatial distribution.Collisions at R col will have a typical energy of ∼ GM 2 ⊙ /R ⊙ ∼ 10 49 erg, which could obviously power only a weak, short flare.On the other hand, collisions much closer to the SMBH could generate a light curve on par with standard supernovae in terms of the total energy, and perhaps similar to TDEs in terms of some of the observables (Balberg et al. 2013;Amaro-Seoane 2023).
Figure 8 shows the fraction of collisions which occur in the different radial bins in the simulations with no injected stars from binaries.The effects of mass segregation are reflected in the comparison between the three cases of N H = 0 (a single-mass population), N H = 4 × 10 3 and N H = 5 × 10 4 .
Recall that if there are no injected stars, the only source of stars for collisions are those which diffuse inward the cluster from r > R col .As long as the relaxation time at R col is longer than the collision time there, collisions tend to be confined close to R col , and this is the inferred distribution of DCs shown in in figure 8. Mass segregation and the presence of black holes does shorten T 2B (R col ) to some extent, as mentioned above, and so some penetration of MS stars to r < R col is possible.The fraction of collisions that occur in radii r < R col thus increases.This is especially evident for the larger abundance of black holes, N H = 5 × 10 4 .For example, the fraction of DCs which occur in the radial shell of 5 × 10 −4 R h < r < 10 −3 R h (1-2 mpc) is 4 times larger for N H = 5 × 10 4 than for N H = 0.
Interestingly, this trend is reversed at the smallest radii (r ≲ 10 −4 R h ), where for N H = 0 and even for N H = 4 × 10 3 the fraction of collisions is about 10 −3 , whereas for N H = 5 × 10 4 it is essentially zero.This When the injection rate is high but the population is limited to a single-mass, the DC distribution is almost uniform in the logarithmic radial bins.This result reflects the assumed log-normal distribution of the separation distances of the disrupted binaries, which carries on to a similar distribution of the orbits of the injected stars (when assuming a max = 1 au; see also in Balberg & Yassur (2023)).It is only the radial bins closest to R T which are naturally underrepresented in this distribution.
A second component of heavier black holes causes some quantitative, but not qualitative, changes in the spatial distribution of the DCs.First, the fraction of collisions closest to R T decreases, because some of the injected stars that travel through this region do experience sufficient gravitational scattering.Recall that the Balberg injected stars have eccentric orbits, so they are mostly scattered by black holes at r ∼ 100R T , where their number is not negligible.These captured stars can undergo a TDE before participating in a collision.The implication is that mass segregation not only enhances the TDE rate, but does so mostly at the expense of the most energetic collisions.The fraction of collisions that involves energies greater than 10 52 erg decreases from about 0.15 in a single-mass population (N H = 0) to about 0.10 for N H = 5 × 104 .Note that this effect comes on top of an overall reduction of about 50% in the total collision rate, so clearly mass segregation can constrain the observational potential of the highest-energy DCs.This preliminary conclusion requires more scrutiny, especially in terms of the assumptions regarding a max and the mass distribution of the stars, and also on the mass of the SMBH when applied to other galaxies.
The suppression of the fraction of collisions at the smallest distances from the SMBH is offset by an increased fraction of collisions at the opposite end, close to R col .The enhancement of the fraction in this particular range for large N H is again a result of the drift induced by the dynamical friction aspect of mass segregation, which enhances the partial probability of a captured star to reach larger distances from the SMBH.

CONCLUSIONS AND DISCUSSION
The violent destruction of stars is one aspect of the dynamics that shape the dense cluster of stars which surrounds an SMBH.While the deep potential of the SMBH dominates the instantaneous motions of the stars, and gravitational interactions among the stars drive the long-term evolution of their orbits, destruction preferentially removes stars with certain orbital properties.Destruction thus also impacts the steady-state profile.Moreover, given the large energies involved in stellar destruction, such events allow for unique observable transients, especially in an era of high-cadence surveys, such as ZTF1 (Bellm et al. 2019), ASAS-SN2 (Kochanek et al. 2017), and the upcoming ULTRASAT3 mission (Shvartsvald et al. 2023), and advanced GW astronomy, such as the Laser Interferometer Space Antenna 4 , (Amaro-Seoane et al. 2023).
There exist three potential channels for the destruction of MS stars near an SMBH.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 (a TDE).Such an event is followed by an optical flare, first through disruption and then as the debris accretes onto the SMBH (Gezari 2021).Stars sufficiently close to the SMBH also evolve their orbits through GW radiation energy loss and might inspiral and disrupt "gradually" as they approach and cross the tidal disruption radius, R T .In this case (an MS-EMRI), a possible electromagnetic flare will be accompanied by GW emission in the millihertz frequency band (Alexander 2017).Finally, a high-velocity collision between two stars may be completely destructive (a DC), initiating a supernova-like event in terms of the energy involved (Balberg et al. 2013;Amaro-Seoane 2023).
Since all destructive channels require stars to pass through the inner part of the cluster, and possibly evolve their orbital paths in this region, the event rates can be sensitive to the details of the stellar profile close to the SMBH.The focus of this work is on a first attempt to self-consistently study the event rate of the stellar destruction channels, including DCs, in a "mass segregated" profile (Bahcall & Wolf 1977;Alexander & Hopman 2009).Equipartition of energy in gravitational scatterings should create a segregated overabundance of heavier objects close to the SMBH.Such objects, and specifically black holes with masses of a few solar masses, shorten the relaxation time scales close to the SMBH, and accelerate diffusion in angular momentum and energy of the MS stars.Furthermore, the accelerated diffusion is also biased, with an average drift that increases the lighter stars' angular momentum and decreases their gravitational energy (following the convention that the potential energy of a bound object is positive).
The analysis is based a simplified numerical Monte Carlo code of stellar dynamics in galactic centers, based on the specific version by Sari & Fragione (2019) and the self-consistent tracking of DCs added by Balberg & Yassur (2023).Here the code was extended to allow for a two-mass population which includes N L and N H stars of masses 1 M ⊙ and 10 M ⊙ , respectively (the latter being black holes).Also added were gravitational scattering terms which reproduce the effects of dynamical friction between populations of unequal mass.The study was performed assuming a Milky Way SMBH of M • = 4 × 10 6 M ⊙ .The focus here is on the steady-state profiles and rates, which are typically reached within a few tenths of the global cluster relaxation time.The time-dependent evolution of the event rates is not analyzed here (and requires some further assumptions not included in the code), but it is noteworthy that rates of various destructive channels may evolve differently in time (see, e.g., the case of TDEs and EMRIs in Broggi et al. (2022)).The estimates derived here should there-fore be applicable to a Milky Way-like galaxy, whose age is approximately one relaxation time; realistic estimates for shorter time scales, and possibly also for longer ones, require further analysis.
It is constructive to begin with a general observation found here ( §4) that for almost any value of N H , gravitational scattering in the inner part of the cluster is dominated by the heavier objects.These objects are sufficiently overconcentrated in this region so that the dynamics are in the regime of "weak" segregation (Bahcall & Wolf 1977).This is the case even if the overall values of N L and N H correspond to the regime of "strong" segregation (Alexander & Hopman 2009;Preto & Amaro-Seoane 2010), in which the lighter objects are much more abundant and dominate the dynamics; see Linial & Sari (2022) for an analytical interpretation of this result.It is in this context that the trends regarding the rates of destructive events should be evaluated.
TDEs are dominated by MS stars orbiting close to R h (the SMBH radius of influence), which are scattered by other stars into very eccentric orbits that cross R T .Correspondingly, in a gravitationally relaxed cluster, the TDE rate is almost insensitive to the presence of a second component of black holes (a few 10 −6 yr −1 for the Milky Way case; see also in Panamarev et al. (2019)) .A minor enhancement of the TDE rate does arise because a small population of black holes near R h shortens the relaxation time there.On the other hand, the total DC rate in such a cluster does depend on the number of black holes in the population.This is the case because in the steady state of a relaxed cluster DCs deplete the MS stars in the inner region r < R col , since the collision time scale there is much shorter than the gravitational time.The presence of a nontrivial number of heavier black holes at r ≳ R col in the segregated profile creates a net drift effect, which lowers the rate that MS stars diffuse inward to R col , and correspondingly reduces the DC rate by as much as tens of percent for larger values of N H (≥ 10 4 ).Finally, the black holes dramatically reduce the number of stars that can inspiral toward the SMBH through a systematic orbital decay driven by GW emission, since it too is countered by the same drift.This effect tends to reduce the MS-EMRI rate even when DCs are ignored in the simulations (by as much as an order of magnitude for large N H ), and makes MS-EMRIs essentially nonexistent once DCs are taken into account.
The qualification of a "gravitationally relaxed cluster" used in the previous paragraph relates to the case in which there is a negligible addition of stars to the inner part of the cluster from binary tidal disruption.Such disruptions "inject" stars into tight, bound orbits close to the SMBH (Hills 1988), with high (e ≈ 0.99) eccentricities.These stars are primary candidates for destruction close to the SMBH.When assuming a singlemass population (Balberg et al. 2013;Balberg & Yassur 2023), in steady state, the collision time of these stars is much shorter than any gravitational time, so they are predominantly destroyed in DCs, at a rate of ∼ η b /2, where η b is the injection rate.However, in a mass segregated cluster, a significant presence of black holes close to the SMBH shortens the gravitational relaxation time experienced by the injected stars, and allows some of them to be scattered into tidally disruptable orbits.As a result, a significant black hole population directs some of the injected stars to the TDE channel, at the expense of DCs.The effect can be as a large as 50% of the DC rate for N H > 5 × 10 4 .Since some of the stars are injected very close to the SMBH, the MS-EMRI rate is never zero in this scenario, although it still remains susceptible to significant suppression (more than an order of magnitude) in the presence of a large black hole population, again due to the inherent drift outward mentioned above.
While TDEs and MS-EMRIs are unique to the immediate vicinity of R T , DCs can occur essentially anywhere in the range R T < r < R col , which in the Milky Way case spans 3 orders of magnitude in radius.The spatial distribution of collisions is thus important in terms of the corresponding energy, and therefore of the observational potential.As mentioned above, in the absence of injected stars representing disrupted binaries, most collisions occur close to R col .When the black hole population is large (N H ≳ 5 × 10 3 ), they shorten the relaxation time at R col and more MS stars actually penetrate inward to r ≲ R col , increasing the fraction of collisions there (while the total rate of collisions is reduced).On the other hand, for large values of N H the fraction of the most energetic collisions, those near R T is zero, rather than ∼ 10 −3 for lower values (and N H = 0).The net drift for large values of N H completely prevents MS stars from evolving their orbits to acquire very small values of r p , thus disallowing the most energetic DCs.
A similar, though only quantitative, effect arises when the rate of injected stars is large.Again the net drift will tend to reduce the fraction of collisions at the smallest radii (note that for eccentric orbits, collisions close to periapse are most likely).For the particular parameters examined here, which include binaries with original separation distances spread in a log-normal distribution between a b = 0.01 au and a b = 1 au, the fraction of the most energetic collisions with energies greater than 10 52 erg drops from 0.15 for N H = 0 to 0.10 for N H = 5 × 10 4 .This reduction is in addition to the to-Balberg tal collision rate dropping by about a factor of two for the same range of N H .A wider parameter survey is required, but it seems that mass segregation could reduce the rate of highest-energy (and most observable) collisions by tens of percents, if the typical fraction of heavy objects is significant.
While the main trends found here appear robust, there are certainly some important points for further consideration.First and foremost, mass segregation should be investigated with a full, realistic mass function.The two-mass system used here probably captures the main qualitative trends, but obviously quantitative issues require further scrutiny, especially if the fractions of other compact objects (white dwarfs and neutron stars) are large.Another issue is the possible outcome of a collision between an MS star and and a black hole.At the highest possible velocities the black hole probably passes through the star without destroying it (see the Appendix in Metzger & Stone (2017)), but close to R col , or even at r ≳ R col partial or full destruction may be possible.This detail should also be expanded to examine the possible role of low-velocity collisions of two MS stars in the outer parts of the stellar cusp.These more likely result in mass transfer, mass loss, and mergers (Freitag & Benz 2005), but multiple collisions may, however, be destructive.In any case, low-velocity collisions may accumulate to quantitative effects on the dynamics of the cluster (Sills et al. 2005;Dale & Davies 2006;Rose et al. 2023).In particular, they might increase the fraction of stellar-mass black holes (Rose et al. 2022), and generate a further complexity in the segregated profile.
Finally, it must be stressed that the Monte Carlo averaged scheme is only approximate, and does not cover the full range of possible effects.The dense stellar cluster generates coherent torques between slowly precessing orbits, namely resonant relaxation processes (Kocsis & Tremaine 2011, 2015).In principle, this effect is mostly relevant at radii beyond the DC regime (r ≳ 1000 au), since closer to the SMBH relativistic precession decouples individual orbits from the residual torques of the background stars (Bar-Or & Alexander 2016).This assessment should be reevaluated when considering a stellar profile which is modified by injected stars and mass segregation.Also not considered here are short-range, large-angle scatterings.While generally rare, such scatterings may be effective enough to alter the orbital evolution with respect to the prediction based on multiple small-angle scatterings.Both issues mentioned above are better addressed by designated N -body simulations, which will be reported in a separate work (G.Yassur and S. Balberg, 2023, in preparation).
7. ACKNOWLEDGEMENTS I wish to thank Pau Amaro-Seoane, Reem Sari, and Gilad Yassur for very useful discussions and comments.I also thank the anonymous referee for valuable comments and observations.

Figure 1 .
Figure 1.Top panel: the steady-state profile at time T2B(R h ) of the m = 1M⊙ MS stars (blue) and the m = 10M⊙ black holes (red), for the case {NL = 3.5 × 10 6 , NH = 5 × 10 4 }.GW losses, collisions, and injected stars are not included.The profile is presented through NL(r), NH (r), which are the number of stars in each radial shell between r/2 and r, as a function of r.Power-law fits are displayed for comparison with mass segregation theory.Bottom panel: the corresponding local value of ∆AH (r) (equation 26) found for the local values of NL(r) and NH (r) shown in the upper panel.

Figure 2 .
Figure 2. Same as figure 1, but for the case {NL = 3.96 × 10 6 , NH = 4 × 10 3 }.Top panel: the steady-state profile at time T2B(R h ) of the NL(r), NH (r).Bottom panel: the corresponding local value of ∆AH (r) (equation 26) found for the local values of NL(r) and NH (r) shown in the upper panel.

Figure 3 .
Figure 3. Top panel: the steady-state profile at time T2B(R h ) of the m = 1M⊙ MS stars for the case {NL = 3.5 × 10 6 , NH = 5 × 10 4 }.Shown is the number of stars in each radial shell, NL(r).GW losses are always included, and the different cases are no collisions (blue circles), collisions but no injected stars (violet squares), and collisions along with injected stars at η b = 10 −5 yr −1 (green diamonds).Bottom panel: the corresponding profiles of the 10 M⊙ stars for no collisions (red circles), collisions but no injected stars (magenta squares), and collisions along with injected stars at η b = 10 −5 yr −1 (maroon diamonds).Also shown for reference are the results for the code test (no GW and no collisions) from figure 1 (black X's).

Figure 4 .
Figure 4. Similar to figure 3 but for the case {NL = 3.96 × 10 6 , NH = 4 × 10 3 }.Top panel: the steady-state profile at time T2B(R h ) of the m = 1M⊙ MS stars.Bottom panel: the corresponding steady-state profile of the m = 10M⊙ black holes.Also shown for reference are the results for the code test (no GW and no collisions) from figure 2 (black X's).

Figure 5 .
Figure 5.The TDE rate averaged over the last 10% of the simulation time, as a function of the number of m = 10 M⊙ black holes, NH .Shown are the results for the simulation with GW losses but no collisions and no injected stars (black), GW losses and collisions but no injected stars (blue), and GW losses and collisions and stars injected at a rate η b = 10 −5 yr −1 (red).

Figure 6 .
Figure 6.The DC rate averaged over the last 10% of the simulation time, as a function of the number of m = 10 M⊙ black holes, NH .Top panel: the rates for the simulations with GW losses and collisions and stars injected at a rate η b = 10 −5 yr −1 (red).Bottom panel: the same with GW losses and collisions but no injected stars (blue).

Figure 7 .
Figure 7.The MS-EMRI rate averaged over the last 10% of the simulation time, as a function of the number of m = 10 M⊙ black holes, NH .Shown are the rates for the simulation with GW losses but no collisions (black), GW losses and collisions but no injected stars (blue), GW losses, collisions with stars injected at a rate η b = 10 −5 yr −1 (red).

Figure 8 .
Figure 8.The spatial distribution of the DCs averaged over the second half of the simulation time, for simulations with GW losses, collisions, but no injected stars.Shown is the fraction of collisions in each radial bin (between r/2 and r) as a function of r, for stellar populations with three values of NH (the number of the m = 10 M⊙ black holes).These are NH = 0 (black; a single-mass population), NH = 4 × 10 3 (red), and NH = 5 × 10 4 (blue).

Figure 9 .
Figure 9.The spatial distribution of the DCs averaged over the second half of the simulation time, for simulations with GW losses, collisions, and injected stars at a rate of η b = 10 −5 yr −1 .The fraction of collisions in each radial bin and the color-coding are the same as in figure 8.