Introduction

Nuclear star clusters (NSCs) are dense regions consisting of stars and stellar remnants near the centers of most galaxies, including our Milky Way (MW). Most NSCs surround a central supermassive black hole (SMBH) with a mass between 10 6 and 10 9 M e (e.g., Kormendy & Richstone 1995;Ghez et al. 2000Ghez et al. , 2008;;Ferrarese & Ford 2005;Gillessen et al. 2009;Kormendy & Ho 2013). Due to its proximity, our Galactic center (GC) can serve as a unique place to investigate the conditions likely to occur at other galactic nuclei.

While the star formation process in the vicinity of an SMBH still remains a mystery, in particular with respect to the prevalence of binary formation, some studies indicate similarities to the field, where most massive stars (OBA spectral types) reside in a binary or higher-order configuration (e.g., Raghavan et al. 2010;Sana et al. 2012;Moe & Di Stefano 2017). Specifically, there are already three confirmed eclipsing binaries in the inner ;0.2 pc of the GC (e.g., Ott et al. 1999;Martins et al. 2006;Rafelski et al. 2007;Pfuhl et al. 2014), with possibly even more candidates (e.g., Gautam et al. 2019;Jia et al. 2019). Observations of the inner 0.02 pc find a dearth of young few million year old binaries, consistent with dynamical interactions (Chu et al. 2023) and suggesting a binary fraction close to 100% at birth for massive S-cluster stars (e.g., Stephan et al. 2016). Furthermore, X-ray observations have detected a large number of X-ray sources, implying a population of X-ray binaries (XRBs) or cataclysmic variables (e.g., Muno et al. 2005;Hailey et al. 2018).

On the theoretical side, Stephan et al. (2016) suggested that as many as 70% of binaries survive after a few million years of dynamical evolution at the GC. The dynamical interaction includes both frequent flybys from single passing stars that tend to unbind the binary (known as an evaporation process; Binney & Tremaine 2008;Rose et al. 2020), as well as interaction with the SMBH via the eccentric Kozai-Lidov (EKL) mechanism (Kozai 1962;Lidov 1962;Naoz et al. 2016). Further, Naoz et al. (2018) suggested that the existence of binaries may explain the peculiar properties of the stellar disk in the GC (Yelda et al. 2014). Moreover, merging binaries were suggested to form the G2-like object population (e.g., Witzel et al. 2014Witzel et al. , 2017;;Stephan et al. 2016Stephan et al. , 2019;;Ciurlo et al. 2020).

The evolution of massive binaries in the GC is affected by natal kicks that neutron stars (NSs), and possibly black holes (BHs), receive at birth (Fragione et al. 2019;Lu & Naoz 2019;Hoang et al. 2022). Observations of pulsar motion have revealed that NSs receive significantly large kick velocities on the order of hundreds of kilometers per second (e.g., Hansen & Phinney 1997;Lorimer et al. 1997;Cordes & Chernoff 1998;Fryer et al. 1999;Hobbs et al. 2004). It has been demonstrated that natal kicks can account for the misalignment between the orbital angular momentum and spin axes observed in pulsar binaries (Lai et al. 1995;Kalogera 1996Kalogera , 2000;;Kaspi et al. 1996;Kalogera et al. 1998). Studies have suggested that hypervelocity stars (e.g., Zubovas et al. 2013;Bortolas et al. 2017;Fragione et al. 2017;Lu & Naoz 2019), as well as extreme mass ratio inspirals (EMRIs) can be produced as a result of natal kicks disrupting massive binaries in the GC (e.g., Bortolas & Mapelli 2019;Lu & Naoz 2019;Hoang et al. 2022).

It is currently debated as to what the underlying stellar and stellar remnant distribution around SMBHs at the center of galaxies is. Theoretical arguments of a dynamically relaxed population yield ρ(r) ∝ r -α , with α = 3/2-11/4 (Bahcall & Wolf 1977;Alexander & Hopman 2009;Keshet et al. 2009). However, detailed measurements of the stars in our GC suggest a shallower distribution of α = 1. 1-1.4 (Gallego-Cano et al. 2018;Schödel et al. 2018). The distribution of compact objects (COs) at the GC, also known as the "dark cusp," has important implications for the dynamics in the vicinity of an SMBH. In particular, the CO distribution strongly affects the rate of gravitational wave (GW) events, tidal disruption events (TDEs), and the fraction of long-lived binaries in the GC (Amaro-Seoane et al. 2007;Alexander 2011;Pfuhl et al. 2014;Rose et al. 2020).

In this work, we study the evolution of binary stars orbiting the GC's SMBH and the resultant distribution of NSs and BHs. In Section 2, we describe the methodology to form single and binary BH and NS systems from massive stellar binaries, as well as the different natal kick prescriptions. In Section 3, we show that varying the initial stellar distribution steepens the postkick CO distribution, and that observations of the unseen mass in the GC allow us to constrain the initial stellar density profile. We also find that numerous high-energy events will be produced in this environment. In Section 4, we study the spatial and velocity distribution of COs near the GC, and suggest that the Nancy Grace Roman Space Telescope (Roman Space Telescope) may be able to distinguish between different kick prescriptions. We close with a discussion and conclusions in Section 5.

Methodology

In Hoang et al. (2022), Monte Carlo simulations of massive stellar binaries within 0.1 pc of the GC's SMBH were implemented to explore the effects of natal kicks on the binaries. In this work, we expand on these earlier simulations and explore the effects that varying the initial stellar distribution has on the overall CO density profile within the central parsec of the GC. See Figure 1 for a schematic of the methodology.

