Probing Intermediate Mass Black Holes in M87 through Multi-Wavelength Gravitational Wave Observations

We analyse triple systems composed of the super massive black hole (SMBH) near the center of M87 and a pair of black holes (BHs) with masses in the range $10-10^3$ $M_{\odot}$. We consider the post Newtonian precession as well as the Kozai-Lidov interactions at the quadruple and octupole levels in modeling the evolution of binary black hole (BBH) under the influence of the SMBH. Kozai-Lidov oscillations enhance the gravitational wave (GW) signal in some portions of the parameter space. We identify frequency peaks and examine the detectability of GWs with LISA as well as future observatories such as $\mu$Ares and DECIGO. We show examples in which GW signal can be observed with a few or all of these detectors. Multi-Wavelength GW spectroscopy holds the potential to discover stellar to intermediate mass BHs near the center of M87.


INTRODUCTION
Black holes (BHs) are ubiquitous in galaxies and form across a wide range of masses and environments. Stellar-mass BHs are most abundant, with observational evidence first detected in X-rays (Webster & Murdin 1972;Remillard & Mc-Clintock 2006), and more recently by gravitational waves (GWs) with LIGO/Virgo (Abbott et al. 2016a,b). Intermediate mass black holes (IMBHs) are expected to form from the accretion of gas in dwarf galaxies (see (Reines et al. 2019) and references therein), mergers of stars in dense stellar clusters (Portegies Zwart et al. 1999;Devecchi & Volonteri 2009;Pan & Loeb 2012;Mapelli 2016), from direct collapse of inflowing dense gas in protogalaxies (Loeb & Rasio 1994), collapse of Pop III stars from early universe (Madau & Rees 2001;Bromm & Larson 2004) or from supermassive stars in AGN accretion disk instabilities (McKernan et al. 2012(McKernan et al. , 2014. Finding unambiguous observational evidence for IMBHs in galactic nuclei is challenging due to the short lifetime associated with mergers and accretion by the supermassive black holes (SMBHs) there (Natarajan 2014;Johnson & Haardt 2016). There is however some evidence for their existence in the local Universe (Mezcua 2017). In particular, Ultraluminous X-ray sources (ULXs) imply BHs with masses above 20 M in some cases (Miller & Colbert 2004). A recent discovery of IMBHs in dwarf galaxies at (z 2.4) was reported in the Chandra COSMOS Legacy Survey (Mezcua et al. 2018), using a sample of 40 active galactic nuclei (AGN) in dwarf galaxies with redshift below 2.4. More recently (Chilingar-ian et al. 2018) used data mining in wide-field sky surveys and identified a sample of 305 IMBH candidates in galaxy centers, accreting gas which creates characteristic signatures of a type I AGN. Most recently, a new sample of Wandering IMBH was discovered using the large-scale and high resolution radio telescopes (Reines et al. 2019) in nearby dwarf galaxies.
New probe of IMBHs are awaited using 30-m class telescopes (Greene et al. 2019). IMBHs may also be able to detect using the GW signals in different frequencies (Kocsis & Levin 2012).
SMBHs, with masses in the range 10 6 − 10 10 M are believed to exist in the core of almost all of massive galaxies (Wang et al. 2015;Broderick et al. 2015). This includes SgrA* at the Galactic center (Ghez et al. 1998;Gravity Collaboration et al. 2018), as well as the nearby elliptical galaxy M87 (Gebhardt & Thomas 2009;Gebhardt et al. 2011;Walsh et al. 2013). Most recently the mass of the M87 SMBH was precisely measured by Event Horizon Telescope (EHT) collaboration (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f) to be 6.5 × 10 9 M . Being at the center of the virgo cluster, M87 is the nearest giant elliptical galaxy which is the product of many galaxy mergers. As a result, M87 should contain a large population of stellar mass BHs and possibly a handful of IMBHs which used to be at the center of dwarf galaxies that merged in. Here we propose to search for this IMBH populations in M87 through their GW signals.
A large fraction of the IMBHs in M87 might result in binaries in eccentric orbits around the SMBH. This system is well described by the hierarchical triple approach (Naoz et al. 2013;Randall & Xianyu 2019;Hoang et al. 2019;Antonini et al. 2019) in which the inner binary contains of IMBHs while the outer binary includes the SMBH. The IMBHs in the inner binary are expected to emit bursts of GWs at any pericenter passage. The frequency of these GW bursts depend on the orbital parameters.
In this paper, we simulate triple systems composed of M87 and a pair of the BHs and consider the dynamical evolution of the binary BHs. The induced Kozai-Lidov oscillations (Naoz et al. 2013 The paper is organized as follows. In Sec. 2 we briefly review the details of the related hierarchical triple system. In Sec. 3 we introduce the stability conditions that must be taken into account. In 4.1 we review the finite time Fourier transformation for computing the GW signal. In Sec. 4.2 we present a variety of different examples with a potentially detectable GW signal as demonstrated in Sec. 5. In Sec. 6 we briefly consider binaries with non-equal BH masses. In Sec. 7 we compute the lifetime of the related systems. we discuss about the detectability of the GW signal. Finally, we summarise our conclusions in Sec. 8.

IMPACT OF SMBH ON THE EVOLUTION OF BBHS
We model the dynamical influence of a SMBH on the evolution of BBHs, accounting for inner pericenter precession at first order, quadruple and octupole terms in the Lidov-Kozai interaction.
The full Hamiltonian of the triple system is given by, where H LK and H 1P N present the Lidov-Kozai Hamiltonian and first order post Newtonian precession, respectively. H LK = H 1 + H 2 + H 12 with H i , i = 1, 2 referring to the Keplerian Hamiltonian for the inner and outer binaries in the system and H 1,2 denoting the interaction term between them. The interaction term is expressed as a series expansion in the separation of two binaries. We make use of Eqs. (5-8) of Rodriguez & Antonini (2018) for modeling the above components of full Hamiltonian. Using Eq.
(1) we derive the equations of motion for the inner and outer semi-major axes. We also take into account GW emission in the inner orbit as an extra term that shrinks the orbit of the BHs.

STABILITY CONDITIONS
Here we present some stability conditions that must be taken into account in our analysis: • The orbital parameters must be selected such that prevent the inner binary from reaching the Roche limit of the outer-binary (Hoang et al. 2018;Fragione et al. 2019). This implies, (2) • The system must possess dynamical stability, implying that the hierarchical secular treatment is satisfied (Naoz & Silk 2014;Hoang et al. 2018), (3) • For eccentric outer orbits, the pericenter distance must be much larger than the event horizon of SMBH,

GRAVITATIONAL WAVES ESTIMATION
Given the dynamical evolution of the orbital parameters, we calculate the GW amplitude from eccentric binary black holes (hereafter EBBHs) and study the detectability of their GW signal. Using this formalism, we present some examples with potentially detectable GW signals in M87.

GW-Computation
Unlike binaries on circular orbits, EBBHs emit GW in a discrete spectrum (Peters & Mathews 1963) with characteristic frequencies f n = n f orb , where n refers to the harmonic index while f orb = (2π) −1 G(m 1 + m 2 )a −3/2 is the orbital frequency in a circular orbit with semi-major a and m 1 and m 2 refer to the mass of BHs in the inner binary. The amplitude of the GW signal is given by, where h n (a, e, f n ) is defined as, with h 0 (a) being the dimensionless strain for a circular orbit, given by (Peters & Mathews 1963;Randall & Xianyu 2019;Hoang et al. 2019), and where D denotes the angular diameter distance to the source.
To gauge the detectability of the GW signal, we divide the stream of GW into time intervals with a duration ∆T and perform a Finite Fourier Transformation (FFT) of the GW strength for the duration ∆T. The Fourier component  of the GW signal can be computed by taking a time integral of Eq. (5) from −(∆T/2) to (∆T/2) yielding, where, The observational interval ∆T plays an important role. For close-in sources, the strength of GW is large enough to allow smaller values of ∆T. This leads to wider frequency bins, as f − f n 1/∆T. The situation is different for wide-separation binary systems, where the observational time ∆T must be larger. Therefore the GW signal is localized in the frequency domain around specific harmonics. We define characteristic detectable frequencies which are aimed to be within the frequency range of LISA or ground based GW detectors, where α 1 and with f n = n f orb .
Next, we define the characteristic strain of the GW signal as, To check the detectability of the GW signal, Eq. (11) must be compared with the characteristic noise of GW observatories. As discussed in the introduction, we will consider several different GW detectors, including LISA, µ-Ares and DECIGO.
Defining the noise generically as S c ( f ), the signal to noise ratio (S/N) is given by,

GW signal in M87
Having presented a formalism to compute the GW signal for a generic system, we now apply it to the case of M87.
Since D = 16Mpc, we choose a relatively large value for the observational time, ∆T. From Eq. (8), a large ∆T has two different effects, though. On the one hand, it enhances the GW amplitude. On the other hand, it diminishes the frequency width of the GW signal as ( f − f n ) 1/∆T. This leads to a discrete spectrum of observable frequencies. The combination of these two effects lead to a range of ∆T that can lead to an observable GW signal at a range of discrete frequencies.
To clarify the above points, we present some examples with potentially detectable GW signals at different frequencies. We choose ∆T differently in each of these examples to both help pushing the GW signals above the noise as well as increasing the amount of observable frequencies.
Throughout our analysis, we neglect the GW bursts from the inspiral phase of IMBHs around M87. In another paper (Emami & Loeb 2019), we estimated the timescale associated with these bursts to be, It is easy to see that τ GW O(Gyr) is much longer than both of the observational time and the lifetime of the orbit of IMBHs. Therefore, we can safely ignore the GW decay time in our analysis.

Probing GW in LISA
First, we consider GWs in the LISA band. We adopt the LISA noise curve (Robson et al. 2019;Emami & Loeb 2019) and present the characteristic GW signal on the top of that. Figure 1 shows the GW signal in the LISA band for few different equal mass binaries with each component having, m(≡ m/M ) = (10, 50, 100, 300, 500, 700, 1000). The rest of the parameters are chosen to be the same:ã in = 0.05,ã out = 4000, e in = 0.7, e out = 0.7, where hereafterã ≡ (a/AU). On the left panel, we use an interpolation for the points with π( f − f n )∆T < 1, whereas the right panel, presents the points satisfying π( f − f n )∆T < 1. We adopt ∆T = (8, 7, 5, 2, 1, 0.5, 0.5) yrs.

Probing GW with the µ-Ares detector
Since the LISA noise increases at frequencies below milli-Hz, it is advantageous to use the newly proposed µ-Ares detectors for detecting GW signals in the frequency range from µHz to milli-Hz. Figure 2 presents the characteristic strength of GW at those frequencies. For BH masses with m = (10, 50, 100, 300). The plot is shown forã in = 0.1,ã out = 4000, e in = 0.8, e out = 0.7. Comparing Figure 2 with Figure  1, reveals that the strength of GW form = 10 is enhanced for larger value of inner semi-major axes. This is due to the Kozai-Lidov oscillations which boost the GW signal above the noise toward larger frequencies. Here we consider ∆T = (10, 8, 5, 2) yrs.
There are clearly overlapping regions in frequency for these three observatories. This is particularly helpful in removing degeneracies between various parameters in the system. Multi-wavelength spectroscopy of GW can provide novel information about the parameters of the BBHs that are otherwise degenerate, and could be potentially used as a way to discover IMBHs in M87.

Time evolution of GW amplitude
Having presented the frequency evolution of the GW signal for a fixed time, we next study how the signal changes with time. In Figure 4, we draw the evolved GW in a wide range of frequencies. Here we adoptm = 300, e in = 0.6, e out = 0.6,ã in = 0.05,ã out = 6000. We considert = (1040,2040,3040,5040,6040,7632) for ∆T = (9, 8, 6, 5, 4, 0.1) month. In this example, increasing ∆T mostly affects the detectability of the GW signal for the DECIGO observatory.

Impact of orbital parameters in GW amplitude
In our simulations, we fix some of the orbital parameters such as the arguments of pericenter for the inner and outer binaries (taken to be 0 • ), longitude of ascending node for inner and outer binaries (chosen to be 0 • ) and mutual inclination (taken to be 90 • ). We have also taken the BHs to have zero spins, but allowed rest of the parameters to vary. This includes the inner and outer semi-major axes and eccentricities. We noticed that changing the outer semi-major axis has very minor impact on the results as long as the BBH is far from the tidal disruption distance, a out a in (M S M BH /M BH ) 1/3 . Similarly, changing e out does not affect the signal significantly. On the other hand, changing a in and especially e in affect the strength of the signal dramatically. Owing to the importance of these parameters, we consider their effect on the GW signal for multiple detectors. Figure 5 presents the influence of e in on the GW signal, assumingm = 300, e out = 0.6,ã in = 0.1 andã out = 4000. The resulting GW signal could be observed by µAres, LISA or DECIGO detectors.
Finally, in Figure 6 we examine the impact of changing a in = (0.03, 0.05, 0.07, 0.1, 0.15, 0.2) on the detectability of GW signal. Here we adoptm = 300, e in = 0.7, e out = 0.6,ã out = 6000. The plot shows that changing a in only affects slightly the GW signal.

DETECTABILITY OF GW SIGNAL
Next, we consider the detectability of GW signals for some of the examples above. We compute the signal to noise ratio (S/N) using Eq. (12), and label a GW signal as detectable    if S/N ≥ 10. Since different GW detectors are focused on different frequency bands, we need to define the frequency bands for different detectors using Eq. (12). Table 1 presents the frequency ranges for different GW detectors. Although for most cases we take the universal upper and lower limits in the integrals, the lower limit for µAres depends on the minimum value of the frequency, which could be slightly larger than 10 −6 Hz. Owing to this, we take the maximum value between 10 −6 Hz and f m , defined to be the minimum value of the frequency of GW signal. Likewise, for the LISA experiment, we take the minimum value of the frequency of signal as the lower limit of the integral. In addition, since the GW frequency is chirped towards the merger, the computation of the SN also depends on the approximate point in the evolution of the system. In other words, the time evolution of the system affects the GW signal and so the S/N. Therefore, lower signal to noise may evolve with time and get enhanced. With this in our mind, in the following we present S/N ratio for one set of the examples. Table 2 presents the signal to noise ratio for the case withm ≡ m/M = (10, 50, 100, 300, 500, 700, 1000) and withã in = 0.1,ã out = 4000, e in = 0.8, e out = 0.7, and ∆T = (9, 8, 5, 3, 2, 1, 1, 1)yr. The GW signal is evaluated at the same time for all cases. Since the dynamical evolution of different BHs differ depending on their masses, the orbital and  peak freak frequency, is different in these examples. This changes the location of the signal and so depending on the details of the evolution S/N ratio changes by changing the evaluation time at different BH masses, it is possible to make the peak frequency and GW signal similar.

NON-EQUAL MASS BHS
So far, we only focused on the BBHs with equal masses.
Here we briefly consider the case with non-equal BH masses. As a test example, we fix the mass of one of BHs to be m 1 = 100M and change the mass of its companion in the range m 2 = (30, 50, 100, 300, 500, 700)M . The rest of the orbital parameters are taken to be a in = 0.1AU, e in = 0.4, e out = 0.6, a out = 6000AU. Figure 7 presents the GW signal for this regime.
Unlike the examples in Sec. 5, we evaluate the GW signal at the time for which the peak frequency (defined in Eq.  14) to be around f p = 3 × 10 −4 Hz. This makes the GW signal behave very similarly. Therefore the peak frequency is a key parameter in our system and leads to an almost universal behavior at different BH masses. In closing, we note that each IMBH may carry a cluster of stellar-mass BHs around it, enhancing the rate of detectable GW signals from its vicinity.

LIFETIME OF ORBITS
Finally, we consider the lifetime of BBHs orbits under consideration. Figure 8 presents the lifetime of BBHs for different values of their masses and as a function of e in and a in . For simplicity, we only focus on equal mass BHs. From the plot it is clear that the lifetime of the BBHs is a strong function of the orbital parameters as well as the BH masses. As expected, the lifetime increases by decreasing the values of e in and increasing the value of a in .

CONCLUSIONS
Multi-Wavelength GW detectors can monitor the spectrum of GW signals from the center of M87 over a wide range of frequencies and orbital parameters. GW spectroscopy enables one to reproduce the shape of the GW signal and get novel information about the physical process behind such signals. Focusing on triple systems in M87 made of an inner binary BHs with different masses, from stellar to IMBHs, we presented a consistent method for detecting GW signals by integrating over the observation time within the lifetime of different GW detectors. We demonstrated that the frequency peaks from various GW sources can be used to entangle the signal from closer in sources with continuum frequency bands. The parameter space where different GW detectors may overlap in their frequency of GWs, could be used as a novel way to break the degeneracy between different orbital parameters. We have neglected the possible encounters between the BBHs and individual BHs. A more detail numerical simulation over a wider time range including all of the possible encounters between different BHs is for a future investigation.