Gravitational-wave Signatures from Compact Object Binaries in the Galactic Center

The recent detection of gravitational-wave (GW) emission from merging stellar-mass black hole (BH) binaries and neutron star (NS) binaries has expanded our ability to understand the universe (e.g., Abbott et al. 2016a, 2016b, 2016c, 2016d, 2017a, 2017b, 2017c, 2017d, 2019a; LIGO Scientific Collaboration et al. 2020). However, it remains challenging to distinguish between the formation channels of these merging sources. Recent studies have emphasized the significant contribution of dynamical formation channels in dense stellar environments to the overall GW signals (e.g., Portegies Zwart & McMillan 2000; Miller & Lauburg 2009; O'Leary et al. 2009; Banerjee et al. 2010; Downing et al. 2011; Antonini & Perets 2012; Rodriguez et al. 2015; Cholis et al. 2016; Gondán et al. 2018; Hoang et al. 2018; Lower et al. 2018; Randall & Xianyu 2018; Rodriguez et al. 2018; Samsing & D'Orazio 2018; Zevin et al. 2019; Hoang et al. 2020). One of the unique places that contributes to this overall GW signal via dynamical formation is the center of galaxies (e.g., Antonini & Perets 2012; Petrovich & Antonini 2017; Rodriguez et al. 2016a, 2016b; Hoang et al. 2018; Stephan et al. 2019).

Almost every galaxy has a supermassive BH (SMBH) at its center (e.g., Kormendy & Richstone 1995; Ferrarese & Ford 2005; Kormendy & Ho 2013). The closest SMBH to Earth is Sagittarius A*, with 4 × 10⁶ M_⊙ at the center of the Milky Way (e.g., Ghez et al. 2005; Gillessen et al. 2009). Hence, detailed observations of the Galactic Center (GC) can provide valuable insights into the environment and dynamics that must exist in other galaxies as well.

Surrounding the SMBH at the center of our Galaxy is a dense environment called the nuclear star cluster, which includes populations of mostly old stars (≳1 Gyr; e.g., Lu et al. 2009; Bartko et al. 2010; Do et al. 2013a, 2013b; Feldmeier-Krause et al. 2015; Nogueras-Lara et al. 2019; Schödel et al. 2020), as well as a subset of young stars (4–6 Myr; e.g., Schödel et al. 2003; Ghez et al. 2005, 2008; Gillessen et al. 2009, 2017). Binaries that may exist in this nuclear star cluster undergo frequent interactions with neighboring stars, as well as gravitational perturbation from the SMBH (e.g., Heggie 1975; Hills 1975; Heggie & Hut 1993; Rasio & Heggie 1995; Heggie & Rasio 1996; Binney & Tremaine 2008; Hopman 2009; Stephan et al. 2016; Hoang et al. 2018; Hamers & Samsing 2019; Rose et al. 2020). These interactions result in astrophysical phenomena such as hypervelocity stars (e.g., Hills 1988; Yu & Tremaine 2003; Ginsburg & Loeb 2007) and stellar and compact binary mergers (e.g., Antonini et al. 2010, 2011; Prodan et al. 2015; Stephan et al. 2016, 2019; Hoang et al. 2018). Of course, all these phenomena require high abundance of binaries in the vicinity of the SMBH, and the survival of such binaries depends on the stellar number density (e.g., Alexander & Pfuhl 2014; Rose et al. 2020).

Binaries are common in our Galaxy, with more than half of KGF stars and more than 70% of OBA stars having a stellar companion (e.g., Raghavan et al. 2010). Thus, it is reasonable to assume that binaries are also common at the GC. Already, there are three confirmed binaries within ∼0.2 pc of the GC observed by spectroscopy. The first case is IRS 16SW, a massive eclipsing binary examined by Ott et al. (1999) and Martins et al. (2006) with ∼50 M_⊙ for each object. Additionally, Pfuhl et al. (2014) reported two binary systems. One is a long-period binary with mass components >30 M_⊙, and the other is an eclipsing Wolf–Rayet binary of ∼20 and ∼10 M_⊙. In particular, Rafelski et al. (2007) proposed that the total mass fraction of massive binaries in the GC is comparable to that of the Galaxy's O stars and is about 7% of the total massive stellar population. Additionally, Stephan et al. (2016) suggested that around 70% of the initial binary population in the GC is expected to remain from the last star formation episode, which occurred 6 Myr in the past.

Moreover, the abundant X-ray sources detected within the GC indicate a potential stellar companion feeding accreting BHs (i.e., X-ray binaries; e.g., Muno et al. 2005; Cheng et al. 2018; Zhu et al. 2018; Mori et al. 2019). On the other hand, Muno et al. (2006, 2009) and Heinke et al. (2008) suggested that these X-ray sources could instead be cataclysmic variables. Additionally, the recent discovery of gas-like objects, the first of which was G2 (Gillessen et al. 2012), suggested the high potential existence of young binaries in the GC (e.g., Witzel et al. 2014, 2017; Ciurlo et al. 2020). Lastly, recent work by Naoz et al. (2018) showed that some puzzling properties of the stellar disks could be explained with the existence of binaries.

Within the vicinity of an SMBH, a stable binary has a tighter orbital configuration than the orbit of its center of mass around the SMBH. In such a system, gravitational perturbations from the SMBH can induce large eccentricities on the binary orbit, known as the "eccentric Kozai–Lidov" mechanism (EKL; e.g., Kozai 1962; Lidov 1962; Naoz 2016). However, we note that the dense environment of a nuclear star cluster surrounding the GC also provides a high chance of encounters. This may lead to the overall unbinding of binaries or lead to the capturing and hardening of compact object binaries (e.g., Heggie 1975; Heggie & Hut 1993; Heggie & Rasio 1996; Binney & Tremaine 2008; Hoang et al. 2020; Rose et al. 2020).

