Supernova Kicks and Dynamics of Compact Remnants in the Galactic Centre

TL;DR

This study uses simulations to explore how supernova explosions affect the orbits of stars and remnants near the galactic center, revealing that black holes tend to stay bound and may form a dense cluster.

Contribution

It provides the first detailed simulation-based analysis of supernova kicks' effects on stellar and remnant dynamics in the galactic center.

Findings

Black holes remain bound in the galactic center after supernova kicks.

Supernovae cause lighter stars to have higher eccentricity and inclination.

A significant fraction of remnants could become X-ray binaries.

Abstract

The Galactic Centre (GC) is a unique place to study the extreme dynamical processes occurring near a super-massive black hole (SMBH). Here we investigate the role of supernova (SN) explosions occurring in massive binary systems lying in a disc-like structure within the innermost parsec. We use a regularized algorithm to simulate 30,000 isolated three-body systems composed of a stellar binary orbiting the SMBH. We start the integration when the primary member undergoes a SN explosion, and we analyze the impact of SN kicks on the orbits of stars and compact remnants. We find that SN explosions scatter the lighter stars in the pair on completely different orbits, with higher eccentricity and inclination. In contrast, stellar-mass black holes (BHs) and massive stars retain memory of the orbit of their progenitor star. Our results suggest that SN kicks are not sufficient to eject BHs from the GC. We thus predict that all BHs that form in situ in the central parsec of our Galaxy remain in the GC, building up a cluster of dark remnants. In addition, the change of NS orbits induced by SNe may partially account for the observed dearth of NSs in the GC. About 40 per cent of remnants stay bound to the stellar companion after the kick; we expect up to 70 per cent of them might become X-ray binaries through Roche-lobe filling. Finally, the eccentricity of some light stars becomes > 0.7 as an effect of the SN kick, producing orbits similar to those of the G1 and G2 dusty objects.

Tables

Table 1. Table 1: Initial conditions and differences between the three sets of simulations.

	set A	set B	set C
$M_{2,\mathrm{ZAMS}}$	TH, $0.1-150M_{☉}$	$0.9\times M_{1,\mathrm{ZAMS}}$	TH, $0.1-150M_{☉}$
CW disc	$0.04-0.13$ pc	$0.04-0.13$ pc	$0.001-0.13$ pc
$N$	10,000	10,000	10,000
$\Delta t_{\mathrm{sim}}$	1 Myr	1 Myr	1 Myr

Table 2. Table 2: Fraction of objects that become unbound from the SMBH.

Object type	set A	set B	set C
all	6%	6%	5%
BHs	0%	0.4%	0.04%
NSs	3%	24%	2%
II star	10%	0.05%	8%

Table 3. Table 3: Fraction of stellar binaries that do not break 1 Myr after the SN kick.

Object type	set A	set B	set C
all	38%	48%	39%
BHs	69%	80%	71%
NSs	2%	11%	2%

Table 4. Table 4: Fraction of stellar binaries that may turn on as X-ray emitters.

	set A	set B	set C
min($R_{\ast }$)	0	0.1%	0.01%
max($R_{\ast }$)	27%	76%	29%

Equations

$$\mathbf{v}=-v_0(\hat{\text{ı}}\sin \beta +\hat{\text{ȷ}}\cos \beta ),$$

$$\sin \beta$$

$$\cos \beta$$

$${\mathrm{E}}_{i+1}=\mathcal{M}+e\sin {\mathrm{E}}_i,$$

$$r=a(1-e\cos \mathrm{E}).$$

$$v_0^2=G{M}_{b,i}\left(\frac{2}{r}-\frac{1}{a}\right),$$

$${\mathbf{v}}_{\mathbf{k}}=v_k(\hat{\text{ı}}\cos \omega \cos \varphi +\hat{\text{ȷ}}\sin \omega \cos \varphi +\hat{\mathbf{k}}\sin \varphi );$$

$${\mathbf{v}}_f$$

$$v_f^2$$

$$=$$

$${\mathbf{h}}_f=\mathbf{r}\times {\mathbf{v}}_f,$$

$$|\mathbf{r}\times {\mathbf{v}}_f{|}^2=G{M}_{b,f}{a}_f\left(1-e_f^2\right),$$

$$|\mathbf{r}\times {\mathbf{v}}_f{|}^2=r^2\left[v_k^2{\sin }^2\varphi +{\left(v_k\cos \omega \cos \varphi -v_{\mathrm{orb}}\sin \beta \right)}^2\right];$$

$$\cos \eta =\frac{v_0\sin \beta -v_k\cos \omega \cos \varphi }{{\left[v_k^2{\sin }^2\varphi +{\left(v_k\cos \omega \cos \varphi -v_0\sin \beta \right)}^2\right]}^{1/2}}.$$

$${\mathbf{v}}_s=\frac{M_{1,f}}{M_{b,f}}{\mathbf{v}}_k-\frac{\Delta M{M}_2}{M_{b,f}{M}_{b,i}}\mathbf{v}$$

Full text

††pubyear: 2016††pagerange: Supernova Kicks and Dynamics of Compact Remnants in the Galactic Centre–A

1 Introduction