Birth Configurations

Each system begins as a hierarchical triple, comprising an inner binary of two main-sequence stars (m 1 and m 2 ) and an outer binary consisting of the orbit around an SMBH. The frame of reference is selected to be the invariable plane and we define the orbital parameters of the inner (outer) binary using the Keplerian elements for the semimajor axis, a 1 (a 2 ), eccentricity, e 1 (e 2 ), inclination, i 1 (i 2 ), argument of periapsis, ω 1 (ω 2 ), longitude of the ascending node, Ω 1 (Ω 2 ), and true anomaly, f 1 ( f 2 ). The inner and outer orbits are inclined to each other by a mutual inclination, i tot = i 1 + i 2 .

We define m 1 to be the more massive stellar binary member, such that it is always the first to undergo a supernova (SN) explosion. The mass distribution of m 1 is chosen from a Kroupa initial mass function ranging from 8 to 100 M e (Kroupa 2001). The mass ratio, defined as q = m 2 /m 1 , is chosen from a uniform distribution ranging from 0.1 to 1 (Sana et al. 2012). We set the mass of the SMBH at m • = 4 × 10 6 M e (e.g., Ghez et al. 2005;Gillessen et al. 2009).

The eccentricity distribution for the inner binary e 1 is uniformly distributed between 0 and 1, while the outer orbit eccentricity e 2 is taken from a thermal distribution (Jeans 1919). The mutual inclination i tot between the inner and outer orbit is distributed isotropically. The argument of periapsis, true anomalies, and the inner binary longitude of ascending node are selected from a uniform distribution between 0 and 2π.

We choose the outer semimajor axis a 2 to follow a powerlaw density cusp, n ∝ r -α , with a minimum semimajor axis of 500 au and a maximum of 1 pc. We vary α across the range of 0-3, in half integer increments and for each value of α, we run 1.5 million Monte Carlo simulations of the stellar binary orbiting around the SMBH.

The semimajor axis of the inner binary a 1 is determined from the period distribution ( ) dn dP P log 0.45 µ - (Sana et al. 2013), with the minimum and maximum value for a 1 selected for each system according to the following conditions.

1. First, we require that the stellar binaries' orbit pericenter be greater than 2 times the Roche limit of the system to ensure the stellar binary is not disrupted prior to the first natal kick:

The Roche limit of the stellar binary defined as

where R j is the radius of the star at mass m j and μ Roche,ji is the approximation of the Roche lobe radius (Eggleton 1983):

2. The upper limit for the a 1 distribution comes from ensuring that the system is hierarchically stable (Naoz 2016):

3. Finally, each triple system must also satisfy the following criteria of the stellar binary system not crossing the Roche limit of the SMBH before m 1 undergoes an SN explosion:

Binary Destruction

The initial stellar binaries can be destroyed either before or after the SN. We track merged and unbounded stellar binary members in our simulation. Therefore, our simulations consist of a population of binary and single-star systems orbiting the SMBH. There are three paths to destroying the binary before either star has gone SN.

1. SMBH Roche limit crossing. 32%-46%, from α = 0 to 3, respectively, of the initial stellar binary distribution (see Figure 2) did not meet the Equation (5) criterion. These evolve independently as single stars orbiting the SMBH.

In the statistical analysis below, we incorporate both the single-star population and the binary star population.

2. Stellar mergers induced by EKL. A fraction of stellar binaries will experience eccentricity oscillations induced by the EKL mechanism (Naoz 2016) and can become a merged stellar product before the first natal kick (e.g., Antonini et al. 2014; Prodan et al. 2015; Stephan et al. 2016; Wang et al. 2021). Following Hoang et al. (2022), we incorporate a simplified condition for which systems that have an EKL timescale shorter than general relativity (GR) precession may merge (or at least undergo mass transfer). We find that roughly 1%-6%, from α = 0 to 3, respectively, of the initial stellar binaries fall into this category and are excluded from undergoing an SN explosion in our simulations. 3. Unbinding via neighboring scattering interactions (evaporation). Weak gravitational interactions with nearby stars can unbind the binary over the evaporation timescale (e.g., Binney & Tremaine 2008): ( ) ( ) ( ) ( ) t r G r a m m m 3 32 ln , 6 p evap 1 1 2 s pr = L + where ( ) ln 5 L= is the Coulomb logarithm (Rose et al. 2020), m p is the average mass of the perturbing star, ( ) ( ) r Gm r 1

• s a =+ , and ρ(r) is defined below. Note that for simplicity, we ignore the eccentricity of the binary about the SMBH, since it will only change the timescale by a factor of a few (Rose et al. 2020).