Recently, Stephan et al. (2019) investigated the dynamical evolution of binary stars in the vicinity of an SMBH subject to the EKL, including tidal interactions, general relativity (GR), and single and binary stellar evolution. They showed that while 75% of stellar binaries in the GC that are interacting with the central SMBH become unbound by the interaction, the remaining 25% will merge after a few Myr. Of the merging binaries, ∼14.6% will become compact objects binaries, while the remaining ∼85.4% will merge while the stellar components are main-sequence, red giant, or helium star phase. The final results of Stephan et al. (2019) show that 1.8% of those ∼14.6% compact objects will form a BH binary (BH–BH), 0.6% will form a BH–NS binary, 1.2% will become white dwarf–NS binaries (WD–NS), 15.1% will become WD binaries (WD–WD), and no NS binaries (NS–NS) can form.

We note that natal kicks that compact binary objects received during their supernova (SN) explosions also affect the number of binaries formed and their orbital parameters (e.g., Kalogera 2000; Bortolas et al. 2017). In the GC, where binaries undergo dynamical hierarchy interactions with the SMBH, natal kicks might eject those compact binaries with a high escape velocity to completely unbind the triple systems (e.g., Antonini & Rasio 2016; Michaely et al. 2016; Bortolas et al. 2017; Parker 2017; Lu & Naoz 2019). On the other hand, Lu & Naoz (2019) recently showed that SN kicks can more often result in shrinking the separation than expanding the orbit, thus contributing to possible GW events.

Currently, terrestrial GW detectors can only observe merging BH or NS binaries during their final inspiral phase. Nevertheless, detections of GW emission from binaries still in orbit are significant in revealing the binaries' formation history (e.g., Breivik et al. 2016; Nishizawa et al. 2017). Those sources can best be resolved via the Laser Interferometer Space Antenna (LISA), which is sensitive to mHz frequencies (e.g., Folkner et al. 1998; Amaro-Seoane et al. 2017; Robson et al. 2018, 2019). Observation of those GW sources can potentially contribute to our understanding of close binary evolution, the distribution of X-ray sources, SN explosions, gamma-ray bursts, and galactic structure (e.g., Yu & Jeffery 2010). Fortunately, binaries close to the SMBH exhibit measurable eccentricity oscillation owing to the EKL cycle and are thus distinctive from isolated field binaries (e.g., Hoang et al. 2019; Randall & Xianyu 2019; Deme et al. 2020; Emami & Loeb 2020).

The paper is organized as follows. We first provide the basic equations for estimating the signal-to-noise ratio (S/N; Section 2). Next, we analyze the potential detectability of compact objects within the vicinity of the GC while being agnostic to the formation mechanism (Section 3). Then, we generate a large population of binaries at the onset of their Roche limit crossing following their EKL evolution (Section 4.1). Next, we evolve these binaries using the COSMIC (Breivik et al. 2020a) code into compact object binaries (Section 4.2) and investigate the resulting GW signatures in terms of the LISA sensitivity curve (Section 4.3). We offer a crude approximation of the LIGO detection rate in Section 4.4 and a final discussion in Section 5.

For completeness we specify the relevant equations for calculating the LISA S/N; see Kocsis et al. 2012; Robson et al. 2019 and Hoang et al. (2019) for a complete derivation. The S/N of a binary, with a semimajor axis a and eccentricity e, is given by

$$\begin{eqnarray}&&\langle {\rm{S}}/{{\rm{N}}}^{2}(a,e)\rangle =\,\displaystyle \frac{16}{5}\int \displaystyle \frac{| \tilde{h}(a,e,f){| }^{2}}{{S}_{n}(f)}{df},\end{eqnarray} \tag{ 1 }$$

where, S_n (f) is the effective noise power spectral density,

$$\begin{eqnarray}&&{S}_{n}(f)=\displaystyle \frac{{P}_{n}(f)}{{ \mathcal R }(f)},\end{eqnarray} \tag{ 2 }$$

and has units of Hz⁻¹. ${ \mathcal R }(f)$ is the dimensionless sky- and polarization-averaged signal response function of the instrument, and P_n (f) is the power spectral density of the detector noise.

The square of the GW strain amplitude in the frequency domain, $| \tilde{h}(a,e,f){| }^{2}$, from Equation (1), is the sum of the strain amplitude at each orbital frequency harmonic, f_n . In other words,

$$\begin{eqnarray}&&| \tilde{h}(a,e,f){| }^{2}\approx \sum _{n=1}^{\infty }{h}_{n}^{2}(a,e,n){T}_{\mathrm{obs}}^{2}{\left(\displaystyle \frac{\sin (\pi ({f}_{n}-f){T}_{\mathrm{obs}})}{\pi ({f}_{n}-f){T}_{\mathrm{obs}}}\right)}^{2}\end{eqnarray} \tag{ 3 }$$

(e.g., Kocsis et al. 2012; Hoang et al. 2019), where

$$\begin{eqnarray}&&{h}_{n}(a,e,{f}_{n})=\displaystyle \frac{2}{n}\sqrt{g(n,e)}{h}_{0}(a),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{h}_{0}(a)=\sqrt{\displaystyle \frac{32}{5}}\displaystyle \frac{{G}^{2}}{{c}^{4}}\displaystyle \frac{{m}_{1}{m}_{2}}{{D}_{l}a},\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}g(n,e)=\displaystyle \frac{{n}^{4}}{32}\left[{\left({J}_{n-2}-2{{eJ}}_{n-1}+\displaystyle \frac{2}{n}{J}_{n}+2{{eJ}}_{n+1}-{J}_{n+2}\right)}^{2}\right.\\ \,\left.+\,(1-{e}^{2}){\left({J}_{n-2}-2{J}_{n}+{J}_{n+2}\right)}^{2}+\displaystyle \frac{4}{3{n}^{2}}{J}_{n}^{2}\right],\end{array}\end{eqnarray} \tag{ 6 }$$

