Dissecting the properties of neutron star–black hole mergers originating in dense star clusters

The detection of gravitational waves emitted during a neutron star–black hole merger and the associated electromagnetic counterpart will provide a wealth of information about stellar evolution nuclear matter, and general relativity. While the theoretical framework about neutron star–black hole binaries formed in isolation is well established, the picture is loosely constrained for those forming via dynamical interactions. Here, we use N-body simulations to show that mergers forming in globular and nuclear clusters could display distinctive marks compared to isolated mergers, namely larger masses, heavier black holes, and the tendency to have no associated electromagnetic counterpart. These features could represent a useful tool to interpreting forthcoming observations. In the local Universe, gravitational waves emitted from dynamical mergers could be unraveled by detectors sensitive in the decihertz frequency band, while those occurring at the distance range of Andromeda and the Virgo Cluster could be accessible to lower-frequency detectors like LISA. The theoretical framework explaining the formation of a black hole–neutron start binary in isolation is well established, but little is known for the same phenomenon occurring in a denser environment as a star cluster. Using N-body simulations, the author models such complex interactions finding that the resulting binary systems may be quite different from isolated mergers and, with potentially no electromagnetic counterpart, they could be observed via gravitational waves with detectors sensitive at in decihertz frequency band.

T he third observational campaign, O3, operated by the LIGO-Virgo Collaboration will enlarge the family of gravitational-wave (GW) sources, currently comprised of 11 confirmed black hole (BH) binary and 1 neutron star (NS) binary mergers 1 . The loot of observations accumulated during O3 will hopefully include the first BH binary with component masses in the so-called upper mass gap, i.e., 50-150 M ☉ , and the first NS-BH merger. In fact, as reported in the gravitational-wave candidate event database (GraceDB, https://gracedb.ligo.org/), both events might have been already recorded. The observation of an NS-BH merger represents a crucial cornerstone in both GW astronomy and stellar evolution and dynamics. The GW signal emitted by this type of sources encodes information about the mass ratio of the binary, the BH spin, the NS compactness, and equation of state [2][3][4] . The next generation of GW detectors will allow us to chase these objects up to the dawn of Universe 5 , offering us the unique chance to follow them from the inspiraling phase to the merger. The formation of an accretion disc during the merger could trigger a kilonova event 6 and power a short Gamma-ray burst (sGRB) 7 , making NS-BH mergers promising multi-messenger sources. The coincident observation of GWs and a kilonova 8 , the detection of peculiar precession in the jets produced during the sGRB 9 , or the anisotropic emission of ejected matter 10 are some of the proposed signatures of a putative NS-BH electromagnetic (EM) counterpart. The properties of the EM emission depend on the binary properties close to the merger. A high eccentricity, for instance, can affect the amount of mass ejected, the mass accreted onto the BH and the angular momentum transferred 11 . The actual development of an EM counterpart is expected to depend likely on the mass ratio, the BH spin, and the NS equation of state 2-4,10-12 . For mass ratios smaller than '1/3 − 1/4, the NS undergoes tidal disruption inside the BH's innermost stable circular orbit (ISCO), thus preventing the EM emission 13 .
The scenarios behind NS-BH formation are still not fully understood. In the case of isolated binaries, binary population synthesis tools predict typical merger rates Γ = 9 − 100 yr −1 Gpc −3 (Gpc, gigaparsec) at redshift zero 14,15 , nearly circular mergers 16 with typical chirp masses 17 M chirp ¼ 3 M , and BHs with masses 18 strongly peaked around m BH~7 M ☉ and, in general, <20 M ☉ . The picture is loosely constrained for NS-BH mergers formed via dynamical interactions in star clusters (globular clusters, GCs), nuclear clusters (NCs), and galactic nuclei 19 , owing to the recent advances in our understanding of the physics governing the formation and retention of BHs and NSs. Indeed, star clusters might be able to retain long-lived BH subsystems [20][21][22] -which could persist at present in a number of Milky Way GCs 22,23 , the detection of NS in star clusters suggest natal kicks lower than previously thought 24,25 , and the discovery of GWs emitted by BHs as heavy as 30 M ☉ revolutionized our knowledge of stellar evolution for single and binary stars. Large number of BHs and NSs and the presence of heavy BHs can impact significantly the probability for NS-BH binary to form and, possibly, merge. Bridging reliable stellar dynamics, up-to-date stellar evolution recipes, and a detailed description of the last phases of binary evolution is crucial to assess the properties of dynamical mergers. Finding significant differences between dynamical and isolated mergers would represent a piece of crucial information to interpret future GW observations.
Here, we study the complex dynamical interactions, involving BHs and NSs in dense clusters, focusing on hyperbolic encounters between a binary, composed of a compact object and a stellar companion, and a single compact object, exploring two configurations: either the compact object in the binary is a NS and the third object is a BH (configuration NSSTBH), or viceversa (BHSTNS). Combining our simulations with observations of galactic GCs and NC in the local Universe, and with Monte Carlo simulations of GCs, we infer for dynamical NS-BH mergers an optimistic merger rate of~0.1-0.2 events per year and Gpc cube in the case of GCs, and~0.01 events per year and Gpc cube for NCs. Despite the small value, we find that dynamical mergers exhibit peculiarities that make them distinguishable from isolated mergers: chirp masses above 4 M ☉ , BH masses above 20 M ☉ , and the absence of associated EM emission if the BH is highly spinning and has a mass above 10 M ☉ . We calculate the associated GW emission showing that these mergers can be observed with detectors sensitive in the decihertz band and even with millihertz detectors like the laser interferometer space antenna (LISA), provided that they took place at distances typical of the Andromeda Galaxy or the Virgo galaxy cluster.

Results
Dynamical formation of NS-BH binaries in star clusters. To investigate this dynamical formation channel we exploit 240,000 direct N-body simulations that take into account up-to-date stellar evolution recipes for natal BH mass 26 and general relativistic corrections in the form of post-Newtonian formalism 27 . Figure 1 shows the trajectories of one of our simulations.
As detailed in the Method section, we vary both scattering parameters (binary semimajor axis and eccentricity, impact parameter and velocity of the third object) and environmental quantities (metallicity and velocity dispersion of the host cluster). To connect our results to real star clusters, we exploit the catalog of Milky Way GCs to find a tight relation connecting the GC velocity dispersion (σ), mass (M GC ), and half-mass radius (R GC ) (see the Method section for more details). To complement and support our simulations, we perform a deep analysis on the MOCCA Survey Database I 28 , a collection of over 2000 Monte Carlo models of GCs that span a wide portion of the phase space and represent GCs with present-day masses M GC~3 × 10 5 M ☉ and half-mass radii R GC~1 -3 pc. This sample allowed us to reconstruct the history of all NSSTBH and BHSTNS in 1298 models, and to derive an average scattering rate of dR sca /dt = 6.3 Gyr −1 for configuration NSSTBH and 245.4 Gyr −1 for BHSTNS. We find that a scattering results in the formation of a NS-BH in~1.27-1.59% of the cases, with the lower(upper) value corresponding to the NSSTBH(BHSTNS) configuration, but none of them merge within a Hubble time. On average, for configuration BHSTNS the scatterings occur at~0.01 times the cluster core radius R c , whereas for NSSTBH the scattering location is broadly distributed between 0.01 and 0.3 R c , but still well inside the cluster interiors. These scatterings occur late in the cluster life, usually several times the cluster half-mass relaxation time t rel . Figure 2 shows the distribution of the scattering time, t sca , normalized to t rel , calculated at 12 Gyr.
Despite the richness of information encoded in the MOCCA database, current models represent GCs and cannot be used to describe more extreme environments like NCs. Moreover, stellar evolution for BHs is not updated yet and the treatment used for close encounters does not include general relativistic corrections. Therefore, we use N-body simulations to perform a thorough investigation of this dynamical channel, using the analysis performed on MOCCA to: (a) compare with the scattering rate derived via N-body simulations, and (b) infer the time at which a scattering can occur.
The probability of NS-BH mergers. To identify potential mergers in our 240,000 N-body models, we need to associate to any newborn binary a formation time, t form . This is calculated through two quantities: the scattering time t sca , which we extract from the distribution of t sca /t rel derived for MOCCA models (Fig. 2), and the cluster relaxation time t rCL , which we extract from the distribution of values calculated for 157 galactic GCs (ref. 29 ) and 228 NCs (ref. 30 ). The NS-BH formation time is thus calculated as t form = t sca /t rel × t rCL .
After formation, we calculate the NS-BH merger time 31 t GW and, for each candidate, we draw 100 different values of t form , retaining only candidates for which the drawings is t form + t GW < 14 Gyr in at least 50% of the cases. Unfortunately, this requirement alone does not ensure that a merger can successfully take place. Indeed further interactions can soften and even destroy it if the NS-BH binding energy is lower than the mean kinetic energy of the environment 32 . The limiting value of the binary semimajor axis above which this can happen is called hard-binary separation a h = G(m 1 + m 2 )/σ 2 .
"Soft" binaries, a > a h , can be disrupted by strong encounters over an evaporation time 33,34 t evap = [σ 3 (m 1 + m 2 )/(32πG 1/2 2m * ρ a ln Λ] 1/3 , depending on the binary properties (m 1 , m 2 , a), the cluster velocity dispersion (σ), density (ρ), and average stellar mass (m * ), and the Coloumb logarigthm (ln Λ). Therefore, to shortlist merger candidates, we require the simultaneous fulfillment of three conditions: (i) t form + t GW < 14 Gyr, to avoid NS-BH binaries with delay times larger than a Hubble time, (ii) a > a h , to avoid soft NS-BH binaries, and (iii) t GW < t evap , to avoid NS-BH binaries that can be disrupted by further interactions.
Note that the delay time calculated this way don't account for the cluster formation time, t fCL , thus among all candidates satisfying simultaneously the three conditions above only a fraction f will satisfy also t fCL + t form + t GW < 14 Gyr. In our calculations, we assume that the majority of clusters form at redshift z~2 (ref. 34 ), corresponding to t fCL~1 0 Gyr. As shown in Fig. 3a, the fraction of merging NS-BH, p GW , increases at increasing the sigma, but depends poorly on the scattering configuration (NSSTBH or BHSTNS) and the metallicity. A rough limit to the merger rate in the local Universe for NS-BH mergers in clusters can be written as 35 : where Γ c is the merger rate per unit of time and cluster, ρ MWEG = 0.0116 Mpc −3 is the local density of galaxies 36 , and N c is the number of clusters in a given galaxy. The merger rate per cluster is given by where dR/dt is the rate of binary-single interactions and can be calculated combining N-body and MOCCA models as detailed in the Method section. To infer the number of binaries that at a given time coexist in the cluster, we exploit the 12 Gyr output of MOCCA models, in the case of NSSTBH configuration a GC hosts up to 4 NS-stellar binaries in 90% of the cases, and up to 7 binaries in the remaining 10%, while for BHSTNS GCs have <4 BH-stellar binaries in the 95% of the cases, and up to 12 in the remaining 5%. In our calculations, we assume N bin = 4 as a fiducial value. As shown in Fig. 3b, Γ c increases at increasing the velocity dispersion, is larger for the BHSTNS configuration at fixed sigma value, and larger for lower metallicities. In all the cases, Γ c is well described by a power law in the form Γ c = (σ/σ c ) α . Configuration BHSTNS displays a larger Γ c values due to the fact that they involve heavier binaries compared to NSSTBH, thus they are characterized by larger cross section and, thus, scattering rates.
Using the Γ c -σ dependence, we can exploit Eq. (1) to calculate the merger rate for Milky Way equivalent galaxies, namely those galaxies that share similar properties with our own, like a population of N c ∼ 200 metal-poor clusters with a relatively low velocity dispersion, σ~5-6 km/s. Under these assumptions we find a NS-BH merger rate with the two extremes corresponding to different metallicities. We find a remarkably well agreement with very recent results based on a sample of~140 Monte Carlo simulations of GCs (ref. 37 ). Regarding galactic nuclei, the mass and half-mass radius of the galactic NC (ref. 38 ) are M GC = 2.2 × 10 7 M ☉ and R GC~5 pc, respectively, thus corresponding to σ = 40-60 km/s. Under these assumptions, the merger rate for NCs in Milky Way analogs is Our estimates rely upon the assumption that all Milky Way-like galaxies harbor an NC, thus they represent an upper limit to the actual merger rate. Note that the merger rates for NCs and GCs are comparable for configuration BHSTNS, thus suggesting that NCs might account for 10-20% of the total population of dynamical NS-BH mergers. Figure 4 shows the variation of Γ c as a function of the cluster mass and half-mass radius. Our results are superimposed to the sample of observed GCs (ref. 29 ) and NCs (ref. 30 ). Only the heaviest and more compact NCs can sustain at least 1 event per Gyr. Table 1 summarizes all the models investigated, highlighting the number of exchanges-an exchange marks the formation of a NS-BH-, the number of mergers, and the number of possible EM counterpart out of 10,000 simulations. We stress that none of the observed clusters in the sample exhibit any evidence of a central massive BH (MBH), neither supermassive-for NCs-nor of intermediate-mass-for GCs.
Inside the so-called influence radius, R inf , it is possible to show that the relaxation time for MBH with masses in the range 10 4 -10 9 M ☉ is similar to t rel of clusters with masses in the same mass range, thus NS-BH formation could proceed similarly to NCs. However, deep into the influence radius, where the mass budget is dominated by the MBH itself and the velocity dispersion scales with R inf −1/2 , the relaxation time will increase as R inf 3δ/2 , being δ > 0 a factor that depends on the matter distribution around the MBH. For δ = 1, this implies that the relaxation time inside 0.1 R inf exceeds a Hubble time if M MBH > 10 6 M ☉ , thus indicating that the NS-BH formation channel explored here could be strongly suppressed in heavy galactic nuclei. The late evolution of a NS-BH binary formed around an MBH will depend on a number of processes. First, due to mass segregation, the binary will migrate inward, passing through regions at increasing density and velocity dispersion. This corresponds to a reduction of the hard-binary separation, meaning a larger probability for the binary to be disrupted if its hardening rate is not sufficiently large. Second, the increasing gravitational torque associated with the MBH can tidally rip apart the binary. Third, if the binary survive to both energetic scatterings and tidal torques, the reduced distance to the MBH could onset Kozai-Lidov oscillations 39,40 , which can excite the binary eccentricity up to unity potentially shortening its lifetime. Quantifying these effects for supermassive black holes (SMBHs) is challenging, owing to the fact that the physics regulating star formation and dynamics around an SMBH is still not fully understood. For intermediate mass black holes (IMBHs) in star clusters, this is even more difficult, owing to the lack of conclusive evidence of their existence and of a well-constrained formation scenario. For instance, recent numerical models suggest that IMBHs forming out of a sequence of stellar collisions are associated with clusters retaining only one or two BHs after the The ratio between the Chi-square and the number of degrees of freedom is χ 2 /NDF < 2.6 for NSSTBH and <0.5 for BHSTNS, assuming a 10% error associated to the measure of the scattering rate and the measure of the merging probability.
IMBH growth, thus limiting the probability for the NS-BH dynamical channel presented here to take place. Besides the formation of NS-BH mergers, we find in the case of configuration BHSTNS that the NS flyby can push the stellar companion on an orbit passing sufficiently close to the BH to trigger the stellar disruption and associated tidal disruption event (TDE). The probability for this to happen increases at increasing the velocity dispersion, being~1% for metal-poor and 1.5% for metal-rich clusters with σ = 5 km/s. The scattering rate for these events is larger for metal-poor systems, as here the BH mass is larger, resulting in a larger cross section and, thus, in a larger scattering rate. For values typical of Milky Way GCs, the resulting TDE rate is Γ TDE~( 2.4-4.2) × 10 −9 yr −1 , compatible with the value expected for TDEs triggered by BH binaries 41,42 . The 90% of disrupted stars have a mass < 0.5 M ☉ , thus possibly representing white dwarfs or low-mass main sequence stars.

Identifying dynamical NS-BH mergers with GW emission.
According to the forefront of binary stellar evolution recipes 18 , BH in isolated NS-BH mergers are expected to feature masses strongly peaked in the range 6.5-8.5 M ☉ and NS masses broadly distributed between 1.4-2 M ☉ , thus corresponding to chirp masses < 4 M ☉ . Figure 5a shows the chirp mass, M chirp , distribution for all our high-velocity dispersion dynamical models. We refer to models with σ = 100 km/s to discuss the general properties of dynamical mergers. This choice is motivated by the larger number of mergers for these models, which allow a more robust statistical investigation of merger mass distribution. Nonetheless, the overall distribution shown in the following does not differ from those at smaller σ values, although the latter are affected by a lower statistics. Mergers forming dynamically in our simulations show a non-negligible probability to have M chirp larger compared to the isolated channel. In configuration NSSTBH, up to 52%(32%) of mergers in metal-poor(rich) clusters have a chirp mass above this threshold. The percentage decreases for BHSTNS configuration but is still not negligible, being 14-17%, with the lower limit corresponding to metal-rich systems. A chirp mass above 4 M ☉ , thus represents the first clear distinctive mark of an NS-BH merger with a dynamical origin.
By definition, a large chirp mass indicates a large binary mass and, in the case of NS-BH binaries, this can indicate a large BH mass. In fact, the second characteristic mark of dynamical mergers is apparent in the BH mass distribution shown in Fig. 5b. For clarity's sake, we overlay to our predictions the same quantity inferred for isolated NS-BH mergers 18 . In metal-poor clusters, we find that >50% of NSSTBH and 17% of BHSTNS simulations lead to a merger involving a BH with mass m BH > 20 M ☉ . The percentage drops to 16 and 4%, respectively, for metal-rich clusters, due to the lower maximum BH mass set by stellar evolution for metal abundances close to solar (see the Methods section for further details about the initial BH mass spectrum).
However, we note that in comparison to isolated mergers, which predicts a narrow peak at m BH~7 M ☉ , dynamical mergers show a broad distribution even in the mass range 10-20 M ☉ , thus suggesting that the dynamical channel could dominate over isolated binaries already in this BH mass range. A BH mass above 10 M ☉ , thus represents the second distinctive mark of a dynamical origin for NS-BH mergers.
EM counterparts. One of the most interesting outcomes of a NS-BH merger is the possible development of an EM counterpart. This is associated with an accretion disc formed from NS debris during the merging phases. The disc can form only if the BH tidal field torns apart the NS before it enters the BH event horizon, a condition fulfilled if the NS tidal radius exceeds the BH's ISCO (ref. 39 ) where Z 1,2 are functions of the BH adimensionless spin parameter χ = a BH /m BH . Therefore, the merger will not feature associated EM counterpart if R tid /R ISCO < 1. Note that the opposite does not represent a conditio sine-qua-non for the development and detectability of an EM signal, as in the case R tid /R ISCO > 1, this depends on the geometry of the merger with respect to the observer and other potential observational biases. Figure 6 shows how the R tid /R ISCO ratio varies at varying m BH , assuming a NS radius 3 R NS = 12 km and mass m NS = 1-3 M ☉ , and different χ values. For mildly rotating BHs (χ~0.5), mergers meet the condition to enable EM emission only if the BH has a mass m BH < 3.8 M ☉ . In this case, neither the isolated channel nor the dynamical are expected to be prone to EM emission, being the mass of merging BHs larger than this threshold. For spin values similar to those inferred from LIGO observations 1 (χ~0.7), the threshold BH mass shifts to 5.2 M ☉ . In this case, we find 15 mergers out of 854 merger candidates, regardless of the configuration, with mass below this threshold, thus the probability to develop an EM counterpart is limited to <1.8%.
For highly spinning BHs (χ~0.9), instead, the BH mass threshold is 9.2 M ☉ . In this case, the vast majority of mergers in the isolated channel, especially for metal-poor environments, will fall in the region where EM counterpart is allowed, whereas dynamical mergers have a probability of 53.4% to fall outside this threshold, thus implying the impossibility for an EM counterpart to develop. Thus, the absence of a clear EM counterpart with a high-spin BH represents the third clear mark of a dynamical origin.
eventually lead to a merger, we find that at formation dynamical mergers are characterized by an extremely narrow eccentricity distribution peaked around unity. To explore whether some residual eccentricity is preserved when the merger enters the frequency bands of interest for GW detection, we calculate the evolution of the GW characteristic strain as a function of the frequency for all mergers, assuming that they are located in the local Universe, at a luminosity distance D L = 230 Mpc (redshift = 0.05). Note that this is compatible with the luminosity distance inferred for the two NS-BH merger candidates reported in the GRACE-DB. Figure 7a shows the eccentricity distribution as binaries cross the 10 −3 , 10 −2 , 10 −1 , 10 0 , and 10 1 Hz frequency bands. Note that a large fraction of binaries have e > 0.1 in mHz, i.e., in the observation band of spacebased detectors like LISA, but none of them have e > 0.1 when crossing the 1 Hz frequency threshold. Nonetheless, dynamical mergers appear to be potential multiband GW sources in the 0.01-1 Hz frequency range. Figure 7b shows the characteristic strain of mergers with a total merger time shorter than 10 5 yr in all our models. We overlap to the simulated sources the sensitivity curvein terms of characteristic strain-for low-frequency GW detectors (LISA, DOs 5 , ALIA 40,41 , and DECIGO 42 ) and high-frequency detectors (LIGO, KAGRA, and the Einstein Telescope). Decihertz observatories would constitute precious instruments to follow the evolution of these sources during the in-spiral phase down to the merger. In the same plot, we show an example for the signal of a merger taking place within the Andromeda galaxy, located at a distance of~779 kpc, or the Virgo galaxy cluster (~20 Mpc). Mergers occurring at distances between Andromeda and the Virgo cluster could spend enough time in the LISA band to be detected several years prior to the merger.   6 Limits on EM counterpart to a NS-BH merger. Ratio between the NS tidal radius (R tid ) and the BH innermost stable circular orbit (ISCO,R ISCO ) and the BH mass (m BH , x-axis) for mergers in all models with a velocity dispersion σ = 100 km s −1 . Colored regions enclose the limiting values of R tid /R ISCO assuming a spin parameter χ = 0.5, 0.7, and 0.9. Arrows and accompanying labels mark value of m BH above which an electromagnetic (EM) counterpart cannot develop, with smaller values corresponding to lower spins. Points represent all mergers in our models, with different symbols representing different sets. Only points lying above R tid /R ISCO < 1 can give rise to an EM counterpart.

Discussion
We modeled the dynamical formation of NS-BH mergers in massive clusters, exploring the phase space in terms of cluster velocity dispersion and metallicity, and assuming different configurations. We infer an optimistic merger rate of Γ GC = 0.1 yr −1 Gpc −3 for GCs and Γ NC = 0.01 yr −1 Gpc −3 for NCs, much lower than the rate inferred after the first two LIGO observational campaigns 1 (<610 yr −1 Gpc −3 ). This might indicate that dynamical mergers bring a little contribute to the overall population of NS-BH mergers. Nonetheless, our models suggest that dynamical mergers can exhibit distinctive marks potentially useful to interpret GW observations. While the isolated channel predict mergers with BH masses strongly peaked~7 M ☉ , and chirp masses <4 M ☉ , a non-negligible percentage of dynamical mergers could be characterized by BH masses above 20 M ☉ and chirp masses above 4 M ☉ . This difference has important implications for the development of an EM counterpart. For highly spinning BHs (spin χ = 0.9), the isolated channel suggests that all mergers have the possibility to produce coincident EM + GW emission. Conversely, in the dynamical channel up to 50% mergers have BH masses > 10 M SUN , sufficiently large to avoid the NS disruption outside the BH ISCO. We conclude that a dynamical merger might be uniquely identified if it fulfills simultaneously the requirements that: (i) the chirp mass exceeds 4 M ☉ , (ii) the BH mass exceeds 20 M ☉ , and (iii) an EM is absent if the BH spin exceeds 0.9. Dynamical NS-BH mergers appear to be promising multiband sources that might be observable with future decihertz detectors. Exceptional cases could be observed even with LISA, provided that the merger occurred at distances~0.7-20 Mpc, like in Andromeda or in the Virgo Galaxy cluster.

Methods
Comparing observations and numerical models. To compare our models with observations, we exploit the catalog of Milky Way GC (ref. 29 ), which provides, among other quantities, the distribution of velocity dispersion (σ), half-mass radius (R GC ), and relaxation time (t rel ), as shown in Fig. 8. Typical values for galactic GCs are, in general, σ~4-6 km/s, R GC~1 -5 pc, and t rel~1 Gyr. As shown in Fig. 9, the GC mass, half-mass radius, and velocity dispersion are connected by a tight relation in the form Using a least square fit we find α b = 1.14 ± 0.03, with an associated ratio between the χ 2 and the number of degrees of freedom χ 2 /NDF = 0.062.
We use Eq. (9) to convert the velocity dispersion, which is an input parameter in our N-body simulations, into GC mass and half-mass radius. We use the same strategy to compare our results with a sample of 228 NCs observed in the local Universe 30 , exploiting published mass and half-mass radius to calculate the velocity dispersion and half-mass relaxation time (Fig. 8). regularization 43 and includes General Relativity effects via post-Newtonian formalism 27 up to order 2.5. The choice of modeling a compact object paired with a stellar companion is twofold. On the one hand, stars constitute 90% of the total stellar population in a star cluster, making probable for them to be captured by a heavier object. On the other hand, since stars are lighter than compact objects it is energetically convenient for a binary to exchange components and increase its binding energy. Heavy compact binaries in star clusters can indeed form via a sequence of such interactions 45,46 , which indeed can contribute to the formation of NS-BH in star clusters 19 . To initialize the BH and NS masses, we sample the zeroage main sequence mass of the three components assuming a power-law mass function 47 , namely f(m * ) ∝ m * −2.3 . We calculate the remnant masses taking advantage of the SSE tool 48 for NSs and state-of-the-art mass spectra 26 for BHs. The latter are used because stellar evolution recipes for massive stars implemented in SSE are outdated. We show the BH mass spectrum adopted in Fig. 10a. Note that at low metallicities, the mass of the BHs extends to up to 60 M ☉ , while being much smaller for metal-rich progenitors. This is at the basis of the differences between the results obtained for different configurations with different metallicity values. We note that a smoother mass function would lead to a larger population of MBHs. This can lead to an increase of the probability for BHs in NS-BH mergers to have a mass larger than the value typical for isolated binaries (~7 M ☉ ). This, in turn, would increase the amount of dynamical mergers that might be clearly distinguishable from the isolated ones.
We assume that the three-body interaction is hyperbolic and in the regime of strong deflection, namely the outer angle between the incoming and outcoming direction of the scattering object is >90 degrees, and the maximum pericentral distance between the binary and the third object equals twice the binary semimajor axis. We restrict our analysis to strong scatterings, as these are the only capable to trigger an exchange between one of the binary components and the third object. To set the maximum semimajor axis a, we follow recent numerical results showing that this quantity divided by the binary reduced mass μ is proportional to the ratio between the host cluster half-mass radius and total mass 49 , namely a/μ = kR GC / M GC . The scaling constant k = 1/54 claimed in literature is typical of dynamically processed binaries, i.e., that underwent several dynamical encounters, while here we focus on binaries not fully dynamically processed. To mimic this assumption we set k = 10, and we calculate the R GC /M GC ratio through σ via Eq. (9). If a calculated this way is larger than the distance below which the star gets torn apart by tidal torques exerted by the companion, we set as maximum value allowed the hardbinary separation a h . The minimum binary separation is set as the maximum c Binary semimajor axis dstribution for model in which a neutron star-star binary scatters over a single black hole (NSSTBH). As in panel b but here we refer to the other configuration explored.
between 100 times the star's Roche lobe or 1000 times the ISCO of the compact object in the binary. This avoids the possibility that the star plunges inside the BH or is immediately disrupted before the scattering takes place. We initialize our scattering experiments basing our assumptions on previously published works focused on Monte Carlo modeling of GCs. To check the consistency of our assumptions, we compare the distribution of the binary semimajor axis normalized to the hard-binary separation in our models with σ = 5 km/s and MOCCA simulations, as shown in Fig. 10b, c. This quantity seems well suited to compare the two approaches, as it contains information about binaries orbital properties, via the semimajor axis and the component masses, and their hosting environment via σ. We find that, in general, the adopted distribution does not deviate dramatically from MOCCA results, thus providing an acceptable compromise that allows us to expand the study beyond the capability of MOCCA models. The initial binary eccentricity is sampled from a thermal distribution. Initial velocities of the binary and the single object are taken assuming a Maxwellian distribution of the velocities characterized by the star cluster velocity dispersion σ. We assume σ = 5, 15, 20, 35, 50, and 100 km/s, and two values for the stellar metallicity, either Z = 2 × 10 −4 , typical of old GCs, or solar values (Z = 0.02). As summarized in Table 1, our models can be divided into two main set, depending on the scattering configuration, each set is divided into two subsets, depending on the metallicity, and each subset is divided into six simulations sample depending on σ. Thus, we gather a total of 24 simulation samples each consisting of 10,000 simulations.
Calculating the scattering rate for N-body models. The interaction rate can be written as dR/dt = nσΣ, where n is the density of scattering particles, σ is the velocity dispersion, and Σ is the binary cross section For each simulation set, we calculate Σ by using the median value of a, e, m 1 , m 2 , and m 3 , calculated from the assumed initial distribution. The number density n of scattering particles depends critically on the amount of NSs and BHs left in the cluster. For BHs, we exploit our recent studies on BH subsystems in GCs (refs. 22,23 ). Using MOCCA models, we find that the typical density of the BH ensemble is comparable to the cluster density, n BH ≡ n ≈ M GC /(m * R GC 3 ). For NSs, instead, we consider the fact that segregation is mostly prevented by the BHs present in the cluster, whereas NS-to-total mass ratio for typical clusters is of the order of 0.01, a limit imposed by the standard initial mass function. Thus, we assume n NS ≡ 0.01n as an upper limit to the actual NS number density. Under these assumptions, we derive an optimistic estimate of the scattering rate dR/dt that, for σ = 5 km/s, results into: dR=dtð5 km=sÞ ¼ 2À4 Gyr À1 for NSSTBH; ð11Þ ¼ 150À400 Gyr À1 for BHSTNS: ð12Þ Note that these estimates fall in the range of values derived from self-consistent MOCCA models, for which we find 6.3 Gyr −1 and 254 Gyr −1 , respectively.

Data availability
The data sets generated during the current study are available from the corresponding author on reasonable request. The updated version of the ARCHAIN code used to carry out the N-body simulations is available from the corresponding author on reasonable request.