We point out that we are testing a wide range of density profiles, α = 0-3, see Equation (8). However, observations of the GC suggest a shallow, core-like profile (α ∼ 1.1-1.4; Gallego-Cano et al. 2018;Schödel et al. 2018). Thus, following Gallego-Cano et al. (2018) and Rose et al. (2020), we adopt the evaporating population distribution to be with α = 1.3. In this case, most binaries have an evaporation timescale longer than the SN timescale for a range of separations about the SMBH (e.g., Rose et al. 2020). Only about 7%-9% of the remaining stellar binaries, from α = 0 to 3, respectively, will evaporate before the first SN. 5The remaining inner binaries can also be destroyed at a later time due to natal kicks or close encounters with the SMBH. Because m 1 is the more massive companion, it will undergo an SN explosion first. The first natal kick can disrupt the binary, leading to the formation of two separate orbits around the SMBH (m 1 -SMBH and m 2 -SMBH). If the binary survives m 1 ʼs natal kick, then we are left with a binary consisting of a CO and star orbiting the SMBH. This scenario may result in the formation of an XRB (Section 3.4). m 2 ʼs natal kick provides an additional way of destroying the binary and for the creation of GW mergers (Section 3.6). Either natal kick can also push the binary onto an orbit inside the SMBH Roche limit, resulting in the destruction of the binary.

Presupernova Evolution

Each star in the inner binary experiences mass loss due to main-sequence evolution. Between birth and the first SN, the inner and outer binary will expand due to mass loss but the outer binary will expand by a negligible amount because of the large mass of the SMBH. Using the rapid single stellar evolution code SSE (Hurley et al. 2000), we determine the time that each star becomes a CO and the mass prior to and following this event. By adopting adiabatic expansion, which conserves angular momentum, the inner binary semimajor axis immediately before the first SNe, -a 1,pre SN is:

where m 1,pre SN and m 2,pre SN are the masses of m 1 and m 2 immediately before the first SN.

Applying Supernova Kicks

We assume instantaneous SN kicks that are isotropically distributed. SN kicks for NSs are selected from a normal distribution with an average of 400 km s -1 and standard deviation of 265 km s -1 (e.g., Hansen & Phinney 1997;Arzoumanian et al. 2002;Hobbs et al. 2004). We adopt two different BH kick prescriptions due to observational uncertainties. In the fast BH kick prescription, BHs have the same kick distribution as NSs. In the slow-BH kick prescription, the BHs receive the same linear momentum kick as NSs (Bortolas & Mapelli 2019).

Recent studies suggest the possibility that NSs receive smaller birth kicks if they are formed through electron-capture supernovae (ECSNe), for a mass range of 6-10 M e (Miyaji et al. 1980;Nomoto 1987;Poelarends et al. 2008;Jones et al. 2016;Leung et al. 2020). However, only a window of approximately ∼0.2 M e in the 6-10 M e mass range actually undergoes electron-capture kicks (Doherty et al. 2017;Hiramatsu et al. 2021;Willcox et al. 2021;Stevenson et al. 2022). We also conducted two separate numerical experiments, considering ECSNe. One for which all the stars between 6 and 6.2 M e underwent ECSN, or ∼5% of the entire population. Thus, the inclusion of ECSNe at face value seems negligible. On the other hand, an extreme case is one in which all NS progenitors undergo ECSN. We chose an ECSN kick distribution that will lead to the maximum variation by taking a Maxwellian with σ = 30 km s -1 . The numerical experiment, in this case, is consistent with the slow-BH kicks, where the postkick density profiles have nearly identical slopes and similar amounts of unbound systems. This result is insensitive to the particular choice of ECSN kick distribution (e.g., Gessner & Janka 2018; Willcox et al. 2021;Stevenson et al. 2022). We, thus, omit the results from these experiments to avoid clutter throughout the paper.

m 1 will undergo an SN explosion first because it is the more massive companion. The SN kick is applied by adding a Cartesian velocity kick vector to the orbital velocity vector of m 1 and changing m 1 to the post-SN mass found with SSE. Following the first SN, there are two main scenarios: the inner binary survives m 1 ʼs SN kick or is disrupted by m 1 ʼs SN kick.

In the scenario where the inner binary survives, it can remain bound to the SMBH on an elliptical orbit or become unbound from the SMBH on a hyperbolic orbit. Just before m 2 ʼs SN kick, we adiabatically expand the orbits due to mass loss from m 2 using Equation (7).

For elliptical orbits, if the timescale between the first and second SN kick exceeds 10 times the orbital period, we randomly select the eccentric anomaly at the time of m 2 ʼs kick from a uniform distribution between 0 and 2π. Otherwise, we determine the eccentric anomaly by iteratively solving the elliptical Kepler equation using Newton's method. For hyperbolic orbits, we solve the hyperbolic Kepler equation using the HKE-SDG package (Raposo-Pulido & Peláez 2018) to find the hyperbolic anomaly at the time of m 2 ʼs kick. With either the eccentric or hyperbolic anomaly, we calculate the true anomaly, and then determine the Cartesian coordinates of the orbit immediately prior to m 2 ʼs kick.

With the calculated Cartesian coordinates of the orbit, m 2 ʼs SN kick is applied by adding the Cartesian velocity kick vector to the orbital velocity vector of m 2 and changing m 2 to the post-SN mass found with SSE.

In the second scenario where the inner binary is disrupted by the first SN kick, m 1 and m 2 form separate binaries with the SMBH. Then m 2 ʼs SN kick is applied by adding the Cartesian velocity kick vector to the orbital velocity vector of m 2 and changing m 2 to the post-SN mass found with SSE.

Interaction with the Supermassive Black Hole

If the separation of the inner binary (either progenitor or postkick binary) is larger than the SMBH's Roche limit, Equation (5), then the binary is disrupted and we follow the individual star's evolution. Further, binaries disrupted by natal kicks form two separate orbits around the SMBH (m 1 -SMBH and m 2 -SMBH). If the binary is disrupted by m 1 ʼs natal kick, then it is possible that m 2 will be on an orbit that will create a TDE (Section 3.5). On the other hand, if the binary is disrupted after the second kick, the result may lead to an EMRI (Section 3.3).