The Galactic Centre (GC) is an overwhelmingly crowded environment: it hosts a large population of late-type stars, few hundreds of early-type stars, several gas structures, and a $\sim{}4\times 10^{6}\,\mathrm{M}_{\sun{}}$ super-massive black hole (SMBH), associated with the radio source SgrA∗ (see, e.g. Mapelli & Gualandris 2016 for a recent review).

Being the closest galactic nucleus to us, the GC has been one of the most studied places in the Universe over the last decades, but still poses a plethora of questions. First of all, star formation is not expected to take place ‘in normal conditions’ within the central few parsecs of the Milky Way, because molecular clouds are disrupted by the extreme gravitational field of the SMBH (Sanders, 1998). Nevertheless, a population of $\sim{}30$ young ($\sim{}20-100$ Myr) B-type stars (the S-stars) are found within 0.04 pc from SgrA∗ (Morris, 1993; Ghez et al., 2003). Moreover, several hundred young ($\sim{}2-6$ Myr, Lu et al. 2013) early-type stars lie in the innermost parsec. About $20\%$ of them belong to the so-called clockwise (CW) disc, a thin-disc structure extended between $\sim 0.04$ and 0.13 pc from SgrA∗. CW disc members are found to follow a top-heavy (TH) mass function, best fitted by a power law with index $\alpha{}\sim{}1.7$ (Lu et al., 2013).

Many scenarios have been proposed to explain the formation of the early-type stars, and this issue is still debated. According to the migration scenarios, the young stars formed out of the central parsec, and then sank to the GC by dynamical mechanisms (e.g. Kim & Morris 2003; Kim et al. 2004; Fujii et al. 2008, 2010; Perets et al. 2009; Perets & Gualandris 2010, but see Petts & Gualandris 2017 for the most recent discussion of the issues connected with this scenario), while in the in situ scenarios star formation can be triggered by the disruption of a giant molecular cloud (e.g. Sanders 1998; Bonnell & Rice 2008; Hobbs & Nayakshin 2009; Alig et al. 2011; Mapelli et al. 2012; Alig et al. 2013; Mapelli et al. 2013b).

Moreover, two faint dusty objects, named G1 and G2, have been observed in the vicinity of the SMBH, on highly eccentric orbits (Clénet et al., 2004a, b; Clénet et al., 2005; Ghez et al., 2005; Gillessen et al., 2012, 2013; Pfuhl et al., 2015). In particular, G2 approached the SMBH to $\sim{}200$ AU in Spring 2014, during periapsis passage, and was not completely disrupted by the SMBH (Witzel et al., 2014). G1 and G2 have been extensively studied since their first discovery: they have been proposed to be gas clouds (e.g. Schartmann et al. 2012), dust-enshrouded low-mass stars (e.g. Burkert et al. 2012; Ballone et al. 2013; Scoville & Burkert 2013; Prodan et al. 2015), tidally-disrupted stars (e.g. Guillochon et al. 2014), or planetary embryos and/or proto-planetary discs (e.g. Murray-Clay & Loeb 2012; Mapelli & Ripamonti 2015; Trani et al. 2016a), but their origin remains unclear.

Another major issue about the GC is the observed cored profile of late-type ($>1$ Gyr) red giant stars (Buchholz et al., 2009; Do et al., 2009, 2013; Aharon & Perets, 2015). If the old stellar population in the GC were relaxed, we would expect it to follow a cusp profile with slope $\sim{}-7/4$ (Bahcall & Wolf, 1976), but red giant stars exhibit a cored profile, thus hinting to a possible missing cusp problem. A genuine cored distribution might be the result of the infall of star clusters (Antonini et al., 2012; Perets & Mastrobuono-Battisti, 2014), or of a recent SMBH binary merger (Merritt, 2010). Alternatively, red giant stars might not be representative of the late-type stellar distribution, as suggested in a series of recent studies considering the population of faint giants and sub-giants (Gallego-Cano et al., 2017; Schödel et al., 2017; Baumgardt et al., 2017): if only these populations are considered, the stellar cusp seem to be present, even if it appears to be shallower than expected. The still apparent lack of bright red giants might be explained as the outcome of stellar collisions, removing mainly red giants (Genzel et al., 1996; Alexander, 1999; Bailey & Davies, 1999; Dale et al., 2009).

Pulsars are also missing in the GC. Based on the large number of massive stars, we expect a large population of neutron stars (NSs) in the central parsec (up to few thousands in the innermost $\sim 10^{-2}$ pc, according to Pfahl & Loeb 2004), but only one was observed: PSR J1745-2900, a magnetar orbiting SgrA∗ with $\lesssim{}0.1$ pc of projected separation (Kennea et al. 2013; Mori et al. 2013; Rea et al. 2013, but see also Wang et al. 2006; Muno et al. 2008; Ponti et al. 2016 for other possible candidates). The observation of the magnetar implies that the lack of other pulsar detections is an intrinsic dearth of pulsars in the central parsec, rather than an effect of strong interstellar scattering screen toward the GC. Moreover, the magnetar formation efficiency might be higher in the GC than in other places (Dexter & O’Leary, 2014).