where D_l is the luminosity distance between the source and detector, T_obs is the observation time of each binary source, G is Newton's gravitational constant, and J_i is the ith Bessel function evaluated at each n, e (Peters & Mathews 1963). The frequency harmonic, f_n in Equation (3) is defined as f_n = nf_orb, where

$$\begin{eqnarray}&&{f}_{\mathrm{orb}}(a)=\displaystyle \frac{1}{2\pi }\sqrt{\displaystyle \frac{G({m}_{1}+{m}_{2})}{{a}^{3}}}.\end{eqnarray} \tag{ 7 }$$

The FWHM of the integral from Equation (3) gives

$$\begin{eqnarray*}&&\int \mathrm{sinc}{\left(\pi ({f}_{n}-f){T}_{\mathrm{obs}}\right)}^{2}{df}\approx 0.885895/{T}_{\mathrm{obs}}.\end{eqnarray*}$$

Thus, combining Equations (1), (3), and (4), the final expression of the S/N of LISA can be approximate as

$$\begin{eqnarray}&&\langle {\rm{S}}/{\rm{N}}(a,e)\rangle \approx \displaystyle \frac{8}{\sqrt{5}}{h}_{0}(a)\sqrt{0.885895\times {T}_{\mathrm{obs}}\sum _{{n}_{\min }}^{{n}_{\max }}\displaystyle \frac{g(n,e)}{{n}^{2}{S}_{n}({f}_{n})}}.\end{eqnarray} \tag{ 8 }$$

In this paper, we calculate the S/N of each binary system based on the entire LISA mission lifetime T_obs = 4 yr, to accumulate the highest signals. However, for systems that merge within 4 yr, we take T_obs as their merging timescale.

Additionally, the characteristic strain of an evolving binary with eccentric orbits can be crudely approximated with the Fourier transform of a stationary binary of Equation (3) (Hoang et al. 2019):

$$\begin{eqnarray}&&{h}_{c}^{2}(a,e,f)=4{f}^{2}{\left|\tilde{h}(a,e,f)\right|}^{2}\times \min \left(1,\displaystyle \frac{{f}_{n}}{{\dot{f}}_{n}}\displaystyle \frac{1}{{\rm{\Delta }}{T}_{\mathrm{obs}}}\right),\end{eqnarray} \tag{ 9 }$$

Here, the factors within the minimum functions are taken as a result of the fact that each signal only accumulated within their own lifetime or the LISA observation timescale (e.g., Cutler & Flanagan 1994; Flanagan & Hughes 1998). If a and e change insignificantly over the observation timescale of LISA, the signal power only accumulates in each frequency bin. This equation is used as an approximation of the binary's strain curve in the LISA parameter space.