Normalization

Throughout this paper, we normalize the density distribution by the M-σ relation (Tremaine et al. 2002):

where M 0 = 10 8 M e and σ 0 = 200 km s -1 . In the rest of this paper, we refer to the numbers of NSs and BHs as expected from this normalization process. Here, we can also recognize a notable quantity called the "sphere of influence," which signifies the radius at which the gravitational potential is dominated by the SMBH. Equation (8) implies that this value is

, in our own GC.

Dark Cusp and High Energetic Phenomena Predictions

The Relationship between Dark Cusp and Stellar Density Distribution

The various dynamical processes described in the section above disrupt a significant fraction of binaries before the first SN. The natal kicks disrupt the majority of the remaining binaries, and by the end of the simulations, only a small fraction of all initial binaries remain bound to their companion (see Table 1 for details). The majority of the systems are single COs that are either orbiting the SMBH or unbound from the SMBH. In general, the COs do not remain in their initial position and are scattered.

There are two significant outcomes for a single or binary configuration postkick. One is if the binary or single remains bound to the SMBH, meaning the configuration postkick has Keplerian energy smaller than zero. The other is to become unbound to the SMBH; in other words, the Keplerian energy is larger than 0. Out of these ;20% are on a trajectory to escape the galaxy.

A schematic description of this result is depicted in Figure 3, where we show an example for α = 1.5, which is a core-like distribution similar to the one observed in our GC (e.g., Gallego-Cano et al. 2018;Schödel et al. 2018). Although only 3% of the NS progenitor population is formed within 0.1 pc, natal kicks move NSs that were originally located at a distance > 0.1 pc toward the GC, and ultimately 9.3% of NSs end up within 0.1 pc. At r = 0.107 pc, the average kick velocity (;400 km s -1 ) is equal to the circular orbital velocity and serves as a critical point for differentiating the behavior of the NS population in the two regions. 26.7% (64%) of the NSs Note. Initial systems in the single configuration are the result of binaries being disrupted by one of the processes that occur before the first natal kick, described in Section 2.2. initially formed with a semimajor axis less (greater) than 0.1 pc are unbound from the SMBH. The combination of this, along with a steepening of the NS number density within the 0.1 pc threshold, leads to a dense concentration of NSs within the 0.1 pc radius and a scarcity beyond it.

Below we highlight a few observational tests that can be used to constrain the CO progenitors' stellar distribution due to the unique nature of the GC and the natal kicks. Natal kicks efficiently move COs closer to the SMBH. Thus, observational constraints of the dark cusp may be used to constrain the initial stellar distribution.

Future observations can be used to constrain the dark cusp. The separation of a young binary at the inner 0.1 pc of the GC is sensitive to the underlying density profile, and measurements of such systems could be used to place constraints on the dark cusp (e.g., Alexander 2005; Rose et al. 2020).

The Effect of Progenitor Distribution on the Postkick Density and Eccentricity Distribution

Below we provide a detailed analysis of the NS distribution. The fast-kick BH distribution follows the NS distribution (only with a different normalization). The slow-BH kick results are described in Appendix B.

In Figure 4, we show the changes in the bound NS semimajor axis due to the kicks for three different density profiles, from extremely shallow (α = 0, left) to extremely steep (α = 3, right), as well as a core-like distribution closer to the observed distribution (α = 1.5, middle). As depicted, NS progenitors formed near 1 pc can move orders of magnitudes away from their birth positions while those formed in the nearby vicinity of the SMBH are scattered by only an order of magnitude or so. The shallowest initial density profiles (i.e., α = 0 and 0.5) contain the majority (;99%) of the NS population outside of 0.1 pc and so are significantly perturbed by the NS kicks and steepen dramatically within a = 0.107 pc. As the value of α increases, a larger fraction of NSs are initially within 0.1 pc of the SMBH, and so the increase in steepness is less susceptible to natal kicks, as further demonstrated in Figure 5.

In Figure 5, we show the NS progenitor (left panel) and bound NS (right panel) density distributions after the natal kicks. The bound NS density profiles are all steeper than their corresponding progenitor profiles. As the initial progenitor profiles increase in steepness, the corresponding amount of steepening in the bound profile decreases. The postkick density profile can be estimated analytically from the number of systems that become unbound to the SMBH. Conservation of particles implies that the main driver of the postkick distribution is the fraction of systems remaining. We provide details in Appendix A.

We apply a density criterion on the NS density profiles to constrain the expected initial stellar profile from observations of the precession of S0-2 caused by the unseen mass within S0-2's orbit (e.g., Do et al. 2019;GRAVITY Collaboration et al. 2022;Heißel et al. 2022). The upper limit is derived by assuming that all of the enclosed mass is NSs. In this case, an initial stellar profile with α < 3 is consistent with this constraint. However, if there are also WDs and stellar-mass BHs in this vicinity, assuming the typical population fraction of 0.26:0.014:2.3 × 10 -3 of WDs:NSs:BHs (Kroupa 2001), this means that about ;5% of the unseen mass is in NSs. Then, an initial stellar profile of α > 2 is incompatible with the mass constraints. Further observational measurements may be able to disentangle the mass fraction of NSs within S0-2ʼs orbit and provide a more stringent test on the initial stellar profile.