In this paper, we investigate the missing cusp problem and the missing pulsar problem by means of a simple argument: supernovae (SNe) occurring in binary systems in the central parsec might reshuffle the orbits of pulsars, other dark remnants, and old red giant stars. Moreover, SN kicks in binary systems might push low-mass stars into highly eccentric orbits, such as the ones of G1 and G2. SN kicks might be very large (Hobbs et al., 2005), and strongly affect the orbit of a binary star, in many cases breaking up the binary.

Thus, we use three-body simulations to study the dynamical effects of SNe occurring in binary stellar systems within a disc structure similar to the CW disc. We show that SN kicks reshuffle the orbit of the binary and contribute to sculpt the density profile of the GC.

2 Methods

To study the evolution of GC binaries after SN explosions, we run 30,000 three-body encounters of systems composed of a SMBH and a stellar binary. The mass of the SMBH is fixed to $4.3\times{}10^{6}\,{}M_{\odot}$ (Gillessen et al., 2009).

The binary centre of mass is initially assumed to sweep a Keplerian orbit around the SMBH, with eccentricity and semi-major axis in agreement with recent observations of the CW disc: the eccentricity $e$ is drawn from a Gaussian distribution centred in $\langle{}e\rangle{}=0.3\pm{}0.1$ (Yelda et al., 2014), while the semi-major axes $a$ follow the distribution $f(a)\propto a^{0.07}$ (Do et al., 2013), and we chose the mean anomaly of the orbit from a uniform distribution between 0 and $2\pi$; this ensures that the surface density profile of stars is $\Sigma(r)\propto r^{-0.93}$, consistent with the surface density distribution of the early type stars in the disc (Do et al., 2013); $a$ can vary in the range [$0.04,0.13$] pc or [$0.001,0.13$] pc depending on the simulation set, as detailed below.

The eccentricity $e_{\ast}$ of the stellar binary is randomly extracted from the thermal eccentricity distribution $f(e_{\ast}{}){\rm d}e_{\ast}{}=2e_{\ast}{}{\rm d}e_{\ast}{}$ (Jeans, 1919), while the semi-major axis $a_{\ast}{}$ distribution takes the form $f(a_{\ast}{}){\rm d}a_{\ast}{}\propto a_{\ast}{}^{-1}{\rm d}a_{\ast}{}$, i.e. follows the so-called Öpik’s law (Öpik, 1924). The semi-major axis of the binary system spans from a minimum of 0.2 AU to a maximum distance equal to $0.4\,{}r_{\rm j}$, where $r_{\rm j}=a\,{}(1-e)\,{}(\tfrac{1}{3}M_{b,i}/M_{\bullet})^{1/3}$ is the Jacobi radius (where $M_{b,i}=M_{1,i}+M_{2}$, $M_{1,i}$ is the mass of the primary before the SN explosion, $M_{2}$ is the mass of the secondary, and $M_{\bullet}$ is the mass of the SMBH)111 The resulting initial orbital period of the stellar binaries ranges from $\sim 0.1$ yr to $\sim 10^{3}$ yr, and about one fourth of them has a period $<1$ yr.. This ensures by construction that the binary star will not break by tidal forces at the closest passage with the SMBH, as shown in Trani et al. (2016b). We also ran a sub-sample of 3000 three-body simulations without SN kicks, to check that the binary does not break if SN kicks are not included in the simulations. We found that $a_{\ast}$ and $e_{\ast}$ change by less than $10^{-5}$ in the test simulations without SN kicks, indicating that our initial conditions are stable against tidal forces.

Since a large fraction of stars (especially massive stars) form in binary systems (e.g. Sana & Evans 2011), many CW disc members are expected to be binary stars. Gautam et al. (2016) found that the binary fraction of early-type main sequence stars within the CW disc is at least 32% (within the 90% confidence level). Stars with mass in excess of $\sim 9\,\mathrm{M}_{\sun{}}$ are thought to undergo SN explosions after about $3-50$ Myr, depending on their initial masses. Given the young age of stars in the CW disc ($\sim{}2.5-6$ Myr, Lu et al. 2013), most SN events will occur in the next Myrs. However, we use the orbital properties of the CW disc as initial conditions, because it is reasonable to expect that a series of CW-disc like star-forming events have occurred in the past and may occur in the future (it is unlikely that we live in a special time for the history of the GC, Loose et al. 1982).

The masses of the two stars are computed as follows: the zero age main sequence (ZAMS) mass of the more massive (primary) binary component ($M_{1,\rm ZAMS}$) is randomly distributed in the range $[9,150]\,\mathrm{M}_{\sun{}}$, sampling a TH mass function (${\rm d}N/{\rm d}m\propto m^{\alpha}$, with $\alpha=-1.7$), in agreement with the observations of the GC (Lu et al., 2013). The ZAMS mass of the less massive (secondary) binary component ($M_{2,\rm ZAMS}$) is either extracted from the TH mass function in the range $[0.1,\,{}150]\,\mathrm{M}_{\sun{}}$ (simulation sets A and C) or assumed to be equal to 0.9 times the ZAMS mass of the primary star (simulation set B), as detailed below. We then evolve each primary star from $M_{1,\rm ZAMS}$ to its pre-SN mass $M_{1,i}$, according to the PARSEC stellar evolutionary tracks (Bressan et al., 2012; Chen et al., 2015); the pre-SN mass of a star as a function of its ZAMS mass is shown in Figure 1. Similarly, we calculate the mass of the secondary star at the time of the SN explosion of the primary ($M_{2}$) using the PARSEC stellar evolutionary tracks for a star of ZAMS mass $M_{2,\rm ZAMS}$.

