Binary Natal Kicks in the Galactic Center: X-ray Binaries, Hypervelocity Stars, and Gravitational Waves

Theoretical and observational studies suggest that stellar binaries exist in large numbers in galactic nuclei like our own Galactic Center. Neutron stars (NSs), and debatedly, black holes (BHs) and white dwarfs (WDs), receive natal kicks at birth. In this work we study the effect of two successive natal kicks on a population of stellar binaries orbiting the massive black hole (MBH) in our Galactic Center. These natal kicks can significantly alter the binary orbit in a variety of ways, and also the orbit of the binary around the MBH. We found a variety of dynamical outcomes resulting from these kicks, including a steeper cusp of single NSs relative to the initial binary distribution. Furthermore, hypervelocity star and binary candidates, including hypervelocity X-ray binaries, are a common outcome of natal kicks. In addition, we show that the population of X-ray binaries in the Galactic Center can be used as a diagnostic for the BH natal kick distribution. Finally, we estimate the rate of gravitational wave (GW) events triggered by natal kicks, including binary mergers and EMRIs.


INTRODUCTION
Most galaxies, including our own, host a MBH at its center (e.g. Kormendy & Richstone 1995;Ghez et al. 2000Ghez et al. , 2008Kormendy & Ho 2013). Around these MBHs are the incredibly dense structures composed of stars and stellar remnants, called nuclear star clusters (NSCs).
Both theoretical and observational evidence suggest that binaries are common in NSCs. On the theoretical side, Stephan et al. (2016) showed that as many as 70% of the binaries formed in an initial star formation episode survive after a few Myrs. Furthermore, Rose et al. (2020) showed that these surviving young (few Myrs) binaries can have a wide range of orbital configurations and separations. Naoz et al. (2018) showed that a high binary fraction can explain many puzzling aspects of the young stellar disk in our own Galactic Center (e.g., Yelda et al. 2014).
On the observational side, recent works have suggested that binaries are prevalent in our Galaxy for massive stars (about 70% for ABO spectral type stars) (e.g. Raghavan et al. 2010;Moe & Di Stefano 2017). Thus, we can reasonably infer a comparably large binary fraction for massive stars in NSCs. Furthermore, the large number of observed X-ray sources in the central 1 pc suggest relatively large numbers of binaries containing BHs and NSs in the Galactic Center (Muno et al. 2005;Hailey et al. 2018). Moreover, there are three confirmed binaries in the inner ∼ 0.2 pc (e.g., Ott et al. 1999;Martins et al. 2006;Rafelski et al. 2007;Pfuhl et al. 2014), with possibly more candidates (e.g., Jia et al. 2019;Gautam et al. 2019).
Recently, Lu & Naoz (2019) explored natal kicks in hierachical triple systems, including stellar-mass binaries near an MBH (see also Hamers et al. 2018). They found that the natal kicks in these systems can lead to the inner binary shrinking as well as expanding, leading to many possible observable phenomena. In this work, we explore the dynamical and observable consequences of two successive natal kicks on a population of stellar binaries in the Galactic Center. Below we give a brief summary of the steps we take in our study and the structure of this paper: 1. We generate a population of massive stellar binaries. This is described in Section 2.1.
2. We evolve the stellar binaries up to the m 1 's natal kick and eliminate any binaries that undergo Roche Lobe overflow before m 1 's natal kick. This process is described in Section 2.2.
3. We apply the natal kick to m 1 and recalculate the orbits. We check the new orbits for any natal kickinduced observable phenomena (e.g. X-ray Binary, hypervelocity star). We then apply the natal kick to m 2 and repeat the process of checking for outcomes. We describe how we apply the natal kicks in more detail in Section 2.3. We describe the criteria we used to classify the outcomes and describe the results of our Monte Carlo simulations in Sections 3 and 4.
A simplified schematic of our steps is given in Figure  1, and a summary of the outcomes of our simulations is given in bar chart form in Figure 4.

METHODOLOGY
We employ simple Monte Carlo simulations to explore the outcomes of two natal kicks on massive stellar binaries orbiting a MBH within 0.1 pc. Inwards to this limit ∼ 1.5 × 10 5 , stars are expected to reside, estimated from the M-σ relation (e.g., Tremaine et al. 2002). Thus, this limit is used to estimate the effects of natal kicks on the system close to the MBH.