Lastly, the kicks may also significantly affect the NS eccentricity distribution, especially for extremely cuspy density profiles. Initially, all CO progenitors begin on a thermal eccentricity distribution. Note that a thermal distribution may not accurately describe the eccentricity distribution at the GC (Geller et al. 2019), but is used here as a proxy. In Figure 6, we display the changes in NS eccentricity due to the kicks for three different density profiles for a shallow (α = 0), intermediate (α = 1.5), and steep (α = 3) distribution. For shallow and intermediate initial stellar distributions (α = 0 and 1.5), the postkick eccentricity distribution follows the initial thermal distribution at lower eccentricities and drops slightly when e > 0.7. When considering the steeper distribution near α = 3, the orbits tend toward circularization, resulting in a higher proportion of orbits characterized by low eccentricities.

Extreme Mass Ratio Inspirals

EMRIs are GW emission events that take place when stellar-mass COs inspirals onto SMBHs. They are one of the prime science motivators of the future Laser Interferometer Space Antenna and other millihertz detectors ). The black line represents the line that NSs would remain on provided that there were no natal kicks. As depicted, for α = 0 and α = 1.5, the initial inner parts of the parameter space are almost entirely devoid of NS progenitors. Postkick, about ;6% (;8%) of the total progenitors' population for α = 0 (α = 1.5) have moved inward of 0.1 pc. For the α = 3 case, where the distribution is initially constant in a log initial , not much changes postkick. This behavior is further illustrated in Figure 5.

(e.g., Amaro-Seoane et al. 2017). A natal kick can drive a CO into the SMBH (e.g., Bortolas et al. 2017;Hoang et al. 2022). To estimate if a kick resulted in an EMRI we compare two timescales. One describes the characteristic GW decay timescale:

where c is the speed of light, e is the eccentricity of the object around the SMBH, postkick, and a is its semimajor axis (Peters 1964). The other timescale is two-body relaxation t relx , which is the result of weak kicks with other neighboring objects. On one hand, these kicks can result in EMRIs by changing the angular momentum of the orbit and driving it into the loss cone. On the other hand, the kicks can increase the angular momentum, yielding a more circular orbit and thus suppressing the formation of an EMRI. Following Amaro-Seoane et al. (2007), we classify an orbit to be an EMRI if t GW,EMRI < (1 -e)t relx is satisfied.

We convert the number of EMRIs in our simulations to the number of EMRIs within the sphere of influence, as expected from the M-σ relation. As shown in Figure 7, we find that EMRI formation is sensitive to the initial stellar distribution surrounding the SMBH. In particular, the expected number of EMRIs range from nearly zero EMRIs for a shallow cusp (α = 0) to 270 EMRIs for a steep cusp (α = 3). Considering a stellar profile that closely resembles the one observed in the GC (Gallego-Cano et al. 2018), we expect fewer than 10 EMRIs driven by natal kicks. For all initial stellar profiles, the majority of EMRI progenitors are formed within 10 -1 pc and are the result of NSs inspiraling onto the SMBH. We find that 98% (92%) are NS EMRIs and 2% (8%) are BH EMRIs, for α = 0 (3). We note that for α 2 the expected number of EMRIs from this channel is lower than the expected number of EMRIs from two-body relaxation (e.g., Hopman & Alexander 2006;Sari & Fragione 2019), and orders of magnitude lower than the expected number of EMRIs in SMBH binaries (Naoz et al. 2022;Naoz & Haiman 2023). For the extreme cusp case, i.e., α 2.5, the expected numbers of combined NS and BH EMRIs are comparable to the lower limit of the SMBH binary case. We suggest that extreme cusp profiles may also contribute to the revised stochastic background estimations presented in Naoz & Haiman (2023). We reserve this calculation for future studies.

X-Ray Binaries

Inner binaries that survive m 1 ʼs natal kick can have their orbital separation decrease. Following Naoz et al. (2016), we classify systems as XRB if the inner binary postkick pericenter drops below a Roche . We note that a binary system crossing the Roche limit is a necessary, yet insufficient condition for its evolution into an XRB. The transformation into an XRB also depends greatly on the specific evolutionary characteristics of the secondary star. Therefore, we can provide an upper limit on the number of XRB created from natal kicks.

We find that 3.3 × 10 -3 NS XRBs form per NS and 1.2 × 10 -3 BH XRBs form per BH in our simulations for all values of α other than α = 3. There is a decrease in the XRB fraction for α = 3 because there is a significant decrease in the number of initial stellar binaries (see Table 1). We find that the formation of XRBs is related to the properties of the inner binary and is independent of the binary's outer orbital parameters, such as distance away from the SMBH. From Figure 7, we expect nearly 400 to be formed within the sphere of influence due to natal kicks. From this, 94% of the XRBs are NS XRBs, and 6% are BH XRBs. The high abundance of X-ray sources observed at the GC (e.g., Hailey et al. 2018;Zhu et al. 2018) might be explained by these XRBs.

Tidal Disruption Events

TDEs occur when m 1 ʼs natal kick disrupts the stellar binary and pericenter of the m 2 -SMBH orbit drops below the SMBH tidal radius

where R * is the radius of the star and m * is its mass. We further require that r t is greater than the SMBH Schwarzschild radius and that m 2 passes within the tidal radius before its own natal kick to classify the system as a TDE. TDEs are a rare outcome of natal kicks acting on binaries. We find that no TDEs driven by natal kicks are expected to occur within the sphere of influence of the SMBH.