We start the simulations at the time of the SN explosion, by giving the binary a SN kick. The mass of the compact remnant is computed using the delayed SN model (Fryer et al., 2012), according to Spera et al. (2015). We distinguish between NSs and stellar-size black holes (BHs) depending on whether their mass is respectively below or above $3\,\mathrm{M}_{\sun{}}$. The mass of the SN relic as a function of the progenitor ZAMS mass is shown in Figure 1. The distribution of SN kicks is uncertain, especially for BHs (see e.g. Gualandris et al. 2005; Fragos et al. 2009; Repetto et al. 2012; Janka 2013; Beniamini & Piran 2016; Mandel 2016). For NSs, we adopt the kick distribution reported by Hobbs et al. (2005), i.e. a Maxwellian distribution with one-dimensional variance equal to $\sigma=265$ kms*-1*, derived on the basis of peculiar motions of 233 pulsars. If the remnant is a BH, the kick velocity is then normalized to the mass of the remnant, assuming linear momentum conservation. The natal kick is given to the binary members as detailed in Appendix A.

We performed three sets of three-body simulations composed of 10,000 runs each: (i) in set A, the ZAMS mass of the secondary star is extracted from the TH mass function between 0.1 and $150\,\mathrm{M}_{\sun{}}$ and the disc inner and outer limits are 0.04 and 0.13 pc, respectively; (ii) set B differs from set A as the ZAMS mass of the secondary is set equal to 0.9 times the ZAMS mass of the primary; we explored this possibility due to the uncertainty in the mass correlation of binary stars (Kobulnicky & Fryer, 2007; Sana et al., 2011); (iii) in set C the ZAMS mass of the secondary is computed as in set A, but the disc is assumed to span from 0.001 to 0.13 pc, according to Chen & Amaro-Seoane (2014), who suggested that the CW disc inner edge was closer to SgrA∗ in the past. The differences among the three sets and the parameters of the simulations are summarized in Table 1. We chose to perform 104 runs per each set of simulations because each set aims at reproducing the population of the CW disc, which has been estimated to be as massive as a few$\times 10^{4}\,\mathrm{M}_{\sun{}}$ (Lu et al., 2013). Moreover, we checked that 104 runs per each simulation set are sufficient to filter out statistical fluctuations.

Each three-body system is evolved from the moment of the SN explosion and for $\Delta t_{\rm sim}$=1 Myr; this timescale ensures that stars still bound to the SMBH generally complete $\sim{}100-1,000$ orbits by the end of the simulation, thus they totally absorb the effect of the kick. On the other hand, our simulations cannot account for the two-body relaxation, which can affect stellar orbits of a CW disc structure within $10-100$ Myr (Šubr & Haas, 2014); thus a choice of $\Delta t_{\rm sim}>1$ Myr would not be beneficial.

The simulations are evolved by means of a fully regularized N-body code that implements the Mikkola’s algorithmic regularization (MAR, Mikkola & Tanikawa 1999b, a). This code is designed for studying the dynamical evolution of few-body systems in which strong gravitational encounters are frequent and the mass ratio between the interacting objects is large. The MAR scheme removes the singularity of the two-body gravitational potential for $r\rightarrow{}0$, by means of a transformation of the time coordinate (see Mikkola & Tanikawa 1999b for details). Our implementation uses a leapfrog scheme in combination with the Bulirsh-Stoer extrapolation algorithm. The code integrates the equations of motion employing relative coordinates by means of the so called chain structure. This change of coordinates reduces round-off errors significantly (Aarseth, 2003): the relative error in the energy conservation of our simulations is always below $10^{-10}$. The code can be used as a stand-alone regularized N-body code or as a module of the direct N-body code HiGPUs-R (Spera, in preparation; see Capuzzo-Dolcetta et al. (2013) for the current version of HiGPUs).

Figure 2 shows an example of the binary orbit before and after the SN kick in set C: two stars of 9 and 0.1 $\,\mathrm{M}_{\sun{}}$ orbit the SMBH with semimajor axis $a\approx 45.0$ mpc and eccentricity $e\approx 0.34$. In this particular case the binary breaks, but both members remain bound to the SMBH. The secondary star settles on a different orbit with respect to the initial one (with semimajor axis $\approx 53.0$ mpc, eccentricity $\approx 0.60$ and change in the orbital inclination of $\approx 32$ degrees) while the 1.3 $\,\mathrm{M}_{\sun{}}$ NS produced by the SN only experiences a small ($\lesssim 10\%$) change in its orbital parameters.

3 Results