Monte Carlo Birth Distributions
We run a total of 300,000 Monte Carlo simulations of a binary star (the "inner binary") orbiting a MBH (the "outer binary") in a hierarchical triple configuration system. The inner binary comprises of main sequence stars with birth masses m 1 and m 2 , and the MBH has a mass of m MBH = 4 × 10 6 M . The inner (outer) orbit is characterized by the following orbital parameters: the semimajor axis a 1 (a 2 ), the eccentricity e 1 (e 2 ), the argument of periapsis ω 1 (ω 2 ), the longitude of the ascending node Ω 1 (Ω 2 ), and the true anomaly f 1 (f 2 ). The inner and outer orbits have a mutual inclination of i tot .
For our Monte Carlo simulations we set m 1 as the more massive star and therefore is always the one that undergoes the first supernova in the system. The birth value of m 1 is chosen from a Kroupa mass distribution ranging from 8 − 100 M (Kroupa 2001), and m 2 /m 1 = q, where the mass ratio q is chosen from a uniform distribution ranging from 0.1 − 1. Note that because of this mass distribution, some m 2 's will be white dwarf (WD) progenitors, and a small percentage will not undergo SN at all. Among the systems that are not eliminated in the pre-SN evolution (see Section 2.2), 53.2% of stars become NSs, 19.5% become BHs, WD 26.2% become WDs, and 1.2% never become a compact object (CO).
The inner (outer) orbital eccentricities are chosen from a uniform (thermal) distribution ranging from 0 − 1. We distribute the mutual inclination i tot isotropically. Further the arguments of periapsis, true anomalies, and the inner binary longitude of the ascending node Ω 1 are chosen from a uniform distribution between 0 and 2π. Since we are working in the invariable plane reference plane, Ω 2 = Ω 1 − π.
The outer semi-major axis a 2 is chosen from a number density distribution that goes as n ∝ r −2 with a minimum value of 1000 au, which is roughly the semi-major axis of the well-studied star S0-2 in the S-cluster in the Galactic Center, and a maximum value of 0.1 pc (see above). Note that current observations of late-type giant stars in the Galactic Center show a shallow core-like profile, rather than a cusp (e.g. Do et al. 2013). It remains debated whether this flat density distribution is representative of the rest of the stellar population (e.g. Genzel et al. 2003;Do et al. 2009;Bartko et al. 2010;Schödel et al. 2020). Theoretical models have found steeper cusps ranging from 1.5 -2.75, depending on the mass component (Bahcall & Wolf 1977;Alexander & Hopman 2009, e.g.). We expect the density distribution of binaries around the MBH to have an effect on some parts of our results. For example, we expect more EMRIs for a steeper cusp since supernovae closer to the MBH are more likely to result in an EMRI (Bortolas & Mapelli 2019b). On the other hand, we do not expect that our results for X-ray binaries will be greatly affected, since their formation depends more on the properties of the inner binary. We leave a comprehensive exploration of the effect of varying the density distribution to a future work.
We choose the inner binary semi-major axis a 1 from the period distribution dn/dP ∝ log(P ) −0.45 (Sana et al. 2013), with the minimum and maximum values of a 1 chosen for each system according to the following criteria. First we require the inner orbit pericenter to be greater than 2 times the Roche limit of the system: (1) Figure 1. A simplified schematic of our analysis process. Note that this schematic does not illustrate every single possible outcome of natal kicks on a binary.
The Roche limit a Roche,ij for either star in the inner binary 1 is defined as: where R j is the radius of the star of mass m j , found by the ZAMS mass radius relation: (Demircan & Kahraman 1991), and µ Roche,ji is defined by: (Eggleton 1983). To find the upper limit of a 1 , we require that each triple system is dynamically stable (Naoz 2016a): a 1 a 2 e 2 1 − e 2 2 < 0.1 Finally, we require that the inner binary does not cross the Roche limit of the central MBH, and thus eliminate systems that do not fulfill the criterion below: (e.g. Naoz & Silk 2014). We show the initial distributions of m 1 , m 2 , a 1 , a 2 , e 1 , and e 2 in Figure 2.