TDEs are expected to result via two-body relaxation processes (e.g., Rees 1988;Hopman & Alexander 2005;Fragione & Sari 2018;Madigan et al. 2018;Akiba & Madigan 2021) and in SMBH binaries (e.g., Chen et al. 2009Chen et al. , 2011;;Melchor et al. 2023;Mockler et al. 2023).

Inner Binary Gravitational Wave Mergers

The natal kicks can also direct the surviving inner binaries into regions of the parameter space where GR triggers a GW merger within a timescale shorter than the evaporation timescale at the GC. The inner binary GW merger timescale due to GR effects is (Peters & Mathews 1963):

~+ -

We label a system as a GW merger if t GW < t evap . In some cases, the EKL-induced eccentricity oscillations play a significant part in inducing a GW merger. If the EKL timescale is shorter than the GR precession timescale, we describe the EKL-induced GW merger timescale as:

where e 1,max is the maximum EKL-induced eccentricity and is estimated following Wen (2003). GW mergers are weakly dependent on the assumed initial stellar distribution and will result in 10-25 GW mergers within the sphere of influence of the SMBH.

Predictions for the Roman Space Telescope

Compact Object Distribution Beyond 1 pc

Consider a 3 Gyr star formation episode within 1 pc of the SMBH (consistent with Chen et al. 2023).foot_2 Within 1 pc, all NS and BH progenitors are initially orbiting the SMBH, but the natal kicks unbind a significant fraction of the COs from the SMBH potential, as described above (see Table 2). As expected, the percentage of COs that remain bound to the SMBH increases for a steeper initial stellar distribution.

As a test case, we focus on the α = 1.5 distribution. This density distribution is close to the GC observed stellar distribution (e.g., Gallego-Cano et al. 2018), and agrees with the constraints in Figure 5. With α = 1.5, 91% of the unbound systems are single COs (average speed of ;575 km s -1 ) and 2% are CO binaries (average speed of ;300 km s -1 ). The remaining 7% are ejected during their stellar lifetime due to their companion's natal kick and will undergo their own SN explosion outside the sphere of influence. These hypervelocity stars (average speed of ;600 km s -1 ) can briefly be observed for 10 6 -10 7 yr before becoming COs and contributing to the CO distributions. The combined gravitational potential of the MW will be significant enough to slow down the majority (∼70%) of systems unbound from the SMBH but bound to the MW potential with orbits scattered around the Galactic plane.

Here we focus on those COs that remain bound to the MW after 3 Gyr, and their potential detection using the Roman Space Telescope (Section 4.4).

We utilize the publicly available Python package for galactic dynamics galpy (Bovy 2015), to model a simple MW potential. We follow the orbits of all (bound and unbound to the MW) COs beyond the inner ∼1 pc. In Figure 8, we present the 3D distribution of a sample of COs ejected from the central parsec of the GC in a Galactocentric coordinate frame. We display the position of COs and a few selected orbits 100 Myr after the star formation event and within radial distances of 500 pc and 5 kpc. In both panels, the orbits cross within the inner regions of the GC, consistent with what is expected for being expelled from this region and falling back into the MW potential. The slow-BH kick prescription results in the BHs being concentrated closer to the GC than the NSs. The fast BH kick prescription results in the same density distribution of BHs and NSs, since by definition the fast BH kick prescription is matched to the observationally determined NS kick distribution.

The Relation between Galactic Latitude and Kick Prescription

Sweeney et al. ( 2022) recently analyzed the distribution of COs, including natal kicks, from the entire Galactic population (thin disk, thick disk, the stellar halo, and bulge). As suggested in Figure 8, the COs originating from the GC may also reach large distances. Below, we compare the GC population to the full Galactic population.

In Figure 9, we depict the Galactic latitude distribution of COs ejected from the central parsec of the GC after 3 Gyr and within a Galactocentric cylindrical radius of 8 kpc. Nearly 70% of NSs (left panels) and fast-kick BHs (right panels) are located at least a degree off the Galactic plane, whereas only 20% of the slow-kick BHs exhibit the same characteristic. Due to the strong natal kicks, the distribution of NSs and fast-kick BHs peaks near 3°off the Galactic plane. Notably, there is a subset of objects (∼21%) expelled from the central parsec that that is completely unbound from the MW. The distribution of slowkick BHs from the central parsec is concentrated within 1°. The decline beyond a few degrees is attributed to the comparatively lower velocities of the natal kicks. We propose that the GC population can be differentiated from the rest of the Galactic population. In Figures 9 and 10, we compare our results to the publicly available simulation data in Sweeney et al. (2022), and note that the COs from the Galactic population are ∼10 4 times more numerous. 7 As shown in Figure 9, the Galactic population's distribution of NSs and BHs are preferentially located at higher Galactic latitudes compared to the GC population. Thus allowing for the potential differentiation of these populations.

In Figure 10, we display the spatial and velocity distributions of the two BH populations ejected from the central parsec of the GC. As expected, the slow-kick BHs are more concentrated toward the GC and remain closer to the Galactic plane compared to the fast-kick BHs population (see left panel). We note that the Galactic population of NSs (and slow-kick BHs) in Sweeney et al. (2022) extends well beyond the GC distribution in both the x and z directions. This is because the natal kicks are occurring throughout the Galaxy and are not localized to the GC. The right panel shows that the Galactic population can reach higher velocities (maximum ∼ 870 km s -1 ) while the GC population attains slightly lower velocities (maximum ∼ 730 km s -1 ).