In this Section, we present the results of our three-body simulations. Figure 3 shows the density profile of the simulated binaries reconstructed from the individual three-body experiments, by considering the distance of each star from the SMBH at time $t=1$ Myr (i.e. the end of the simulations). Figure 4 shows the orbital parameter distribution of objects that are still bound to the SMBH at the end of the simulation. The semi-major axes, eccentricities and pericentre radii are computed with respect to the SMBH, while the inclination is computed with respect to the CW disc. In the plots, we distinguish between BHs ($\sim{}54\%$ of the remnants), NSs ($\sim{}46\%$ of the remnants) and secondary stars. We stress that the large number of BHs is due to the assumed TH mass function. The distribution of BHs resembles the pre-SN distribution, implying that most BHs are anchored to their initial orbit and SN kicks are not efficient in scattering them in different regions of the phase space. In contrast, NSs are more easily scattered by the SN explosion and spread in the innermost and outermost regions of the GC. This is valid for all the sets of simulations.

In both sets A and C, $75-80\%$ of the mass within $0.1-1$ pc is in the form of BHs. From Figure 3 it is also clear that in set A and C BHs dominate the final mass distribution for a large range of radii, i.e. within $\sim 0.02-0.2$ pc in set A and $\sim 5\times 10^{-4}-0.2$ in set C. We stress that our simulations were run only for 1 Myr and that they only account for three bodies at a time: we cannot say anything about the effect of two-body relaxation on the BHs confined in the central pc. However, it is likely that the BHs will relax their distribution near the SMBH within $\sim 100$ Myr; the relaxation time within the disc is computed as described in Šubr & Haas (2014)222 As a matter of fact, the relaxation time-scale of the old and almost isotropic stellar component within the SMBH influence radius has been estimated to be $\sim 10$ Gyr or even longer (Merritt, 2010, fig. 3). However our paper focuses on the young population of stars in the CW disc: due to their coherent motion, stars in a disc self-relax on a much shorter time-scale $T_{\langle{}e^{2}\rangle{}}$, that can be computed as the time over which the mean-square-eccentricities of stars grow from zero to a given $\langle{}e^{2}\rangle{}$ (Stewart & Ida, 2000). If we assume a CW-disc like structure and we use ${\langle{}e^{2}\rangle{}}^{1/2}\approx 0.3$, it follows that $T_{\langle{}e^{2}\rangle{}}\sim 100$ Myr (Šubr & Haas, 2014, equation 4). .

The behaviour of the secondary stars strongly depends on their masses: in set A and C, the secondary binary member is a low-mass star in most cases ($\sim 90\%$ secondary stars have masses lower than $3\,\mathrm{M}_{\sun{}}$) and is easily scattered on a completely new orbit after the SN kick, while in set B the secondary is a high-mass star ($\gtrsim 9\,\mathrm{M}_{\sun{}}$) and keeps memory of its initial orbit. In set B, secondary stars are much more massive than the remnants of primary stars, thus they dominate the mass distribution within $\sim 0.02-0.5$ pc. This would be true only within a short timescale, as most of the companion stars in set B will undergo a SN explosion within a few Myr and they will turn into compact remnants.

Figure 4 shows that NSs notably change their orbit in set B, with respect to set A and C. We detail the case of eccentricity: in set A and C only $\sim{}16\%$ of NSs have eccentricities $e>0.5$, while in set B the fraction of NS with $e>0.5$ is as high as $\sim{}55\%$. On the other hand, $\sim{}25\%$ of secondary stars have eccentricities $e>0.5$ in set A and C, while their fraction drops to $\sim{}5\%$ in set B. All the other orbital parameters show a similar behaviour, and show that SN kicks are distributed between the two binary components, with weights that anti-correlate with the mass of the bodies.

Figure 5 shows stellar pericentre distance versus eccentricity for the three sets of runs, distinguishing between compact remnants and secondary stars. This Figure highlights that heavy objects are spatially more confined. The bulk of stars and remnants still show orbital parameters that resemble the ones of the CW disc; however, several stars are scattered to highly eccentric orbits.

Most of stars and compact remnants continue orbiting the SMBH by the end of the simulation, and only $\sim{}5\%$ objects become unbound from the SMBH. In set B, NSs absorb most of the SN kick and about one fourth of them escapes the SMBH potential, while secondary stars are strongly anchored to their initial orbit and all but five of them stay on a bound orbit by the end of the simulation. We find an opposite trend in sets A and C, where secondary binary members are generally lighter than NSs: about 10% companion stars end up unbound, while only a few per cent of NSs escape. Only few BHs become unbound from the SMBH: 0%, 0.4%, and 0.04% in sets A, B and C, respectively. Table 2 is a summary of these possible outcomes.

Figure 6 shows the cumulative distribution of radial velocities of stars and compact remnants that escape the SMBH potential well. One to a few per cent of these objects attain escape velocities $>900$ km s*-1*, thus they may potentially be detected as hyper-velocity stars escaping the galactic potential (a more accurate investigation was already done by Zubovas et al. 2013). In set A and C, hyper-velocity stars are secondary binary members (except for one NS in set A). In set B only five secondary stars escape from the SMBH, while the NSs are the fastest escapers. In four cases, both binary members escape the SMBH potential and they stay bound in a binary system.

Stellar binaries break after the SN explosion in the $50-60\%$ of runs. Table 3 lists the fraction of binaries that survive the SN kick. A higher fraction of binaries do not break in set B, because secondary stars are more massive than compact remnants; however secondary stars will undergo SN explosion in a short timescale, thus even these binaries might be disrupted in few Myr.

4 Discussion

4.1 The population of dark remnants