Pre-SN Evolution
For each star in the inner binary, we find the time at which it becomes a compact object and its mass before and after this transition with SSE (Hurley et al. 2000). The stars experience mass loss due to main sequence evolution between birth and the first supernova, and between the first supernova and the second supernova. Between birth and the first supernova the inner binary will expand due to mass loss (the outer binary will also expand, but by a negligible amount due to the large mass of the MBH). We find the inner binary semi-major axis immediately before the first supernova, a 1,pre−SN by adopting adiabatic expansion, which conserved specific angular momentum.: We note that we do not perform a numerical integration of the triple system either before the first kick or between the first and second kick. Eccentricity oscillations in the inner binary driven by the Eccentric Kozai-Lidov effect (EKL, Naoz 2016b), coupled with stellar evolution, can cause the inner stellar binary to become tidally locked, collide, or overflow its Roche Lobe, long before either of the stars undergo a supernova (e.g., Stephan et al. 2016Stephan et al. , 2019. We cannot account for all of these systems without a numerical simulation. However, we use an analytical criteria to identify systems where EKL causes the stellar binary to overflow its Roche Lobe. The EKL timescale on which the inner binary experiences eccentricity oscillations is approximately (e.g., Antognini 2015): This timescale needs to be shorter than the inner binary General Relativity (GR) precession timescale t GR in order for EKL to excite the inner binary eccentricity 2 . The GR precession timescale of the inner binary is defined as: where c is the speed of light. For systems with t ELK < t GR,inner and mutual inclination between 40 • and 140 • the maximal eccentricity e 1,max achieved by the inner binary during one Kozai-Lidov cycle can be approximated analytically (e.g. Miller & Hamilton 2002;Wen 2003;Fabrycky & Tremaine 2007;Liu et al. 2015). We approximate e 1,max using the method outlined in Wen (2003).
If the minimum pericenter distance achieved during t EKL (we also require that t EKL is less than the timescale to the first SN) is less than the binary Roche Limit a 1 (1 − e 1,max ) < a Roche , the stellar binary will undergo Roche Lobe overflow and mass transfer. We found that roughly 10% of initial stellar binaries fit these criteria. The subsequent binary evolution is out of the scope of our methodology. Thus, we eliminate these systems from our simulations and do not evolve them to the supernova stage. We note that the approach here uses the three-body Hamiltonian up to the quadrupole order, and cannot account for extreme eccentricities excited by the octupole term (e.g. Naoz 2016a). Moreover, for inner and outer eccentric orbits, the EKL mechanism can excite high eccentricity for wide range of inclinations (e.g., Lithwick & Naoz 2011;Naoz et al. 2013a;Li et al. 2014b,a). Therefore, our approach here eliminates conservative fraction of systems. However, we do expect that the binary orbital distribution from this approach will capture the distribution of the binaries post EKL, just before the first SN (as was highlighted in Rose et al. 2019).
We note that even though we eliminate these systems in the pre-SN stage, the statistics that we give later in the paper are with respect to the initial triple population we generated in Section 2.1. This is so that the event formation fractions we give later on in the paper can be simply multiplied by a star formation rate.

Supernova Kicks
We then apply the supernova kick to m 1 . We assume supernova kicks that are isotropically distributed in the inner binary orbital frame. We also assume instantaneous mass loss at the moment of the SN, using post-SN masses found with SSE (Hurley et al. 2000). For NS progenitors, we assume kicks drawn from a normal distribution with an average of 400 km/s and standard deviation of 265 km/s (e.g., Hansen & Phinney 1997;Arzoumanian et al. 2002;Hobbs et al. 2004). For WDs, which are present in the second set of kicks, we assume a normal distribution with an average of 0.8 km/s and standard deviation of 0.5 km/s (e.g., El-Badry & Rix 2018; Hamers & Thompson 2019). There is much uncertainty about whether BHs receive natal kicks, and if they do, what the distribution of kick velocities would be (e.g., Gualandris et al. 2005a;Willems et al. 2005;Fragos et al. 2009;Repetto et al. 2012). Thus, we adopt three different kick distributions for BH progenitors: 1. "Fast" BH kicks that are drawn from the same distribution as NS kicks. (2019b), "slow" BH kicks that assume BHs receive the same linear momentum as NSs (i.e. we draw a kick velocity from the NS distribution, then multiply this by (1.4 M /m BH ), where m BH is the mass of the BH progenitor, and 1.4 M is the typical NS mass).

No BH kicks.
We run 100,000 simulations for each scenario, for a total of 300,000 simulations. We present statistics from all three scenarios in Sections 3 and 4, and in Figure 4.
To apply the the natal kick, we simply add the Cartesian velocity kick vector to the velocity vector of m 1 , and change m 1 to the post-SN mass found with SSE. We then use new velocity vector and mass to recalculate new orbital parameters.
The two main scenarios resulting from this first kick are the "inner binary survived" and the "inner binary disrupted" cases (see Figure 1). In the case that the inner binary survives the first kick, we have two further branching scenarios. The inner binary can either stay bound to the MBH (ellitical orbit), or become unbound from the MBH (hyperbolic orbit). In the case that the inner binary is disrupted by the first kick, the component masses form separate binaries with the MBH (i.e. m 1 -MBH and m 2 -MBH). These binaries can either be bound or unbound. We then evolve m 2 up to its own Percentage of inner binaries disrupted/ survived after m1's and m2's kick. As expected, the percentage of inner binaries that survive both kicks is inversely correlated to the magnitude of BH natal kicks.
natal kick, and recalculate the orbits after m 2 's kick. Between m 1 's and m 2 's kicks we evolve the orbits as follows: 1. We adiabatically expand the orbits, due to m 2 's mass loss, using Equation 7 2. 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, from which we calculate the true anomaly. Using the newly calculated true anomaly with the Keplerian elements of the orbit, we calculate the Cartesian coordinates of the objects in the orbit.
3. For elliptical orbits, how we calculate the true anomaly at the time of m 2 's kick depends on the timescales of the system. If the timescale between the first and second natal kick is longer than 10 times the orbital period, we simply choose the eccentric anomaly at the time of m 2 's kick randomly from a uniform distribution between 0 and 2π. Otherwise, we solve the elliptical Kepler's Equation iteratively using Newton's Method to find the eccentric anomaly. From this we calculate the true anomaly, and then the Cartesian coordinates of the orbit.
After each kick we check whether the natal kick had induced any observable phenomena (e.g., X-ray binary, GW merger, etc.). We describe our criteria in Section 3 for the cases where the inner binary is disrupted, and in Section 4 for the cases where the inner binary survived. We also present results and statistics from our simulations for each outcome in the the following two sections. We include a simplified schematic of our process in  In the top (bottom) panels we show the results from after the first (second) natal kick. In the left (right) panels we show the orbital configurations where the inner binary survived (is disrupted by) the natal kicks. In the left panels, we use the crossed pattern to represent the binaries that survived the natal kick, but will be disrupted by crossing the MBH Roche limit upon pericenter passage (Section 4.3). We note that for the "Single Bound to MBH" cases, we have numbers greater than 1 because this number is taken relative to number of initial binaries, and each binary contains two singles.

INNER BINARY DISRUPTED BY NATAL KICKS
In the case that the inner binary is disrupted by either m 1 's or m 2 's natal kick, the result is two separate binaries with each composed of a stellar mass object orbiting the MBH (m 1 − MBH and m 2 − MBH). These orbits can either be bound (elliptical orbit), or unbound (hyperbolic orbit). However, natal kicks are less likely to result in unbound objects, as we find that the majority of objects remain bound to the MBH. See Table 2, for percentages of bound and unbound orbits as a function of stellar type.
The orbital configurations and observables resulting from these systems is described below. We also give rates of each orbital configuration and observables from our Monte Carlo simulations. These rates are summarized in bar chart form in Figures 3 and 4, which also break down the rates by BH kick distribution.

Single objects bound to MBH
Single objects bound to the MBH is the most common resulting orbital configuration in our simulations. We find 0.9-1.2 single objects per initial binary bound to the MBH after the first natal kick, and 1.1-1.3 after both natal kicks (see Figure 3). The vast majority of these orbit are long term stable orbits, with a very small number becoming EMRIs (Section 3.1.1), direct plunges (Section 3.1.2), and TDEs (Section 3.1.3).
In Figure 5 we show the distribution of NSs, BHs, and WDs with respect to the MBH after two successive natal kicks. Note that each triple is evolved in its own simulation, and each system undergoes natal kicks at different times, thus Figure 5 does not accurately show the distribution of a contemporaneous population of COs. However, Figure 5 is a good approximation of NS and BH distributions between a few 10 Myrs and ∼ 1Gyr after the star formation episode. This is because NS natal kicks occur between ∼ 10 − 45 Myr and BH natal kicks occur between ∼ 4 − 10 Myr, and the two-body relaxation timescale in the Galactic Center is ∼ 0.1 − 1 Gyr (Rose et al. 2020). WD natal kicks occur between 45 . For details about how we classify the different observables, see the following sections: X-ray binaries (4.1.1); HVB candidates (hypervelocity binaries, 4.2); GW candidates (4.1.2); HVS candidates (hypervelocity stars, 3.2.1); EMRIs (extreme mass ratio inspirals, 3.1.1); direct plunges (3.1.2); and TDEs (tidal disruption events, 3. and 1.4 × 10 4 Myr, so Figure 5 cannot be taken as representation of a comtemporaneous population of WDs at any given time. Figure 5 shows that after two successive natal kicks, the cusp of single NSs is steeper than the initial stellar binary cusp. This is also true for the BHs, if we assume fast BH kicks. WDs and BHs with no kicks more or less follow the initial stellar binary distribution. In fact, a steeper NS cusp is already present after just one natal kick. In Figure 6 we show the distribution of single NSs, BHs, stars, and binaries bound to the MBH after the first natal kick, for the slow BH kicks scenario. We see a steeper disribution of NSs than for Figure 5. Distribution of bound single compact object orbits after two successive natal kicks in the (a, 1 − e) plane (a and e are the semi-major axis and eccentricity with respect to the MBH). These objects were components of the stellar binaries that were disrupted by the natal kicks. For the NSs and WDs, we show data only from the "fast" BH kicks simulation set as the NS and WD kick distributions do not change between simulation sets and the results from any one set is fairly representative of the others. For the BHs, we show data from both the "fast" (in blue) and no BH kicks (in red) simulation sets to demonstrate two extremes. We plot the distribution (with respect to the MBH) of the initial stellar binaries that contained the compact object progenitors in yellow for comparison. In addition, we plot the distribution of EMRI progenitors (empty star symbol) and direct plunge progenitors (empty circles). BHs and stars, and a steeper distribution of binaries with a NS primary than for binaries with a BH primary. Thus, for a relatively young stellar population, before mass segregation has taken place, we expect a steeper cusp of NSs relative to the stellar and BH cusp. This is of course excepting the case in which BHs receive fast kicks similar to NSs, in which case the BH cusp will be as steep as the NS cusp after the first natal kick (not shown in figure).

Extreme Mass Ratio Inspirals (EMRIs)
Extreme Mass Ratio Inspirals (EMRIs) are gradual GW inspirals of stellar-mass COs onto MBHs, and are one of the main targets for the future Laser Interferometer Space Antenna (LISA) (Danzmann & et al. 2017).
The timescale on which a CO inspirals onto an MBH is (Peters 1964): In a dense environment such as the Galactic Center, we also have to take into account the effects of two-body relaxation on the potential EMRI. While two-body relaxation from other stars in the cluster can create EMRIs by pushing COs onto high eccentricity orbits, it can also push COs onto more circular orbits, such that t GW,EMRI increases to longer than a Hubble time. Two-body relaxation changes an orbit's angular momentum by an order of itself on a timescale of ∼ (1 − e)t relax , where t relax is: where ln Λ = 15 is the Coulomb algorithm, σ(r) is the velocity dispersion, given by (Kocsis & Tremaine 2011): and ρ(r) is the Galactic Center stellar density (Genzel et al. 2010): Thus, we require that t GW,EMRI < (1 − e)t relax to be classified as an EMRI (e.g., Amaro-Seoane et al. 2007a). EMRI formation is a rare event in our simulations. In particular, we find (3 − 8) × 10 −5 EMRIs per initial stellar binary 3 . This is roughly consistent with results from Bortolas & Mapelli (2019b), which studied EMRIs triggered by natal kicks on single stars in the Galactic Center and and found a rate of 10 −7 − 10 −4 EMRIs per SN. Assuming a star formation rate of 10 3 M Myr −1 in the Galactic Center (e.g. Lu et al. 2009), and an average binary mass of 30M (from our Monte Carlo simulations), we obtain an EMRI rate of (1 − 2.7) × 10 −9 yr −1 . Assuming a galaxy density of 0.02 Mpc −3 (Conselice et al. 2005), and that half of galaxies contain a MBH at the center, we calculate an EMRI volume rate of 0.01 − 0.027 Gpc −3 yr −1 from our simulations 4 . However, we note that this rate is likely a conservative one because of the minimum value of 1000 AU that we have chosen for the binary distribution around the MBH. As shown by Bortolas & Mapelli (2019a), stars closer to the MBH are much more likely than stars further away to end up as an EMRI after undergoing an SN. Thus, EMRI rates are very sensitive to the inner edge of the initial binary distribution. Currently, there several observed stars in the S-cluster with semi-major axes smaller than 1000 AU: S0-102/S55 with a semi-major axis of ∼ 900 AU (Meyer et al. 2012), S62 with a semimajor axis of ∼ 700 AU (Peißker et al. 2020a), and S4711, with a semi-major axis of ∼ 600 AU (Peißker et al. 2020b). Bortolas & Mapelli (2019a) found that a population of single stars located between 0.001 pc (∼ 200 AU) and 0.004 pc, modelled after the S-cluster stars, can generate up to a few 10 −4 EMRIs per SN. Thus, if S0-102/S55, S62, and S4711 are indicative of a larger population of stars and stellar binaries inside of 1000 AU, then the EMRI rate we have found here may underestimate the true rate from this mechanism by about an order of magnitude.
Interestingly, the three different BH kick cases give very similar EMRI formation rates. This is because in our simulations most EMRIs involve a NS inspiraling into the MBH (see Figure 5). Out of EMRIs from all three BH kick cases, ∼ 88 % are NS-EMRIs, ∼ 6 % are BH-EMRIs, and ∼ 6 % are WD-EMRIs.
Note that the WDs will not behave like a NS or BH EMRI because while we have a GW signal detectable by LISA while the WD is inspiralling, this signal will stop when the WD is tidally disrupted by the MBH and and replaced with an EM counterpart (e.g. Sesana et al. 2008).

Direct Plunges
Natal kicks can results in orbits with such a high eccentricity that they cause the orbiting object to plunge directly into the MBH instead of gradually inspiralling like EMRIs. Plunging orbits satisfy: where a and e are respectively the semi-major axis and eccentricity of the post-kick orbit, and R s = 2Gm MBH /c 2 is the Schwarzchild radius of the MBH (e.g. Amaro-Seoane et al. 2007b). Unlike EMRIs, which emit detectable GWs for thousands of orbits, creating a coherent signal, plunges emit a brief burst of GW radiation (e.g., Yunes et al. 2008;Berry & Gair 2013). Therefore direct plunges are not expected to be detected by LISA unless originating from the Milky Way Galactic Center (Hopman & Alexander 2005). We find, in our simulations, (1−5)×10 −5 direct plunges per initial stellar binary. The distribution of direct plunge-progenitors is shown in Figure 5.

Tidal Disruption Events (TDEs)
Tidal disruption events (TDEs) will happen when the pericenter of the m 2 − MBH orbit drops below the MBH tidal radius as a result from natal kick of m 1 . We also require that m 2 pass within the tidal radius before it undergoes its own natal kick to qualify as a TDE.
The tidal radius is defined as: where R * is the radius of the star, given by Equation (3), and m * is its mass (e.g. Rees 1988). We also require that r t lies outside of the MBH Schwarzchild radius so that the star is not directly swallowed instead of being tidally disrupted. TDEs are a a rare outcome in our simulations. We find 2 × 10 −5 TDEs per initial stellar binary for the "fast" and "slow" kicks scenarios.

Single objects unbound from MBH
Natal kicks can result in objects that are unbound from the MBH (i.e.) on hyperbolic orbits. In our simulations we find 0.27-0.36 single objects per initial stellar binary are unbound from the MBH after the first kick, and 0.35-0.45 after both kicks, making this the second most common orbital configuration resulting from our simulations (see Figure 3). Most objects are ejected from the inner 0.1 pc with velocities between 1000-10,000 km/s in our simulations, as shown by Figure 7. Despite deceleration by the combined gravitational potential of the MBH, nuclear cluster, and galaxy, objects ejected from the Galactic Center with velocities in excess of 1000 km/s should reach r > 8 kpc with ease, and retain their high velocities (e.g. Kenyon et al. 2008).
Unfortunately, hypervelocity NSs and BHs are not observable. However, a large number of m 2 's are ejected by m 1 's natal kick (see Figure 7). Since m 2 < m 1 and thus all the m 2 's are still stars at the time of m 1 's kick, some m 2 's may be observed as hypervelocity stars. We discuss this possibility further in Section 3.2.1.

Hypervelocity stars
If either the two members of the binary or one member of the binary becomes unbound to the MBH, it may be observed as a hyper-velocity star. A large number of stars have been detected with velocities larger than the the escape velocity from the Galaxy and with trajectories that are consistent with a Galactic Center origin (e.g., Brown et al. 2005Brown et al. , 2007Brown et al. , 2009bBrown et al. ,a, 2010Brown et al. , 2012Brown et al. , 2018Hirsch et al. 2005;Kollmeier et al. 2010;Boubert et al. 2018).
The leading channel to explain these systems has involved a binary star crossing the tidal limit and breaks up, ejecting in the process one of the members (known as the Hills mechanism, Hills 1988). The consequences of this mechanism have been explored in the literature in great details (e.g., Gualandris et al. 2005b;Bromley et al. 2006;Perets et al. 2007;Kenyon et al. 2008;Antonini et al. 2010;Sari et al. 2010;Kobayashi et al. 2012;Zhang et al. 2013;Rossi et al. 2014). Additionally, it has been suggested that SN kicks could also contribute to hypervelocity population (e.g., Zubovas et al. 2013;Lu & Naoz 2019;Bortolas & Mapelli 2019a). In particular, Zubovas et al. (2013) predicted a characteristic spatial anisotropy for hypervelocity stars created by SN kicks. As highlighted in Figure 7, we expect to have a population of hypervelocity stars as well as hypervelocity COs after the first natal kick. The first natal kick takes place between 4 and 45 Myr after the star formation episode, and as shown in the bottom row of Figure  7, over half of the ejected stars will not become COs until 100 Myr to 10 Gyr after star formation. Combined with the fact that the majority of these stars are ejected from the inner 0.1 pc with velocities between 1000 and 10,000 km/s, a large fraction of them can be observed as hypervelocity stars outside of the Galactic Center.

INNER BINARY SURVIVED NATAL KICKS
In the case that the inner binary survives m 1 's or m 2 's natal kick (or both), we have a stellar-mass binary orbiting the MBH. The orbital configurations and observables resulting from these systems is described below. We also give rates of each orbital configuration and observables from our Monte Carlo simulations. These rates are summarized in bar chart form in Figures 3 and 4, which also break down the rates by BH kick distribution.

Binaries bound to MBH
We find 0.1-0.29 binaries per initial binary remain bound to the MBH after one natal kick, and 0.04-0.17 after both natal kicks. As highlighted in Figure 3, assuming no kicks for the BHs increases the survivable rate of binaries around the MBH. Furthermore, the number of BH-BH binaries after two kicks is inversely correlated with the strength of BH kicks. For the "fast" BH kicks case, BH-BH comprise only 4.2% of the survived binaries. This percentage increases to 39.5% and 58.1% for the "slow" and no BH kicks case, respectively. As expected, due to high natal kick magnitudes, NS-NS binaries are relatively rare in our simulations. They comprise ∼ 1% of all surviving binaries for the fast BH kicks case, and ∼ 0.2% for the two other BH kick cases. For the unbound stars and binaries, we also show the time at which m2 will become a CO (bottom row). Note that this time is measured with respect to the initial stellar formation and not with respect to m1's SN. We also show the population of unbound NS X-ray Binaries in the third column.
In Figure 9 we show the distribution of orbital parameters for the surviving CO binaries for the no BH kicks case, compared to the initial binaries. We choose to highlight the no BH kicks because this is the case with the most surviving binaries. In general we find that the surviving binaries have larger semi-major axes (a 1 ) than the initial binaries. Surprisingly, despite the large NS kick magnitudes, the a 1 distribution for binaries containing at least one NS (i.e NS-WD, NS-NS, and BH-NS), peaks at a smaller value than for binaries that don't (i.e. BH-WD, BH-BH). This may be because binaries that start with a smaller a 1 have a better chance at surviving high NS natal kicks. Furthermore, the binaries containing at least one NS generally have much large eccentricities (e 1 ) than binaries that don't. We find similar trends with the semi-major axis and eccentricity of the binaries' orbits around the MBH (a 2 and e 2 ). In general, NS containing binaries have slightly larger a 2 and e 2 than the other types of CO binaries.
These surviving binaries can undergo EKL eccentricity excitations from the MBH which may result in GW mergers (e.g., Hoang et al. 2018;Stephan et al. 2019). We explore this effect in Section 4.1.2. We also find X-ray binaries among the surviving binaries, which we discuss in Section 4.1.1.

X-ray Binary
Natal kicks can actually tighten the inner binary, causing its pericenter to drop below the binary Roche Limit, creating a X-ray binary. Note that we will only get X-ray binaries after m 1 's kick and not after m 2 's kick, since one member of the inner binary needs to be still a star. We classify a system as a potential X-ray binary if the inner binary post-kick pericenter is below a Roche , given by Equation (2). In our simulations we find (6.3 − 9.9) × 10 −3 X-ray Binaries per initial stellar binary, for the no BH kicks and "fast" BH kicks, respectively.
In Figure 10 we show the distribution of X-ray binaries for the different kick models. Unlike Figure 8 here we divide the notation to BH X-ray binaries and NS Xray binaries. As highlighted in Figure 10, the no BH kicks case has zero BH X-ray binaries. Thus, at face value, this implies that detecting BH X-ray binaries at the center of the galaxy may imply that BH may have some kicks 5 .
Interestingly, our simulations predict some X-ray binaries binaries ejected from the GC (see Figure 7). Figure 8. Distribution of triples that remain after natal kicks in the (a1,1 − e1) plane (top plots, a1 and e1 are the inner binary semi-major axis and eccentricity), and the (a2,1 − e2) plane (bottom plots, a2 and e2 are the outer binary semi-major axis and eccentricity). We show X-ray binaries created by m1's natal kick in red, stable triples in blue, and binary disruptions (these are binaries that survive the natal kicks but will be disrupted by the MBH at the next pericenter passage) in purple. A fraction of the binaries in the stable triples will become GW sources.
These comprise about 20% of all X-ray binaries, for all BH kick scenarios. See Section 4.2 for further discussion.

GW Candidates
Natal kicks can shrink the inner binary orbit or increase the eccentricity, or nudge the binary into a part of the parameter space where EKL can excite the binary eccentricity, both scenarios leading to a shorter GW merger time. However, before a binary can merge, it can become unbound due to cumulative encounters with other stars in the cluster (Binney & Tremaine 2008). This takes place on an "evaporation" timescale of: where ln Λ = 15 is the Coulomb algorithm, m b = 1 M is the average mass of background stars, and σ(r) and ρ(r) are given by Equations 12 and 13 respectively. Note that our orbits are very eccentric, however, as shown by Rose et al. (2020), the eccentricity may only change the evaporation timescale by a factor of a few. There are two pathways to merger, depending on whether EKL plays a significant part. In the case where t EKL > t GR,inner (given by Equations 8 and 9), i.e. systems where GR effects dominate over EKL effects and the inner binary does not experience EKLinduced eccentricity excitations, we calculate the inner binary gravitation wave merger timescale following Peters (1964): Figure 9. Distributions of CO binaries bound to MBH after both natal kicks for the no BH kicks case. The initial triple probability densities are shaded in gray for comparison. Notably, CO binaries containing at least one NS tend to be more eccentric than binaries with a BH primary due to the large NS kicks.
If t GW < t Evap , we denote the system as a GW merger candidate.
In the case where t EKL < t GR,inner , we analytically estimate the maximal EKL-induced eccentricity, e 1,max using the method detailed in Wen (2003), and estimate the EKL-induced GW merger time as: (e.g. Liu & Lai 2018;Randall & Xianyu 2018). Using the above method, we find (0.17 − 1.2) × 10 −3 GW merger candidates per initial stellar binary (see Figure 4 for break down by BH kick scenario). Note that this is likely an underestimate because the calculation of e 1,max in Wen (2003) uses the three-body Hamilton up to only the quadrupole order, and may miss octupole effects, which increases the merger rate (Hoang et al. 2018).
Using the same method we used to calculate the volume rate of EMRIs as in Section 3.1.1, we find a volume rate for GW mergers facilitated by natal kicks and EKL of 0.06-0.4 Gpc −3 yr −1 .

Binaries unbound from MBH
Natal kicks can leave the inner binary bound while unbinding it from the MBH. In our simulations we find 0.007-0.01 binaries ejected from the MBH per initial stellar binary after one natal kick, and 0.005-0.006 after both natal kicks, making this the least common orbital configuration resulting from our simulations (see Figure  3). After both natal kicks, these unbound CO binaries are not observable. However, similarly to the ejected single objects from Section 3.2, binaries ejected from the Galactic Center after m 1 's natal kick (these binaries contain one CO and one star), are observable. In the third column of Figure 7 we show the post-kick positions and velocities of these binaries after the first natal kick for the "slow" BH kick case. Most of these binaries are ejected from the GC with velocities between 1000 -10,000 km/s. These binaries are ejected from the Galactic Center between 4 and 45 Myr after the star formation episode, and ∼ 38% of these binaries contain a star that will not evolve to the CO phase until 100 Myr -10 Gyr after the star formation episode, giving up to a few Gyr's window to observe them.
Interestingly, ∼ 15% (17 %) (22 %) of the ejected binaries are X-ray binaries, for the no BH kicks ("slow" BH kicks) ("fast" BH kicks) scenario. For the "fast" BH kicks scenario, 86% (14 %) of the ejected X-ray binaries have a NS (BH) primary. For both of the other BH kicks scenarios all the ejected X-ray binaries have a NS primary.
Of the ejected binaries that are not X-ray binaries, ∼ 76% (89%)(96%) have an NS primary, for the no BH kicks ("slow" BH kicks) ("fast" BH kicks) scenario. The rest have a BH primary.

Binary Disruption by MBH
In our Monte Carlo runs, we make sure that the inner binary birth distribution does not cross the Roche Limit of the MBH (Equation 6). However, sometimes the natal kick will push the inner binary onto an orbit around the MBH that does violate the condition given in Equation (6). In this case the inner binary will be disrupted at the next pericenter passage with the MBH. In our simulations we find (0.68 − 3.9) × 10 −2 binaries disrupted by the MBH, per initial stellar binary.
Typically, the heavier member of the inner binary is captured by (ejected from) the MBH for a bound (unbound) outer orbit, and vice versa for the lighter member (Hills 1988;Yu & Tremaine 2003;Sari et al. 2010;Kobayashi et al. 2012). Previous works have shown that captured COs and stars from these binary disruption events can results in EMRIs, plunges, or TDEs (e.g. Sari & Fragione 2019;Miller et al. 2005). We check whether this is the case for any of our binary disruptions. The inner binary is tidally disrupted by the MBH when it passes inside the tidal radius: The newly captured object's orbit has a pericenter equal to this tidal radius. Its new semi-major axis is given by: We show the probability density of the post-kick inner binary semi-major axis, and the the probability density of the captured objects in Figure 11. Because the binaries disrupted by the MBH are wide by definition (their distribution peaks at ∼ 100AU), the captured objects have relatively large semi-major axes with respect to the MBH (their distribution peaks at ∼ 1 pc). As a result, we find that none of the captured objects in our simulations result in EMRIs, plunges, and TDEs. However, they might have an effect on the stellar cusp in the inner parts of the Galactic Center (e.g. Fragione & Sari 2018).
The other binary member is typically ejected from the MBH at high velocities (e.g. Hills 1988), and may add to the number of hypervelocity stars in our simulations. However, this out of the scope of this paper.

DISCUSSION
We have conducted an analysis of the consequences of natal kicks in massive binaries at the Galactic Center. When a star undergoes supernova to become BH or NS it is expected to receive a large natal kick (e.g. Hansen & Phinney 1997;Lorimer et al. 1997;Cordes & Chernoff 1998;Fryer et al. 1999;Hobbs et al. 2004Hobbs et al. , 2005Beniamini & Piran 2016).
We run three large sets of Monte-Carlo simulations (100, 000 for each set), calculating the dynamical outcome after supernova kick. The three sets correspond to three different kick distributions for BH progenitors: fast kicks (for which the kicks are similar to NS kicks), slow kicks (for which the kicks are normalized by the BH mass), and no kicks. The application of the kicks are done by simple vector analysis. The consist form of the equations are given in Lu & Naoz (2019).
Generally, after each kick there are two main outcomes, where the inner binary becomes unbound or stays bound. The second kick may also unbind a binary in the latter case. After two kicks took place we find that most of the binaries were disrupted, e.g., between 94.5% − 80.1%, for the fast to no kicks models respectively. See Table 1 for the full details. However, even in the cases of inner-binary disruption, the majority of the systems remain bound to the MBH (see Table  2). These consequences yield the following predictions and observational signatures: • Natal kicks create X-ray binaries, at a rate of (6.3 − 9.9) × 10 −3 per initial stellar binary. In general we find that faster kicks create more Xray binaries. In particular, we find that no BH X-ray binaries are created if we assume no BH kicks. This suggests that detecting BH X-ray binaries in the Galactic Center implies some degree of BH kicks. However, other formation channels can create BH X-ray binaries in the absence of kicks (e.g., Naoz et al. 2016;Stephan et al. 2019).
• Hypervelocity stars and COs are a common outcome. We find that after the two kicks have taken place, a fairly significant fraction of the systems became unbound to the MBH (slightly more than 30% of NSs, and between 32.3% − 1.7% of BHs for fast to no kicks, respectively). See Table  2, for details. After the first kick, a number of stars are ejected from the inner 0.1 pc with velocities between 1000 -10,000 km/s. Some of these Figure 10. Distribution of X-ray binaries in the (a1, 1 − e1) plane (bottom three panels), and the kernel density estimation (KDE) of a1 (top three panels). X-ray binaries with a NS (BH) primary is shown in blue (black). We show X-ray binaries from all three BH kick distribution cases. We find that the set receiving "fast" BH kicks resulted in the most BH-star X-ray binaries, with the set receiving no BH kicks resulting in no BH-star X-ray binaries at all. Furthermore, "Fast" and "slow" BH kicks result in different semi-major axis distributions for BH-star X-ray binaries, as shown by the KDEs in the top panels. Thus, the type and semi-major axis distribution of a statistically significant population of X-ray binaries in the GC can be a diagnostic for the hitherto uncertain distribution of natal kicks received by BH progenitors.
stars may be observed as hypervelocity stars outside of the Galactic Center.
• Hypervelocity binaries are created in addition to the ejected COs and stars, about 20% of which are X-ray binaries. These binaries are ejected from the inner 0.1 pc after the first kick with velocities between 1000 -10,000 km/s. In particular, after the first natal kick, we find ∼ (7.4 − 9.5) × 10 −3 binaries (CO + star) ejected from the inner 0.1 pc. These may be observed as hypervelocity binaries and hypervelocity x-ray binaries outside of the Galactic Center.
• Kicks result in a steeper distribution of the bound single NSs and BHs about ∼ 0.04 pc from the MBH compared to the initial distribution. The progenitors began on Bahcall & Wolf (1976) density profile, however, after the two kicks both the NS and the BHs (with fast kicks) result in a steeper distribution inwards to 0.04 pc. At this distance from the MBH, orbital velocity around the MBH is comparable to the average kick velocity of ∼ 200 km s −1 . Beyond this radius, we find that the distribution of the NS and BHs (with fast kicks) have a shallower density distribution compared to their birth distribution. BHs with no kicks and WD (without kicks) follow the initial density distribution. This behavior is depicted in Figure 5.
• Kicks result in a slightly shallower binary density distribution about MBH. Unlike the steeper segregation of singles toward the MBH, the bound inner binary post kicks result in a somewhat shallower density distribution, with a long tail of systems reaching large densities. Additionally, the inner binary semi-major axis distribution after two natal kicks has a greater spread and peaks at a larger value, relative to the initial distribution.
• Exotic events from our simulation include binary GW mergers, EMRIs, direct plunges, and TDEs. Of these, the most common are binary mergers, which occur in our simulations at a rate of ∼ (0.17 − 1.2) × 10 −3 per initial binary. EMRIs, direct plunges, and TDEs are much rarer events in our simulations. For each of these, we find a few times 10 −5 events per initial binary. However, we note that our estimate of the EMRI rate is likely a conservative one due to assumptions about the inner edge of binary distribution around the MBH, see Section 3.1.1 for further discussion.
In summary, we showed that natal kick have a significant effect on the density distribution of COs, and of binaries in the galactic center. For example, it produces a steeper cusp of NSs about the MBH. Additionally, natal kicks naturally give rise to hypervelocity stars, hypervelocity binaries and X-ray binaries. Interestingly, we found that X-ray binaries can serve as a discriminator between BH kicks distribution models. Lastly, natal kicks can also result in a non-negligible rate of exotic events such as TDEs, LIGO, and LISA sources.
BMH and SN acknowledge the partial support from NASA No. 80NSSC19K0321, ATP-80NSSC20K0505 and partial support from the NSF through grant No. AST-1739160. BMH thanks the University of California Office of the President Dissertation Year Fellowship. SN thanks Howard and Astrid Preston for their generous support.