Compact Objects Unbound to the Milky Way

COs with velocities exceeding the escape velocity of the MW are unbound to the MW. 21% of all NSs within the central parsec are unbound to the MW by 3 Gyr (average speed of ;800 km s -1 , at 100 kpc from the center). As expected, the percentage of BHs unbound from the MW depends on the underlying kick prescription. The fast-kick BHs follow the NS percentage, while the slow-kick BHs only result in 2% of BHs Figure 9. Galactic latitude distribution of COs ejected from the central parsec of the GC with an initial stellar profile of α = 1.5 (solid lines) after 3 Gyr. Also displayed is the distribution of postkick COs from the entire Galactic population (dashed lines), as analyzed by Sweeney et al. (2022). In the top row, we plot the cumulative distribution functions for the NSs and BHs with the different kick prescriptions. In the bottom row, we plot the probability density functions for the COs. In both CO population sets, we limit the sample of COs in the distribution to be within a cylindrical radius of 8 kpc from the GC. The teal shaded region shows the range of Galactic latitudes expected to be observed by the Roman Space Telescope (e.g., Penny et al. 2019).

Figure 10. Spatial and velocity distribution of the BHs ejected from the central parsec and are bound to the Galaxy after 3 Gyr. In the left panel, we show the spatial distribution of the two BH populations in the Galactocentric frame. In the right panel, we show the heliocentric speeds of the BHs for both BH populations. Similar to Figure 9, we compare the spatial and velocity distributions to the results obtained by Sweeney et al. (2022). We only display the COs that will remain bounded to the Galaxy in Sweeney et al. (2022) for consistency.

being unbound to the MW (average speed of ;1650 km s -1 , at 100 kpc from the center).

Distinguishing between Kick Prescriptions with Gravitational Microlensing

The different BH natal kick prescriptions predict different distributions of COs as a function of Galactic latitude. Fast kicks result in an increasing number of BHs at increasing latitudes up to about 2°-3°off the Galactic plane, while slow kicks result in a decreasing number of BHs at increasing latitudes (Figure 9). Thus, if the number density of BHs as a function of latitude can be mapped, it would provide a way to determine the type of natal kicks BHs receive.

Gravitational microlensing can be used to measure the masses and velocities of dark massive objects in our Galaxy; for a detailed explanation, please see Hog et al. (1995), Miyamoto &Yoshii (1995), andWalker (1995). In brief, when a foreground object (such as a BH) aligns by chance with a background star along the observer's line of sight, the gravitational field of the foreground mass deflects the background star's light. The observer sees a transient brightening (photometric microlensing) and positional deflection (astrometric microlensing) of the background star. These two signals can then be used to measure the mass, velocity, and distance of the unseen lens. Gravitational microlensing has been proposed as a method to measure the mass distribution of COs toward the Galactic bulge (Gould 2000;Lam et al. 2020).

An isolated stellar-mass BH has recently been detected and characterized with microlensing, using ground-based survey photometry and Hubble Space Telescope follow-up astrometry (Lam et al. 2022;Mróz et al. 2022;Sahu et al. 2022;Lam & Lu 2023). This BH lens has been used to constrain the properties of natal kicks (Andrews & Kalogera 2022) as well as whether the progenitor system was binary or single (Vigna-Gómez & Ramirez-Ruiz 2023).

The Roman Space Telescope, NASA's next flagship mission scheduled to launch by 2027, will conduct several wide-field infrared surveys. Its Galactic Bulge Time Domain Survey (GBTDS) is designed to discover thousands of cold exoplanets via gravitational microlensing (Spergel et al. 2015;Penny et al. 2019). The notional design of the GBTDSfoot_4 will observe an area of ;2 deg 2 around 1.5°off the Galactic plane, avoiding regions within a degree of the GC.