Our results show that in all simulation sets SN explosions are not effective in scattering massive objects on new orbits. If we assume continuous star formation in the inner pc, our simulations predict that a dark cluster of BHs may have formed within the influence radius of the SMBH, because natal kicks do not modify the BH orbital properties significantly.

In set B, secondary stars (which are particularly massive) behave in a similar way as BHs, dominating the mass distribution within $\sim 0.02-0.5$ pc. However, almost all these high-mass stars will undergo a SN explosion within a few Myr and will get a further kick, changing their density distribution around the SMBH. The study of their subsequent evolution is beyond the aim of this paper; however we can speculate about their possible final distribution. It is reasonable to expect that a considerable fraction of the heavier BHs generated from secondary stars will keep moving in the same region, as SN kicks are expected to be weaker if the remnant is a BH. As a consequence, we expect the BH population of set B to acquire further members after the SN explosion of secondary stars. Thus, our simulations clearly show that an excess of dark remnants may hide inside the central parsecs of the Galaxy.

The angular distribution of objects resembles the initial one, i.e. all stars and remnants are preferentially found in the plane of the CW disc, even if light stars spread more their angular distribution, as shown in Figure 7; thus the final distribution of stars and compact remnants is flattened in our simulations. However, several star formation episodes might have occurred in disc-like structures (such as the CW disc) around the SMBH, with possibly different orientations. This might have contributed to build a cluster of dark remnants along the history of the GC. Several episodes of star formation within randomly oriented discs would result in a relatively isotropic distribution of remnants; in addition, secular processes as Kozai-Lidov oscillations (Kozai, 1962; Lidov, 1962) or resonant relaxation (Rauch & Tremaine, 1996) may have assisted the isotropisation of orbits.

The presence of a cusp within the inner parsec of our Galaxy can also give new hints on the origin of S-stars: Hills (1988) proposed that both the S-star cluster and hyper-velocity stars within the galactic bulge may have originated by the tidal breakup of binaries passing close to the SMBH; Šubr & Haas (2016) investigated this possibility via $N$-body simulations, and they reproduce fairly well the orbital parameters of S-stars when they take in account the presence of a possible stellar or dark cusp.

4.2 The NS distribution

In our simulations, NSs are easily scattered far away from the SMBH, especially when they are bound to a massive companion. SN kicks are believed to have isotropic orientations; our simulations show that NSs have a greater probability to end up further away from the SMBH rather than closer in, as expected from geometrical arguments333 If we take a sphere $S$ of radius $r$, centered on the SMBH and passing through the position $P$ of the exploding SN, and we consider the family of straight lines through $P$ (i.e. all possible directions of the SN kick), half of them is found to avoid the sphere $S$; the other half will enter $S$, but it has to exit the sphere after covering a distance in the range ($0,2r$). This means that if we consider all the possible directions of the kick, in most cases the NS will end up farther away from the SMBH.. NS natal kicks may attain high velocities according to Hobbs et al. (2005); in our simulations, we implicitly assumed that SN kicks have a constant distribution in momenta rather than in velocities, thus it is not surprising that NSs are heavily affected by their natal kick. We also expect that SN kicks will deeply affect NS orbits, even when they are not part of a binary system. Thus, we suggest that the dearth of NSs observed in the GC may be due to SN kicks.

Assuming a star formation rate in the GC of $\sim 10^{-3}$ stars yr*-1* (Pfuhl et al., 2011) and introducing a correction factor to account for the fact that we did not simulate the low-mass end of the initial mass function (see Section 2), we can compute the number of NSs that may remain in the GC within 10 Myr or 1 Gyr, i.e. the typical lifespan of respectively a canonical pulsar (CP) and a millisecond pulsar (MSP). Considering only sets A and C, we find that $\sim 500$ ($\sim 800$) NSs would be located within the innermost 0.1 pc (1 pc) after 10 Myr, and $\sim 5\times 10^{4}$ ($\sim 8\times 10^{4}$) would be found in the same region after 1 Gyr.444 Note that the estimates we present are to be considered upper bounds on the number of NS that would remain inside 0.1 or 1 pc; as a matter of fact, we neglected NSs without a companion by assuming the binary fraction to be 1, but isolated NSs cannot share their natal kick with a companion and we expect them to spread farther away from the SMBH even more than NSs in binaries.

Rajwade et al. (2016) recently estimated the observable fraction of NSs in the inner pc of the Galaxy, and found that only the $\sim$2% of CPs and $0.1\%$ of MSPs are nowadays detectable, even neglecting free-free absorption and multi-path scattering. They conclude that the maximum number of pulsars pointing toward us compatible with the null result of present-day surveys is 50 CPs and 1430 MSPs, suggesting that we still do not have a definitive proof of missing pulsars. Following the calculations of Rajwade et al. (2016) and the results in the survey by Macquart & Kanekar (2015), if we assume that all our simulated NSs behave as CPs (MSPs), we would expect to observe $\sim$5 ($\sim 250$) of them at most; in addition, NSs can behave as MSP only in case they stay bound to the companion star or they dynamically acquire a massive companion; only 2% of NSs in sets A and C do not break after the SN explosion, thus the number of candidate MSP in our sample could be reduced by a factor $\sim 10^{2}$. We will investigate these aspects in more detail in a forthcoming paper.