While EKL is one of the processes that can induce compact object binaries, other mechanisms can also drive the formation of compact object binaries. These processes include, but are not limited to, captures (e.g., O'Leary et al. 2009; Hoang et al. 2020), hardening by weak interactions (e.g., Heggie 1975; Alexander & Pfuhl 2014; Rose et al. 2020), three-body interactions (e.g., Binney & Tremaine 2008), etc. Thus, as a first step, we offer a proof of concept that is agnostic to the binary formation mechanism.

We investigate six representative examples of different compact objects involving BHs, NSs, and WDs, as listed in Table 1. Those examples aim to provide a proof of concept of the S/N of the potential compact object binaries within the GC viewing via the LISA detector.

In Figure 1 we depict the S/N for the six representative examples using Equation (8), considering a wide range of separation and eccentricity. ⁶ As expected, the S/Ns of WD binaries are low (≲100, but still larger than 5) for a wide range of the parameter space. The S/Ns for NS binaries are mostly above ∼100, while the S/Ns for binaries with BHs can be as high as 10⁵, with the majority of the parameter space yielding S/Ns > 100.

Zoom In Zoom Out Reset image size

Those compact object binaries will eventually merge via the timescale estimated as

$$\begin{eqnarray}&&\begin{array}{l}{t}_{\mathrm{GW}}=1.602\times {10}^{5}\mathrm{yr}\times {\left(\displaystyle \frac{{m}_{1}}{10{\text{}}{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{m}_{2}}{10{\text{}}{M}_{\odot }}\right)}^{-1}\\ \,\times \,{\left(\displaystyle \frac{{m}_{1}+{m}_{2}}{20{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{a}{0.01\mathrm{au}}\right)}^{4}\times f(e){\left(1-{e}^{2}\right)}^{7/2}\end{array}\end{eqnarray} \tag{ 10 }$$

(e.g., Blaes et al. 2002). In Figure 1, we overplot three representative merger times of 10, 10⁴, and 10⁷ yr.

As represented in Figure 1, the BH–BH systems (top row) with the larger S/Ns will merge in less than 10 yr, thus appearing as a LIGO signal (not shown here). However, still, a large part of the parameter space, with merger timescales between 10 and 10⁴ yr, has S/N above 1000. A strong signal is also depicted in the NS–BH examples (middle row), with S/N above ∼500. The majority of high S/N is expected to have a lifetime larger than 10 yr. However, with weaker signals, the NS binaries (bottom left panel) still have long-lived systems (merging in 10–10⁴ yr) with relatively large S/N (∼100). The S/Ns of WD binaries (bottom right panel) are the lowest and thus span the smallest regions of the parameter space. However, they are still detectable (with S/N of ∼50).

We note that the signal of a system that merges within T_obs < 4 yr cannot simply be described as the sum of the harmonics (e.g., Barack & Cutler 2004), as presented in Figure 1. In addition, in these parameter spaces, some large EKL-induced eccentricity oscillations can possibly take place. This situation is explored by Hoang et al. (2019), which showed that EKL eccentricity oscillations can be observed in LISA and may even infer the existence of an SMBH. In Figure 1, most of the systems with large S/N lay above the t_GW = 10 yr, which renders our use of Equation (8) valid (while the observation time is 4 yr, we use the 10 yr line to guide the eye).

Here we outline the procedure of population synthesis of binaries at the GC. We begin with adopting binaries' distribution properties from Stephan et al. (2016, 2019) of post-EKL systems at the onset of Roche limit crossing. We then generate 10⁶ systems at this stage and use COSMIC (Breivik et al. 2020a) to evolve them with time. Finally, we calculate the GW signal of the resulting compact object binaries via the LISA detector and estimate their LIGO detection rate.

4.1. Post-EKL Population

As mentioned above, recently Stephan et al. (2016, 2019) investigated the stellar binary evolution in the vicinity of an SMBH, including the EKL mechanism, tides, GR, and stellar evolution. For stellar binaries overflowing their Roche limit, the gas flow captured by a stellar companion will result in mass transfer, which renders subsequent stellar evolution different than a single-star evolution (e.g., Hurley et al. 2002; Shappee & Thompson 2013; Toonen et al. 2016; Antonini et al. 2017). Stephan et al. (2019) followed this population using the COSMIC stellar evolution code (see Section 4.2) and predicted a high abundance of compact object binaries within the GC. This formation channel seems promising to generate GW sources within the GC (e.g., Antonini et al. 2010; Antonini & Perets 2012; Hoang et al. 2018; Stephan et al. 2019).

We employ the same numerical setup and initial distribution as Stephan et al. (2016, 2019) and generate 10⁶ systems at the onset of Roche limit crossing. Here, we define the Roche limit of a binary system with Equation (11) (taking m₂ < m₁ and η = 1.6; e.g., Naoz 2016),

$$\begin{eqnarray}&&{R}_{\mathrm{oche}}=\eta {R}_{2}{\left(\displaystyle \frac{{m}_{2}}{{m}_{1}+{m}_{2}}\right)}^{-1/3}.\end{eqnarray} \tag{ 11 }$$

It is worth describing the Stephan et al. (2016, 2019) initial conditions that led to the Roche limit crossing we analyze. In those studies, the primary stellar initial mass function (IMF) is taken from Salpeter (1955), with α = 2.35 and a mass limit between 1 and 150 M_⊙. The mass ratio to the secondary mass in the binary uses that of Duquennoy & Mayor (1991). In other words, the mass ratio m₁/m₂ was taken from a Gaussian distribution with a mean of 0.23 and a standard deviation of 0.42. The SMBH of Sagittarius A* was set to 4 × 10⁶ M_⊙ (e.g., Ghez et al. 2005; Gillessen et al. 2009). The inner binary semi-axis distribution uses the same as that of Duquennoy & Mayor (1991), while the outer orbital period was distributed uniformly in log with a maximum of 0.1 pc. The eccentricity of the inner binary is uniformly drawn from 0 to 1 (Raghavan et al. 2010), and the eccentricity of the outer binary uses a thermal distribution (Jeans 1919). Both the inner and outer arguments of periapsides were taken from a uniform distribution between 0 and 2π. Additionally, those binaries generated must satisfy orbital and analytical stability (see Naoz & Fabrycky 2014; Naoz 2016; Stephan et al. 2016, 2019, for a complete specification). Using these initial conditions, Stephan et al. (2016, 2019) evolved 1500 binaries via the EKL mechanism, tides, GR, and single-stellar evolution. When the binaries crossed their Roche limit, Stephan et al. (2019) used COSMIC to follow their binary stellar evolution.

We use the Stephan et al. (2016, 2019) distribution, as stated above, at the onset of Roche limit crossing, to generate a population of 10⁶ binaries. We show the distribution of the systems' eccentricity and separation in Figure 2. The two peak distributions of the semi-major axis and eccentricity are expected owing to the EKL and tidal evolution (e.g., Fabrycky & Tremaine 2007; Naoz & Fabrycky 2014; Rose et al. 2019).

Zoom In Zoom Out Reset image size

We note that the additional distribution of wide binaries that do not cross their Roche limit may become unbound owing to interaction with single stars in this dense environment (e.g., Stephan et al. 2016; Rose et al. 2020) and thus will not contribute to the compact object binary population. We also note that other formation channels such as BH–BH collisions (e.g., O'Leary et al. 2009) or NS–BH capture (e.g., Hoang et al. 2020) will have only a small contribution to the overall population, but they explore different parts of the parameter space (e.g., Section 2).

4.2. Binary Stellar Evolution and Compact Object Population

As mentioned previously, COSMIC is a stellar evolutionary code that follows not only individual stars beyond their main-sequence evolution after they exhaust their hydrogen fuel but also interactions with a binary companion via mass transfer, the common envelope. Breivik et al. (2020a) presented the COSMIC code and evolved isolated stellar binaries all the way to their compact object stage. COSMIC uses a modified version of BSE that includes line-driven winds dependent on metallicities (e.g., Vink et al. 2001; Meynet & Maeder 2005; Gräfener & Hamann 2008; Vink et al. 2011), mass-loss transfer via Roche limit overflow (e.g., Hurley et al. 2002; Belczynski et al. 2008; Toonen et al. 2014), SN and natal kicks (e.g., Fryer & Kalogera 2001; Kiel et al. 2008; Fryer et al. 2012), etc. The output of BSE provides users the types of binary, binary state, and separation during each evolutionary stage. Breivik et al. (2020a) found around 10⁸ isolated compact binaries within the Milky Way as potential sources of electromagnetic and GW sources. Out of those, 10⁴ systems may be resolvable via LISA (Breivik et al. 2020a).

Since the publishing of Stephan et al. (2019), COSMIC has been modified to include the effects of ultra-stripped SNe that experience reduced natal kick strengths owing to the lack of a hydrogen envelope (Tauris et al. 2015). In comparison to standard core-collapse SNe, which have natal kicks drawn from a Maxwellian distribution with σ = 265 km s⁻¹ (Hobbs et al. 2005), COSMIC assumes that natal kicks in ultra-stripped SNe are drawn from a Maxwellian distribution with σ = 20 km s⁻¹. This inclusion increases the local merger rate of BNSs such that it is consistent with the 90% credible bounds of the LIGO/Virgo empirical rate (Zevin et al. 2020).

As a first step, we investigate the effect our post-EKL initial conditions have on the final compact object binary population, that is, comparing the evolution of isolated field binaries and binaries within the GC. Thus, we generate and evolve 10⁶ isolated stellar binary systems following Breivik et al. (2020a). We adopted an IMF from Kroupa et al. (1993), ⁷ a thermal eccentricity model, and a log uniform orbital period. We assume that all binary parameters are initialized independently and assume a binary fraction of 50%. ⁸ All binaries are evolved for 13.7 Gyr to capture their behavior over the age of the universe. As a proof of concept we focus on WD binaries. The final distribution of those isolated WD binaries is shown in Figure 3 (blue).

Zoom In Zoom Out Reset image size

We also considered the binaries that are on the onset of Roche limit crossing after undergoing the EKL mechanism (e.g., Figure 2). These binaries are evolved following the default settings defined in V3.3 of COSMIC (Breivik et al. 2020b). Note that high eccentricity is excited owing to the EKL mechanism induced by the SMBH. Based on the expectation of the eccentricity distribution after EKL (in the absence of tides; e.g., Teyssandier et al. 2013; Li et al. 2014a; Naoz 2016; Rose et al. 2019), we adopt a uniform distribution for binaries close to the SMBH as a proof of concept here.

In Figure 3 we focus on WD binaries, which represent the largest population in the system, and compare their population near an SMBH (green) to that of isolated binaries (blue, which their eccentricity does not excite). We also explore two metallicity populations: one solar metallicity (Z = 0.02), and one that is 15% of solar (Z = 0.003). As depicted, regardless of which metallicity we adopt, the overall distribution of the semimajor axis (bottom row of Figure 3) and the specific angular momentum (a(1 − e²), top row) of isolated binaries are much wider than those of binaries within the GC. The suppression of the widely separated binary population at the GC is a consequence of frequent interaction with neighboring objects in this dense environment, which tends to unbind the binary (e.g., Stephan et al. 2016; Rose et al. 2020).

We show the semilatus rectum of populations of WD–WD, BH–BH, BH–NS, and NS–NS binaries for metallicities Z = 0.02 and Z = 0.003 in Figure 4. Filtering out all systems that have either unbound, merged, or filled their Roche limit, the percentages of each compact binary type evolved out of 10⁶ systems are described in Table 2. The discrepancies between our results and those of Stephan et al. (2019) are potentially due to different assumptions for wind and ultra-stripped SNe in COSMIC as specified previously. In particular, the inclusion of ultra-stripped SNe explains the high abundance of NS–NS binaries in our simulations compared to Stephan et al. (2019).

Zoom In Zoom Out Reset image size

Table 2. EKL Merger Date and LISA's Detectability

Z	Population	Fraction	EKL Merger	LISA Band
			Rate†	Visibility
		(%)	(Gpc−3 yr−1)	(#)
0.003	BH–BH	1.00	5	2–19
	NS–NS	0.03	0.15	0.2–2
	NS–BH	0.05	0.25	0.025–0.25
	WD–WD*	4	19	15–148
0.02	BH–BH	0.06	0.3	0.3–3
	NS–NS	0.06	0.3	0.4–4
	NS–BH	0.06	0.3	0.2–2
	WD–WD*	3	15	14–137

Note. The third column presents the fraction of compact objects from 10⁶ systems after evolving with the BSE code in COSMIC. The fourth column is the EKL-induced LIGO detection rate of those sources with f_SMBH = 1 and SFR = 10⁻³ yr⁻¹. The fifth column gives the estimated number of detectable events within the LISA band during its lifetime per metallicity. †We note that the actual merge rate may be large because hard binaries will merge simply as a result of GW emission, with no need for EKL assistance. The last column presents the number of detectable sources of compact binary per galaxy via LISA. *We note that WD–WD merges are not visible in LIGO/Virgo Detector.

Download table as:  ASCIITypeset image

As apparent from Figure 4, metallicity significantly affects the abundance of compact objects. Metallicity as a function of stellar wind strongly affects the mass loss and final evolutionary outcome of massive stars (e.g., Kudritzki & Puls 2000; Nugis & Lamers 2000; Heger et al. 2003). In particular, subsolar metallicity yields a higher population for BH and WD binaries but a lower population for NS binaries. For massive BH–BH binary progenitors, high metallicity enhances the solar wind, leading to smaller helium core, smaller hydrogen envelope, and larger mass loss (e.g., Langer 1989; Hamann et al. 1995). Hence, BH–BH binary progenitors with masses >40 M_⊙ and high metallicity might cause them to fallback until only NSs are made (Heger et al. 2003). For less massive stars (≲10 M_⊙) evolving into WDs, lower metallicity allows AGB stars to develop into more massive cores with higher luminosity, which increases the survival chances of a WD binary. In comparison, higher metallicity reduces the mass of the helium core and increases the timescale for main-sequence stars to form WDs (e.g., Heger et al. 2003; Yu & Jeffery 2010; Romero et al. 2015). For compact binaries involving NSs (bottom row of Figure 4), their population fluctuation is less sensitive to metallicity. This is due to the lower mass-loss rate of their lower-mass NS progenitors relative to BH progenitors. In addition, their natal kicks, which might unbind the binary, correlate little with their metallicity (e.g., Giacobbo & Mapelli 2019; Neijssel et al. 2019).

4.3. LISA Detections

The post-evolutionary results from COSMIC presented here can be used in understanding the census of compact object binaries in the GC. Particularly, those compact object binaries have a chance to produce observable GW signals (see Section 2). Those eccentric binaries emit a GW signal that peaks approximately at

$$\begin{eqnarray}&&{f}_{p}(a,e)=\displaystyle \frac{{\left(1+e\right)}^{1/2}}{{\left(1-e\right)}^{3/2}}\times {f}_{\mathrm{orb}}(a),\end{eqnarray} \tag{ 12 }$$

where f_orb is defined in Equation (7) (e.g., Hoang et al. 2019). We assume that the outer semimajor axis of the triple hierarchy is much wider than the inner orbit, so that any GW back-reactions on the outer orbit are disregarded. In addition, we neglect oscillations in other orbital parameters due to EKL and GR since we expect them to only have peripheral effects on the f_p and the modulation of the inclination i and precession of pericenter are small near the high-eccentricity spike (e.g., Naoz et al. 2013). f_p , combined with the chirp mass of the binaries, are used to explain the distribution of strain curvers within the LISA parameter space.

As mentioned before, since the inner orbit of the triple hierarchy system is expected to undergo the EKL mechanism, we reassign a uniform distribution of eccentricity between 0 and 1 to the output binaries of COSMIC. The luminosity distance is set to be D_l = 8kpc = 1.65 × 10⁹ au. As previously discussed, we adopt an observation time that is the minimum of LISA's 4 yr mission duration or the binary's merger timescale. Here we explore the two evolved systems of binaries via COSMIC in Section 4.2 with Z = 0.02 and Z = 0.003. The detectability of each system is visualized by overplotting each system's GW strain curve, square root of Equation (9), onto the LISA sensitivity curve (red, Figures 5 and 7). ⁹ We note that out of 10⁶ binaries, no binary with S/N > 5, of any types in either metallicity, merges within 10 yr, which renders the S/N calculation from Equation (10) an appropriate approximation.

Zoom In Zoom Out Reset image size

We note that previous studies on compact object binary populations (e.g., Belczynski et al. 2010; Yu & Jeffery 2010; Liu et al. 2014; Lamberts et al. 2018; Lau et al. 2020; Sesana et al. 2020; Breivik et al. 2020a) have focused on simulating isolated binary sources of LISA in the Milky Way, including disk, bulge, and halo (e.g., Nelemans et al. 2001; Nissanke et al. 2012; Korol et al. 2017), while Kremer et al. (2018) study the population from Milky Way globular clusters. Here our simulation of eccentric compact object binary sources within the GC provides similar yet distinctive results for the population potentially detectable by LISA. Since the eccentricities of each source change under the EKL mechanism, we expect the actual number of observable sources to fluctuate up and down slightly from our current result.

In Figure 5, we consider BH–BH binaries with the orbital parameters shown in Figure 4 and color-code the S/N, with gray being below 5. The main difference between the solar and subsolar metallicity is the number of systems that are potentially detectable, where the subsolar calculation yields larger abundance of detectable systems. This result is consistent with isolated binaries in the field (e.g., Breivik et al. 2020a) and can be further understood with the aid of Figure 6. Figure 6 depicts the chirp mass, defined as

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{{\left({m}_{1}\times {m}_{2}\right)}^{3/5}}{{\left({m}_{1}+{m}_{2}\right)}^{1/5}},\end{eqnarray} \tag{ 13 }$$

plotted against the specific angular momentum of the systems. The mass of a system, which determines the overall amplitude of its GW signal (h₀, Equation (5)), is visualized via the chirp mass. As shown in Figure 6 and Table 2, not only is the overall population of the solar-metallicity systems 5% of that of the subsolar ones, but they also have smaller chirp masses. In other words, the subsolar-metallicity systems have higher BH masses that result in higher peak frequency, thus rendering them in the LISA detectable region.

Zoom In Zoom Out Reset image size

For NS–BH systems (top row in Figure 7), we found that the total number of systems did not change much between the two metallicity values we adopted (e.g., Figure 4). But due to the larger chirp mass and lower orbital frequency, their overall GW amplitude is smaller. Note that if NS–BH is formed via GW capture (e.g., Hoang et al. 2020), their separation will be much smaller than the post-EKL and binary stellar evolution prediction that is presented here. Therefore, these systems may still be detectable (see Figure 1). Additionally, stronger natal kicks in higher metallicity might lead to lower separation (e.g., Lu & Naoz 2019), which results in more visible systems than presented in Figure 7(b).

Zoom In Zoom Out Reset image size

For NS–NS systems, the simulated observable systems have S/N closer to 10 in either metallicity (Figure 7, middle row). While more systems are observable with higher S/N compared to NS–BH systems, they are less abundant than the BH–BH systems. This behavior is due to the concentration of systems within the low-frequency region $\mathrm{log}{f}_{p}\lt -4$, of which only higher-eccentricity pumps due to EKL might render some systems observable. Furthermore, note that as can be seen in Figure 4, the solar metallicity and subsolar metallicity differ only by about a factor of 2, resulting in a similar signal.

Maybe the most interesting feature we find is for WD–WD systems, depicted in the bottom row of Figure 7. First, it is worth noting that the WD binaries are predicted to be at the highest abundance. That is not surprising owing to the Salpeter IMF adopted by Stephan et al. (2016, 2019). In this IMF, the low-mass stars are at high abundance, which is also consistent with the Stephan et al. (2019) results. Additionally, the semimajor axes of the WD binaries are relatively short (as can be seen in Figure 2), which raises their orbital frequency to that of LISA (log f_p > −3). Lastly, as mentioned above, we adopt a uniform eccentricity distribution that yields a smaller specific angular momentum. Therefore, the combination of all of these factors results in a large abundance of single detections of WD binaries at the GC.

We also offer an approximation for the number of visible sources per galaxy as

$$\begin{eqnarray}&&\begin{array}{l}{\rm{N}}={N}_{\mathrm{steady}\_\mathrm{state}}\times {f}_{\mathrm{Roche}}\\ \qquad \times \,{f}_{{\rm{S}}/{\rm{N}}\ \gt \ 5}\times {f}_{\mathrm{binary}\_\mathrm{fraction}},\end{array}\end{eqnarray} \tag{ 14 }$$

in which

$$\begin{eqnarray}&&{N}_{\mathrm{steady}\_\mathrm{state}}=\mathrm{SFR}\times {t}_{\mathrm{steady}\_\mathrm{state}},\end{eqnarray} \tag{ 15 }$$

where we assume a star formation rate ¹⁰ of SFR = 10⁻³ yr⁻¹ and take f_Roche ∼ 0.25 as the fraction of binaries that merge via Roche limit crossing (Stephan et al. 2019). Additionally, N_{steady_state} is the number of stars in steady state at the galactic nuclei. Motivated by our own Galaxy, we adopt N_{steady_state} ∼ 10⁵–10⁶. Lastly, f_{binary_fraction} is the fraction of each type of binary specified in Column (3) of Table 2, and f_S/N>5 is the fraction of visible binaries in LISA out of each type of compact object binary. Using this estimation, we find that, per galaxy, the EKL mechanism may yield about 14–150, 0.02–2, 0.2–4, and 0.25–20 of WD–WD, NS–BH, NS–NS, and BH–BH binaries, respectively, that could be visible within the LISA band (the full breakdown can be found in Table 2).

4.4. LIGO Detections

Additionally, massive mergers such as BH–BH and NS–BH are also potential signal sources of LIGO (e.g., Fragione & Kocsis 2019; Fragione et al. 2019a, 2019b). To allow for quantitative comparison with the results of Stephan et al. (2019), we adopted their parameters. In what follows, we assume a similar galactic condition to that of the Milky Way to estimate potential signals from compact object binaries of other galaxies via LIGO up to ∼1.5 Gpc (Abbott et al. 2018). We assume a galaxy density of 0.02 Mpc⁻³ (Conselice et al. 2005) and a star formation rate of SFR = 10⁻³ yr⁻¹. The latter is estimated based on the Milky Way properties (e.g., Genzel et al. 2003; Schödel et al. 2003). The fraction of galaxies with an SMBH in the center (f_SMBH) is approximated between 1 and 0.5 (e.g., Ferrarese & Ford 2005; Kormendy & Ho 2013), where we adopt a unity value. Additionally, Stephan et al. (2019) predicted that the fraction of binaries that will merge via crossing the Roche limit is f_Roche ∼ 0.25. Combining, the LIGO detection rate, Γ, of each type of binary within 1 Gpc³ via EKL merger channel is

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Gamma }}=\left(\displaystyle \frac{0.02\times \mathrm{Galaxies}}{{\mathrm{Mpc}}^{3}}\right)\left(\displaystyle \frac{{10}^{9}\times {\mathrm{Mpc}}^{3}}{{\mathrm{Gpc}}^{3}}\right)\times {f}_{\mathrm{SMBH}}\\ \qquad \times \,\mathrm{SFR}\times {f}_{\mathrm{Roche}/\mathrm{merger}}\times {f}_{\mathrm{EKL}}\times {f}_{\mathrm{binary}\_\mathrm{type}}\,,\end{array}\end{eqnarray} \tag{ 16 }$$

where, f_EKL is the fraction of the compact object binaries that will merge owing to EKL, conservatively estimated as f_EKL ∼ 0.1 (e.g., Hoang et al. 2018). The binary type, f_{binary_type}, is specified in the third column of Table 2. We provide the summary of the possible detectable events in the fourth column of Table 2. For the BH binaries, those results are within a factor of ∼3–4 from the Stephan et al. (2019) result using Z = 0.02 but fairly consistent with their estimates for NS–BH. As mentioned above, this is a result of the updated version of COSMIC and the assumptions of the natal kicks in ultra-stripped SNe. We note that in Stephan et al. (2019) there were no NS–NS present, but with this updated COSMIC package (see above), more NS–NS are present. In addition, the binary stellar evolution prescription here creates a large uncertainty when estimating the possible rates.

In Figure 8 we depict the distribution of BH binaries mass ratio after undergoing the COSMIC evolution. As shown in this figure, both distributions have wings of mass ratio deviating from 1:1, with the subsolar population possessing a larger variance. Interestingly, previous LIGO/Virgo observations have found a mass distribution that indicated similar mass components (e.g., Abbott et al. 2016a, 2016b, 2016c, 2016d, 2017a, 2017b, 2017c, 2017d, 2019b, 2020a, 2020b). However, two events GW190412 and GW190814 (e.g., LIGO Scientific Collaboration et al. 2020; Abbott et al. 2020c) possess mass ratio that differs significantly from unity. As mentioned above, the mass ratio distribution in our model, following the Stephan et al. (2019) setting, was taken from a Gaussian distribution with a mean of 0.23 and a standard deviation of 0.42. Thus, while the deviation from the mass ratio of unity is an imprint of the initial conditions, it is interesting to note that the subsolar-metallicity population has a larger range of mass ratios.

Zoom In Zoom Out Reset image size

Almost every galaxy, including the Milky Way, hosts an SMBH in its center, surrounded by crowded nuclear stellar clusters. The proximity of the Milky Way's SMBH provides a unique opportunity to explore dynamics and phenomena that ought to exist in other galaxies. It has been suggested that stellar binaries are of high abundance around the nuclear cluster of the GC (e.g., Ott et al. 1999; Martins et al. 2006; Pfuhl et al. 2014; Stephan et al. 2016, 2019). Within the vicinity of an SMBH a binary has to be on a tighter configuration compared to its orbit around the SMBH, and hence it undergoes the EKL mechanism (e.g., Antonini et al. 2010; Antonini & Perets 2012; Prodan et al. 2015; Stephan et al. 2016, 2019; Hoang et al. 2018). Recently, Stephan et al. (2019) showed that combining EKL evolution with single and binary stellar evolution yields a high abundance of compact object binaries and has the potential to become GW sources for the LIGO/Virgo and LISA. Here, we simulate the potential observable compact object binaries in the vicinity of an SMBH under the effect of EKL via the GW interferometers.

Various merger channels within the GC, besides the EKL-assisted merger channel, also contribute to the formation of compact object binaries. Those include binary–single and binary–binary mediated interactions (e.g., Rodriguez et al. 2018; Arca Sedda 2020) and single–single GW captures (e.g., O'Leary et al. 2009; Tsang 2013; Hoang et al. 2020). Figure 1 is agnostic to the merger channel of those compact objects and highlights that a vast part of the parameter space can be detectable via LISA.

However, we note that for NS–BH and BH–BH binaries the EKL mechanism may be a significant contributor (e.g., Hoang et al. 2018, 2020). Focusing on this merger process, we adopt the initial distribution from Stephan et al. (2019) after EKL distribution and generate a large population at the onset of Roche limit crossing. Using COSMIC, we then evolve these binaries to their compact object stage. As expected, the compact object binary population near an SMBH has different properties, as highlighted in Figure 2, than field binaries. In particular, these binaries possess shorter semimajor axis compared to isolated field binaries, due to interaction with passing objects at this dense environment (Rose et al. 2020). Additionally, we expect the binaries' eccentricity distribution to be excited owing to the EKL mechanism (e.g., Naoz & Fabrycky 2014).

We have tested both solar and subsolar metallicities. As expected, the subsolar metallicity produces more BH–BH binaries, by a factor of ∼18. The population of NS–NS binaries for solar metallicity is about double that of subsolar metallicity. We note that although recent studies reported that the GC may have a supersolar population (e.g., Feldmeier-Krause et al. 2017; Do et al. 2018), the hierarchical formation of galaxies (e.g., White & Rees 1978) suggested that subsolar metallicities should exist at high abundance. Adopting subsolar and solar metallicity allows us to also extrapolate the conditions of these sources to other galactic nuclei. The exploration of supersolar metallicity is beyond the scope of this paper.

From Figures 5 and 7, we see that over the observation interval of 4 yr, BH–BH binaries are the dominant observable sources with the highest S/N of GW signals (S/N > 100, Figure 5), while WD–WD binaries produce the most abundant number of systems mostly clustered around S/N ≲ 10 (Bottom row of Figure 7). Those two types of compact object provide the most promising GW sources of LISA via the EKL merging channel. Particularly, the WD–WD binaries, with their low separation and specific angular momentum (Figure 3), will lead to a strong overlap of orbital frequency and the LISA sensitivity curve. Not only will they contribute to noise confusion, but they might also be detected individually. In comparison, NS–NS and NS–BH binaries are less likely to be observed within this formation channel but possibly abundant through other means, which we leave to explore in future studies.

We also estimate the number of binaries, per galaxy, that are expected to be visible within the LISA band. While the details are highly uncertain, we are motivated by our own Galaxy, for example, which has about a million stars within the inner parsec. We note that this is a conservative estimation because further from the SMBH unbending of the binary due to interactions with neighboring stars is less efficient (Rose et al. 2020); thus, we expect an even larger population than these conservative numbers. Overall we find, depending on the metallicity, about 140–150 WD–WD, 0.2–2 NS–BH, 2–4 NS–NS, and 2–20 BH–BH to be visible within the LISA band.

While we adopted a D_l ∼ 8 kpc for our GC, since the amplitude of the GW signal is linear with the luminosity distance (Equation (8)), a detection of compact binaries in Milky Way–like galaxies with distance >1 Mpc is possible with a longer observation timescale (T_obs ≥ 4 yr). Hence, our results can also be extended to binary sources marginally observable via LISA at this luminosity distance. Particularly, Hoang et al. (2019) have shown that those systems have a chance to be detected via the eccentricity oscillations through the gravitational perturbation from the SMBH in the galactic nuclei, i.e., EKL. For a galaxy that contains a more massive SMBH, some systems might be visible via the eccentricity oscillation with only T_obs ∼ 1 yr (e.g., Emami & Loeb 2020).

The compact object binary population we found may also merge via either GW emission or EKL as time goes by. We roughly estimate the rate that the LIGO/Virgo can detect (see Table 2). These rates are sensitive to the SFR, where we assume a very conservative rate. However, E + A galaxies, or starburst galaxies, may undergo star formation episodes that could possibly increase the stellar population by ∼10%, as well as increase the star tidal disruption events (e.g., Dressler & Gunn 1983; Swinbank et al. 2012; Arcavi et al. 2014; Stone & van Velzen 2016).

We have shown that the GC and other galactic nuclei are potentially significant sources for LISA and LIGO/Virgo. Most importantly, compact object binaries at the GC can have an extremely large S/N in the LISA band. Furthermore, the GW signal from the GC predicted here will have a preferential direction, compared to other detections in LISA, which may assist in disentangling the various signals.

We thank Salvo Vitale, Mark Morris, and Tuan Do for useful discussion. H.W. thanks the UCLA-summer REU program. A.P.S. acknowledges partial support by the Thomas Jefferson Chair Endowment for Discovery and Space Exploration and partial support through the Ohio Eminent Scholar Endowment. H.W., A.P.S., S.N., and B.-M.H. acknowledge the partial support of NASA grant Nos. 80NSSC20K0505 and 80NSSC19K0321 and NSF through grant No. AST- 1739160. S.N. thanks Howard and Astrid Preston for their generous support.