In addition to exoplanets, the Roman Space Telescope could also detect and characterize hundreds of BHs via photometric and astrometric microlensing, as well as a comparable number of NSs if the astrometric precision is sufficient (Lam et al. 2020(Lam et al. , 2023; although see Sajadian & Sahu 2023 for a more conservative estimate based on more stringent characterization criteria). With its photometric precision, the Roman Space Telescope could also be used to study the population of COs in a statistical manner with photometric microlensing (Rose et al. 2022).

A detailed study is beyond the scope of this work, but we suggest that the Roman Space Telescope has the ability to study BH natal kicks and distinguish between slow and fast kicks. In particular, including an additional pointing toward the GC in the GBTDS would enable the measurement of the BH density as a function of latitude, and enable the determination of BH kick speed. We note that a broad range of other science cases would also be enabled by a field at the GC (Terry et al. 2023).

Discussion and Conclusion

NSs and perhaps even BHs receive large natal kicks during birth, with an expected average speed of 400 km s -1 (e.g., Hansen & Phinney 1997;Arzoumanian et al. 2002;Hobbs et al. 2004). Here we consider a GC population of massive stars (both single and binary), with different initial density distributions ρ ∼ r -α , with α ä [0-3]. The GC offers a unique opportunity to study the conditions surrounding SMBHs that probably take place in other galactic nuclei. Focusing on the postkick density distribution and comparing it to observations allows us to infer the initial stellar distribution at our GC.

The kicks in the vicinity of the SMBH may redistribute the orbital configuration of the COs around the SMBH, as well as unbind the binary itself. Adopting a kick distribution with an average kick velocity of 400 km s -1 implies that at ∼0.107 pc from the SMBH, the velocity dispersion around the SMBH is similar to that of the average kick magnitude. Thus, overall, we expect that kicks beyond this distance will more likely be unbound COs from the SMBH (see Figure 3), while those that remain bound (based on their initial orbital configuration; Lu & Naoz 2019), will migrate closer to the SMBH potential.

The natal kicks at the central parsec significantly affect the CO density distribution, i.e., the dark cusp. Here, we find that natal kicks steepen the resulting CO density profiles, with most of the steepening occurring within 0.1 pc for NSs and fast BH kicks. The natal kicks are efficient at driving stellar remnants from an initial semimajor axis beyond 0.1 pc, where the majority of the progenitor population is located, to bound orbits within 0.1 pc from the SMBH (Figures 4 and 5). This result goes beyond the previous studies by Bortolas et al. (2017) and Hoang et al. (2022), which were limited to values of 0.13 pc and 0.1 pc, respectively.foot_5 Even when considering slow-BH kicks, the resulting BH distribution still exhibits a steepening trend, although to a lesser extent (see Appendix B, Figure 11).

Using the predicted postnatal kick CO distribution, we constrained the initial stellar profile from limits on the unseen mass within S0-2ʼs orbit. Specifically, observations suggested that about ∼4000 M e resides inward to S0-2ʼs orbit (1000 au; Do et al. 2019;GRAVITY Collaboration et al. 2020;Heißel et al. 2022). Assuming that this unseen cusp is composed of stellar remnants such as stellar-mass BHs and NSs, we infer the initial stellar density distribution. Considering the standard population proportions of 0.26:0.014:2.3 × 10 -3 for WDs, NSs, and BHs (Kroupa 2001) within S0-2ʼs orbit, an initial stellar profile with α 2 leads to a CO density distribution that is incompatible with the mass constraints, as depicted in Figure 5.

This result is consistent with current observations of the stellar density distribution being close to unity. We note that if we adopt the unseen mass to be smaller than ∼3000 M e inward to S0-2ʼs orbit (e.g., GRAVITY Collaboration et al. 2020), we find a stronger constrain of the initial stellar density of α 1.5. The relation between the initial and final distribution is possible because two-body relaxation and collision effects have negligible effects on the final distribution at these stages (e.g., Rose et al. 2022Rose et al. , 2023)). Also, note that some theoretical arguments suggested that the unseen mass inward to S0-2ʼs orbit is consistent with the existence of intermediate-mass BHs (e.g., Generozov & Madigan 2020;Naoz et al. 2020;Strokov et al. 2023;Will et al. 2023;Zhang et al. 2023). In this case, the inferred initial stellar distribution may be even shallower.

In addition to the steepening of the CO density profiles, natal kicks naturally lead to the creation of EMRIs, XRB, TDEs, and binary GW mergers. From these, EMRIs are the most sensitive to the initial stellar profile, with a few hundred EMRIs expected for the steepest stellar profiles, as depicted in Figure 7. TDEs and binary GW mergers are less sensitive to the initial stellar profile, and we would only expect a handful of them. The numbers of EMRIs and TDEs expected from natal kicks are largely negligible compared to the two-body relaxation processes around a single SMBH (e.g., Alexander 2005;Hopman & Alexander 2005, 2006;Fragione & Sari 2018;Sari & Fragione 2019), and both are much lower compared to the expectation in SMBH binaries (e.g., Mazzolari et al. 2022;Naoz et al. 2022;Melchor et al. 2023;Mockler et al. 2023;Naoz & Haiman 2023). Unsurprisingly, XRB are unaffected by their distribution around the SMBH because the orbital properties of the inner binary directly affect the occurrence rate of XRBs.

A significant fraction of COs are unbound from the SMBH due to their natal kicks and may be potential microlensing events detectable by the Roman Space Telescope. As a proof of concept, we follow the unbound COs formed from an initial distribution of α = 1.5. This distribution is consistent with our aforementioned findings as well as with the observed GC stellar distribution (e.g., Gallego-Cano et al. 2018;Schödel et al. 2018). We follow these COs as they migrate throughout the Galaxy for 3 Gyr (see Figure 8). The adopted kick prescription is reflected in the spatial distribution of the COs in the Galaxy. Specifically, slow-kick BHs ejected from the GC are concentrated closer toward the Galactic plane, while fast-kick BHs and NSs are preferentially located at higher Galactic latitudes.

Lastly, we compared the GC CO distribution to the expected Galactic CO distribution and found that these two populations are potentially distinguishable. Particularly, the GC population is slightly slower (Figure 10) and presents a longer tail toward low Galactic latitude (Figure 9). The GBTDS expected field of view for the Roman Space Telescope is located in a Galactic latitude range that could possibly untangle the true underlying kick prescription for BHs. highest vertical black line in Figure 11. From this we can conclude an initial stellar profile of α < 3 is consistent with this criterion. If there are also WDs and NSs that make up a portion of the mass fraction within S0-2's orbit, with the typical population fraction from Kroupa (2001), then the upper limit is denoted by the lower vertical black line in Figure 11. Here an initial stellar profile with α < 2 is allowed from the mass constraint. With a mixed population of COs, both the NS and BH density profiles converge on an upper limit, regardless of the kick distribution.