4.3 The S-cluster, G1 and G2

The third column of Figure 5 and the three panels in Figure 8 and 9 compare the orbital elements of secondary stars in our simulations with those of the S-stars, G1 and G2. The scatter plots show the pericentre radii ($r_{\rm p}$), semi-major axes ($a$) and orbital inclinations ($i$) versus the orbital eccentricity ($e$), respectively. From these Figures it is apparent that the orbits of the S-stars with $e\lesssim{}0.5$ may originate from the SN kick of a massive companion, especially in Set C, where the inner edge of the disc is closer to the SMBH. However, most S-stars with eccentricities $\lesssim 0.5$ are compatible with a CW-disc origin even without invoking a SN event.

A second bunch of S-stars have eccentricity $e\gtrsim{}0.7$; the orbits of these S-stars (as well as those of G1 and G2) are marginally consistent with our distributions in sets A and C, but not in set B. In particular, the 2.5% (3%) of secondary stars have $r_{\rm p}<10^{-2}$ pc and $e>0.7$ in our simulation set A (set C).

The orbital properties of a small fraction of the simulated stars in set A and C match those of G1 and especially G2. The fraction of secondary stars with eccentricity and semi-major axis compatible with those of G2 is 0.8% and 1.0% in set A and C, respectively. The fraction of secondary stars with eccentricity and semi-major axis compatible with those of G1 is equal to 0.15% for set C, while no stars match G1’s orbit in set A. Low-mass stars experience a long pre-main sequence phase, and by the time they exit this phase a massive star born in the same star formation episode may experience the SN explosion; we propose that G1 and G2 are pre-main sequence stars that were orbiting a massive star in the CW-disc, and they were kicked on a highly eccentric orbit by the SN kick of the companion.

4.4 X-ray binaries

Chandra X-ray observations of the GC and Galactic bulge show that X-ray emitters are highly concentrated in the GC, in excess with respect to the expectations from stellar distribution models (Muno et al., 2005; Hong et al., 2009). This excess is mainly attributed to magnetic cataclysmic variables and primarily intermediate polars (Hong et al., 2009). At least part of this X-ray excess might originate from an accreting BH or NS binary.

In our simulations, a significant fraction ($40-50\%$) of stars is still part of a binary after the SN kick (Table 3). To provide an order-of-magnitude estimate of how many of these simulated binaries might become X-ray emitters, we consider the minimum (min($R_{\ast}$)) and maximum (max($R_{\ast}$)) radius of each secondary star during its life, as obtained from the PARSEC stellar evolutionary tracks (Bressan et al., 2012). We then compare the minimum and maximum radii with the Roche lobe limit, obtaining a lower and upper limit for the number of binaries that may undergo Roche lobe overflow in our simulations. Our calculations show that $0-30\%$ of surviving binaries may undergo a phase of Roche-lobe overflow in sets A and C. The fraction of estimated Roche-lobe overflow systems is $0-76$% in set B: if the secondary star is almost as massive as the primary, a non-negligible fraction of binaries in our simulation may become X-ray binaries (but we stress that our back-of-the-envelope calculation does not account for the lifetime of the companion star). Table 4 shows these results.

Hong et al. (2009) found that the density of X-ray sources within a region of $40\times 40$ pc2 surrounding SgrA* is $\sim{}10^{-6}\,\mathrm{M}_{\sun{}}^{-1}$, compatible with our generous upper limits. However, Hong et al. (2009) take into account all observed X-ray sources within an area much larger than the central pc, thus a direct comparison with our results cannot be done. Moreover, to compare our results with observations we would need to account for binary evolution processes (such as Roche lobe overflow and wind accretion), which are not included in our simulations. A more accurate study of this aspect will be the subject of a forthcoming paper.

4.5 Caveats

In our simulations, we did not take in account any secular processes. Mass precession may occur due to the spherical distribution of late-type stars in the GC; however this process would only cause pericentre advance and would let our results unaffected (Šubr & Haas, 2012). Kozai-Lidov resonances (Kozai, 1962; Lidov, 1962) may have a role especially for stars with large semi-major axes; however Kozai cycles are expected to occur for stars misaligned to the CW disc; in addition, the spherical distribution of old stars would damp the cycles even for stars with inclined orbits. The net effect we may expect is a change in the longitude of the ascending node and argument of the pericentre, but again our results would probably be unchanged (Mapelli et al., 2013b; Trani et al., 2016b). Moreover, two-body encounters can significantly affect the distribution of eccentricities and semi-major axes (e.g. Trani et al. 2016b).

A second possible issue of our simulations is the dependence of our results on the adopted kick distribution. The SN kick distribution is still highly debated (e.g. Beniamini & Piran 2016, and references therein). We aim at investigating this issue in a future paper. However, we notice that our main result (i.e. the formation of a dark cusp) remains mostly unaffected by the prescription we adopt for the kick. In fact, recent studies (Fryer et al., 2012; Mapelli et al., 2013a) suggest that BH kicks might be close to zero for all direct-collapse BHs, while our model attributes quite large kicks to BHs, even if they come from direct collapse. In this respect, our results represent a lower-limit for the building-up of a dark cusp.

5 Summary

In this paper, we investigate the possibility that SN kicks occurring in binary systems affect the orbits of stars and stellar remnants in the GC, by means of regularized three-body simulations. We performed three main sets of simulations, changing the mass distribution of the secondary component of the binary and the inner edge of the initial disc.

Our results show that BH remnants do not change their orbits significantly. Thus, several episodes of star formations in the GC may have built-up a cusp of BHs within the radius of influence of the SMBH. In contrast, NSs are generally scattered away from SgrA∗ and they may end up sweeping highly eccentric orbits; thus the missing pulsars problem may be a result of SN natal kicks, especially if a great number of NS progenitors are initially bound to a companion with a similar mass. Finally, we propose that G1 and G2 might have been low-mass stellar companions to massive CW-disc stars. When their massive companions underwent a SN explosion, G1 and G2 attained their highly eccentric orbits around the SMBH. This scenario deserves further investigation in a forthcoming paper.

Acknowledgements

We warmly thank the anonymous referee for his/her useful comments and suggestions. We also thank Federico Abbate, James Petts, Alessia Gualandris, Maureen van den Berg, Melvyn B. Davies, Mark Gieles, John Magorrian, Alessandro Trani and Fabio Antonini for useful discussions and suggestions. We acknowledge the CINECA Awards N. HP10CP8A4R and HP10C8653N, 2016 for the availability of high performance computing resources and support. We acknowledge financial support from the Istituto Nazionale di Astrofisica (INAF) through a Cycle 31st PhD grant, from the Italian Ministry of Education, University and Research (MIUR) through grant FIRB 2012 RBFR12PM1F, from INAF through grant PRIN-2014-14, from the MERAC Foundation and from the Fondazione Ing. Aldo Gini.

Appendix A Prescriptions for the natal kick on the binary members

Here we describe the procedure we adopted to distribute the SN kick to the two members of the stellar binary system, according to the prescriptions detailed in appendix A1 by Hurley et al. (2002).

If we consider a binary stellar system composed of two stars with initial masses $M_{1,i}>M_{2}$ and relative velocity vector $\bf r$, we can express their relative velocity as

[TABLE]

where $\beta$ is the angle between the position and velocity vectors; we use the notation $\hat{\textbf{\i}},\hat{\textbf{\j}},\hat{\bf k}$ to denote the unit vectors relative to a reference frame $X,Y,Z$ where the mass $M_{2}$ is in the origin and the specific angular momentum ${\mathbf{h}}$ of the binary lies along the positive $Z$ axis. The system geometry is shown in Figure A1 of Hurley et al. (2002). The angle $\beta$ can be expressed in terms of the orbital elements of the binary, i.e. the semimajor axis $a$, eccentricity $e$, eccentric anomaly E and the distance $r$ between the stars:

[TABLE]

We chose the distance at which the SN explosion occurs by drawing a mean anomaly $\mathcal{M}$ (that is uniformly distributed in time) from a uniform distribution between 0 and 2$\pi$. For a given $\mathcal{M}$, we compute the associated eccentric anomaly by recursively solving

[TABLE]

where Ei denotes the $i$-th iteration for the computation of E; the tolerance for the convergence of E was set to $10^{-9}$. Then the distance $r$ between the stars can be computed as:

[TABLE]

The orbital speed $v_{0}$ can as well be expressed as a function of the orbital elements:

[TABLE]

where $M_{b,i}$ denotes the total mass of the binary before the SN explosion.

We assume that the mass of the envelope, $\Delta M$, is instantaneously ejected by the primary star during its explosion, and its mass becomes $M_{1,f}=M_{1,i}-\Delta M$; thus the total mass of the binary after the SN is $M_{b,f}=M_{b,i}-\Delta M$. During the explosion, the distance between the stars does not vary, but the primary star experiences a velocity kick whose modulus $v_{k}$ is randomly drawn from a Maxwellian with one-dimensional variance equal to 265 km s*-1* (Hobbs et al., 2005); if the remnant is a BH (i.e. $M_{1,f}>3\,\mathrm{M}_{\sun{}}$), the kick velocity is then normalized to the mass of the remnant, assuming linear momentum conservation. The kick velocity vector can be expressed as

[TABLE]

the direction of the kick is assumed to be isotropic in space, thus $\phi$ is distributed between $-\pi/2$ and $\pi/2$ as $P(\phi)\propto\cos\phi$, while $\omega$ is uniformly distributed between 0 and 2$\pi$. After the SN explosion took place, the new velocity between the stars is:

[TABLE]

and its modulus is, according to equation (6),

[TABLE]

Note that we use the subscript $f$ to denote the new binary orbital parameters after the SN kick. The previous equation can be solved to obtain the new semimajor axis $a_{f}$ of the binary system. The system specific angular momentum after the SN is equal to

[TABLE]

and its squared modulus is

[TABLE]

where

[TABLE]

the previous equation can be solved to compute the new value for the eccentricity $e_{f}$. The new angular momentum vector is shifted with respect to the initial one by an angle $\eta$ such that

[TABLE]

Since an amount of mass $\Delta M$ is ejected from the primary star, the centre of mass of the system now exhibit a velocity shift equal to

[TABLE]

with respect to the initial reference frame.