Detecting gravitational-wave bursts from black hole binaries in the Galactic Center with LISA

The Laser Interferometer Space Antenna (LISA) is a planned space-borne gravitational-wave (GW) observatory (Amaro-Seoane et al., 2017; Colpi et al., 2024), currently set to launch in the late 2030s. With sensitivity in the $10^{-4}\,\,\mathrm{Hz}$ to $10^{-1}\,\,\mathrm{Hz}$ frequency range, LISA will unlock the source-rich mHz band of the GW spectrum, potentially detecting GWs from stellar-mass black hole binary (BHB) inspirals, massive BHB mergers, extreme mass-ratio inspirals, and millions of ultra-compact galactic binaries composed of white dwarfs and neutron stars (Amaro-Seoane et al., 2023). The ability of LISA to observe compact binaries long before they merge will present unique opportunities to study the dynamics of hierarchical multi-body systems (Amaro-Seoane et al., 2023), where the influence of the perturbing body can have a measurable effect on the gravitational waveform, thus revealing key insights into compact binary formation channels in dense stellar environments (Sigurdsson & Hernquist, 1993; Portegies Zwart & McMillan, 2000; McMillan & Portegies Zwart, 2003; Rodriguez et al., 2016; Belczynski & Banerjee, 2020; Gerosa & Fishbach, 2021; Bartos et al., 2017; Tagawa et al., 2020).

Galactic nuclei are expected to contain an abundance of stellar-mass black holes (BHs), many of which may exist in binaries. Two-body relaxation naturally causes heavier masses to migrate inwards, resulting in a dense compact object core where BHBs are readily assembled and hardened through repeated BH-BH encounters (Freitag et al., 2006; O’Leary et al., 2009, 2016; Antonini & Rasio, 2016; Rodriguez et al., 2018a, b; Samsing, 2018). The Milky Way is no exception, potentially hosting ${\sim}2\times 10^{4}$ BHs within the inner parsec of the Galactic Center (Morris, 1993; Miralda-Escudé & Gould, 2000; Freitag et al., 2006; O’Leary et al., 2009; Naoz et al., 2018; Rose et al., 2022). If the galaxy hosts a central supermassive black hole (SMBH), as is believed to be true for most large galaxies (Kormendy & Ho, 2013) including our own (Ghez et al., 2005; Gillessen et al., 2009), stellar-mass BHBs may become bound to the SMBH on a relatively wide outer orbit, forming a hierarchical triple (Antonini & Perets, 2012; Hoang et al., 2018; Grishin et al., 2018; Rose et al., 2020). We illustrate this special orbital configuration in Fig. 1. Gravitational torques exerted by the SMBH can subsequently induce secular eccentricity and inclination oscillations through the eccentric Kozai-Lidov (EKL) mechanism, driving highly inclined BHBs towards eccentricities of nearly unity (Kozai, 1962; Lidov, 1962; Naoz, 2016) and resulting in a GW source relevant for LISA. These novel sources possess a unique GW signature, as their extreme eccentricity concentrates the GWs into a train of bursts occurring once every periastron passage (O’Leary et al., 2009; Gould, 2011; Kocsis & Levin, 2012; Xuan et al., 2023b). Changes in the morphology and timing of the bursts reflect the secular evolution of the BHB due to the SMBH, allowing one to infer the orbital parameters of both the binary and its perturber (Romero-Shaw et al., 2023). Eccentricity notably enhances GW emission and accelerates the inspiral timescale (Peters & Mathews, 1963), possibly contributing to the compact binary merger rate observed by ground-based GW detectors (Hoang et al., 2018; Stephan et al., 2019; Lim & Rodriguez, 2020; Martinez et al., 2020; Abbott et al., 2023; Fragione et al., 2019).

Previous works have established that SMBH-induced BHB eccentricity oscillations are detectable by LISA out to a few Mpc (e.g., Randall & Xianyu, 2018; Hoang et al., 2019; Wang et al., 2021)¹¹1This effect can also take place in a wide range of masses and separations, see e.g., Randall & Xianyu (2019); Emami & Loeb (2020); Deme et al. (2020)., provided the oscillation timescale is comparable to the mission lifetime. The loudest BHBs are expected to be orbiting the $4\times 10^{6}\,\,M_{\odot}$ SMBH at the Galactic Center (Sgr A^∗), with an estimated few tens of detectable sources in the LISA band (Wang et al., 2021; Xuan et al., 2023b). Due to their close proximity and heavy masses, these binaries can accumulate signal-to-noise ratios (SNRs) between $10^{2}$–$10^{4}$ over the course of a $4\,\,\mathrm{yr}$ LISA mission (Hoang et al., 2019; Wang et al., 2021), necessitating efforts to model and subtract the bursts from LISA data (Littenberg et al., 2020; Littenberg & Cornish, 2023). In this Letter, we study the GW signature of BHBs orbiting Sgr A^∗, prioritizing the initial task of detection. We calculate the GW signal in the time domain, including the full instrument response of LISA, by incorporating simulated background noise to emulate a detection scenario as realistically as possible. As a proof-of-concept, we explore the use of a wavelet decomposition, based on the multi-resolution $Q$-transform (Chatterji et al., 2004; Chatterji, 2005), to resolve bursts from eccentric BHBs in LISA data and constrain their time-frequency properties. Additionally, we show that this time-frequency information provides a means to infer the orbital period and eccentricity of the BHB.

We begin by examining the GW signal originating from a stellar-mass BHB orbiting around a Sgr A^∗-like SMBH, as observed by LISA. In order to generate simulations of the expected signal, we calculate the time-domain GW polarizations $h_{+,\times}$ for the BHB using a leading-order adiabatic approximation described in Barack & Cutler (2004), which expresses the waveform as a sum over harmonics of the orbital frequency:

where $m_{1}$ and $m_{2}$ are the component masses, $D$ is the distance from the detector, and $\boldsymbol{\theta}$ denotes the standard (time-varying) Keplerian orbital elements. The mathematical definitions of the waveform harmonics $h_{+,\times}^{(n)}(t,\boldsymbol{\theta})$ are given in Barack & Cutler (2004). We evolve the system in time by numerically integrating the orbit-averaged equations of motion governing a hierarchical triple, including gravitational quadrupole and octupole-level contributions (Naoz et al., 2011, 2013), general relativistic precession (Naoz et al., 2013), and gravitational radiation (Peters, 1964). The explicit dynamical equations can be found in Naoz (2016). We then evaluate Eq. 1 by sampling the time-evolved orbital elements. We ensure stability by imposing that the BHB cannot exceed its Hill radius (Naoz & Silk, 2014; Hoang et al., 2018; Grishin et al., 2017; Tory et al., 2022), as well as to avoid the breakdown of the secular approximation (Luo et al., 2016; Grishin et al., 2018). Further, the binaries we consider have hardened enough such that we expect them to resist disruption by the local stellar environment (see e.g., Hoang et al. (2018); Rose et al. (2020) for discussion on the survivability of wide binaries in galactic nuclei).

In LISA, the output data streams are the time-delay interferometry (TDI) channels, e.g., the noise-orthogonal $\{A,E,T\}$ variables. These channels are synthesized by combining time-shifted one-way phase measurements recorded by each of LISA’s six laser links (Armstrong et al., 1999; Tinto & Dhurandhar, 2021), and are designed to suppress the otherwise dominant laser noise. We use the instrument model implemented in LISA GW Response (Bayle et al., 2022) to evaluate the interferometric phase shifts for each laser link given an incident GW signal. We use PyTDI (Staab et al., 2022) to convert these measurements into the various TDI channels. Finally, we inject the TDI signal into mock LISA data from the Spritz edition of the LISA Data Challenges (LDC-2b, Baghi, 2022). For simplicity, we use a “cleaned” version of the $1\,\,\mathrm{yr}$ dataset, which is free of artifacts such as noise glitches and data gaps but retains various instrumental and environmental background noises.

Fig. 2 shows $1\,\,\mathrm{yr}$ of simulated LISA observations for a representative $30\,M_{\odot}+20\,M_{\odot}$ BHB perturbed by a Sgr A^∗-like SMBH and undergoing eccentricity oscillations. The signal is characterized by a train of repeating, short-duration pulses lasting a few minutes each, which map onto the frequency domain as broadband bursts of excess power. The burst cadence is equal to the orbital period of the BHB, which is about $3\,\,\mathrm{days}$ for the example shown in Fig. 2. The majority of the GW power is radiated near periastron, causing the spectrum of each to burst to peak at a frequency of approximately (O’Leary et al., 2009)

where $(a,e)$ are the BHB semimajor axis and eccentricity at the time the burst was emitted, and $f_{\mathrm{orb}}(a)=(2\pi)^{-1}\sqrt{G(m_{1}+m_{2})/a^{3}}$ is the orbital frequency. The burst amplitude and frequency track the oscillating eccentricity, as shown in the lower right panel of Fig. 2. The bursts also experience complicated amplitude modulations attributed to both the source’s evolution and the motion of LISA, the latter of which depends on the source sky position and polarization. On account of the light travel time across LISA’s heliocentric orbit (i.e. Roemer delay), the burst arrival times measured by LISA deviate from the true orbital period by $\pm 500\,\,\mathrm{s}$, depending on the orbital phase of LISA.

Throughout Sec. 2–4, we assume the outer orbit of the triple is always coplanar with the sky (i.e. face-on/off inclination). For an inclined outer orbit, the GW signal will be Doppler shifted due to the projected motion of the inner binary along the line-of-sight as it orbits around the SMBH (Laeuger et al., 2024). This effect, known as the Roemer delay, will cause the time between bursts to oscillate about the orbital period, potentially biasing the estimate of the orbital period if the time delay represents a significant fraction of the inner orbital period. By assuming a coplanar outer orbit with the sky, we can ignore this Doppler shifting, which simplifies the data analysis problem. We provide a more detailed discussion of how the Roemer delay impacts our analysis in Appendix A. We do not consider relativistic corrections from the motion of the BHB (Robson et al., 2018; Xuan et al., 2023a) or the gravitational potential of the SMBH (Kuntz & Leyde, 2023; Sberna et al., 2022), opting to focus on a simple scenario in this proof-of-principle study and leave further generalization to future work.

Motivated by their transient and burst-like nature, we adapt the $Q$-transform (Brown, 1991) to detect highly eccentric BHBs in LISA data. Analysis pipelines based on the $Q$-transform search for arbitrary GW transients by filtering the data against windowed sinusoids (wavelets) in place of waveform templates (Chatterji et al., 2004; Chatterji, 2005), and have been applied extensively to process data from ground-based GW detectors and characterize transient noise (Davis et al., 2021; Acernese et al., 2023). We leverage the GWpy (Macleod et al., 2021) implementation of the $Q$-transform, which is itself derived from earlier burst search pipelines (Chatterji et al., 2004; Chatterji, 2005; Rollins, 2011; Robinet et al., 2020).

The $Q$-transform projects time series data onto overlapping time-frequency planes covered by an array of rectangular tiles, with each tile representing a wavelet. Each plane is parameterized by a quality factor, $Q\sim f/\Delta f$, which fixes the time-frequency aspect ratio of the individual tiles for that plane. The significance of any tile is quantified by its energy, given by (Chatterji et al., 2004; Chatterji, 2005)

where $(\tau,f,Q)$ are the time, frequency, and $Q$ of the tile, $x(t)$ is the whitened strain data, and $w(t,f,Q)$ is a window function (specifically, a bisquare window). We normalize the energy such that, in white noise, $|X(\tau,f,Q)|^{2}$ is exponentially distributed with unit mean and variance. The tile density is tuned to guarantee a fractional energy loss (mismatch) between adjacent tiles no larger than a desired amount, with smaller mismatches resulting in denser time-frequency tilings. The analysis outputs all tiles with SNR exceeding a pre-determined threshold, called triggers, where we define the trigger SNR as

Because a single burst typically produces multiple triggers, we group significant triggers that are overlapping or adjacent in time into clusters, and report the maximum-SNR trigger for each cluster. In our implementation, the $Q$-transform is carried out in two stages: the first stage performs a “quick pass” of the data using a low-resolution $Q$-tiling at $15\%$ mismatch; after clustering triggers with SNR greater than $8$, we search $5\,\,\mathrm{hr}$ of data around each trigger with a high-resolution tiling at $1\%$ mismatch to further refine the trigger parameters.

Prior to calculating the $Q$-transform, we whiten the data by its power spectral density (PSD), which is estimated empirically via a Welch median (Welch, 1967) with $2^{19}\,\,\mathrm{s}$ Hann-windowed and $50\%$ overlapping segments. To minimize over-whitening, we calculate the PSD on a “gated” version of the data (Zweizig & Riles, 2021; Usman et al., 2016). This gating is done by first dividing the data into short segments of duration $T_{\mathrm{gate}}=300\,\,\mathrm{s}$; if any of the strain values exceed some threshold, the data are multiplied by an inverse Tukey window, which is zero during the segment containing the excursion, and smoothly transitions from zero to one over a duration $T_{\mathrm{gate}}/4$ on either side. Through trial and error, we find that a gating threshold of $4.5$ times the root-mean-square strain (averaged over the full time series) is sufficient to remove the worst-offending bursts and prevent over-whitening.

Fig. 3 shows all triggers with SNR above $8$ after analyzing simulated LISA data with the $Q$-transform. For simplicity, we consider just a single TDI channel²²2Our analysis could be improved by using additional TDI channels to perform glitch rejection, as true GWs will propagate through LISA differently from instrumental glitches (Robson & Cornish, 2019)., TDI-$A$. We search over time-frequency planes with low $Q$-values³³3The lower bound of $Q=\sqrt{11}$ is an anti-aliasing condition., between $\sqrt{11}$ and $32$, specifically targeting bursts with short durations (maximizing timing accuracy) and high bandwidths. At $15\%$ ($1\%$) mismatch, this gives a total of $4$ ($14$) $Q$-planes. The searched frequency range is from $10^{-4}\,\,\mathrm{Hz}$ to $10^{-1}\,\,\mathrm{Hz}$, covering the nominal LISA observing band. The data contains three sources representing BHBs orbiting a SMBH at the Galactic Center, amidst background LISA noise from LDC-2b (Baghi, 2022). We also simulate a fourth source identical to the example from Fig. 2, but placed at a distance of $50\,\,\mathrm{kpc}$ (roughly the distance to the Large Magellanic Cloud) instead of $8\,\,\mathrm{kpc}$, showing how these more distant bursts are still individually resolvable at sufficiently large eccentricities. Three of the injected sources (red and green markers) achieve maximum eccentricity at the midpoint of the simulated observation period. We note that aligning the time of maximum eccentricity with the midpoint of the observation period is done merely to simulate an “optimal” detection scenario. Shifting the source in time will reduce the number of detected bursts. As an example, we include in Fig. 3 another source (blue markers) which achieves maximum eccentricity at the end of the observation period, and so the source is only observed while its eccentricity is increasing. Longer mission lifetimes not only increase the chances that LISA will see one of these putative sources during a high-eccentricity state, but will also improve LISA’s ability to study the source evolution over time and determine whether the dynamics are consistent with EKL-induced oscillations.

An important limitation of the $Q$-transform is that it is only sensitive to well-localized excess power, and thus cannot combine the power from multiple bursts far apart in time. Even if the summed SNR from all bursts is above our threshold, the source will not be detected if the SNR of individual bursts is consistently below the detection threshold. Our approach could thus be improved upon by stacking power from multiple bursts assuming some burst timing model (Romero-Shaw et al., 2023; Arredondo & Loutrel, 2021), allowing one to filter against a train of wavelets instead of isolated wavelets.

Even though the $Q$-transform only provides generic time-frequency information about the bursts, it is possible to make rough inferences about some of the source’s orbital parameters based on two key factors: the timing of the bursts, and the characteristic “peak” frequency where most of the excess power is concentrated. We outline here a simple method of estimating the orbital period and eccentricity evolution of BHBs using only the trigger time-frequency parameters.

A detection of binary eccentricity oscillations would offer conclusive evidence that the binary existed in a triple, providing key insight into the formation and evolution of stellar-mass compact binaries in dense cluster environments. In this work, we have studied the GW signature of stellar-mass BHBs in the Galactic Center that are driven to high eccentricities by a perturbing SMBH, and outlined techniques for detecting and identifying these sources in LISA data. BHBs in our own galaxy offer the best chances to detect EKL-driven dynamics, though we expect such systems to exist in other galaxies hosting SMBHs as well. We showed that an unmodeled time-frequency analysis can detect GW bursts emitted by highly eccentric BHBs in our galaxy, and characterize their time-frequency properties well enough to provide a rough estimate of the BHB orbital period and eccentricity evolution. This information can be passed to downstream analyses which use wavelets or waveform models to fit the signal, reducing the parameter space to be explored by these more computationally expensive techniques. Although we have focused on eccentric BHBs perturbed by a SMBH, our analysis techniques could be straightforwardly applied to generically eccentric compact binaries which do not undergo EKL oscillations. We worked with $1\,\,\mathrm{yr}$ of simulated LISA data for this study, but emphasize that longer observing times provide better opportunities to detect and study the long-term behaviour of EKL triples.

As this is a proof-of-concept exploration, there is much room to build upon and optimize various aspects of our methodology. Our burst search pipeline is based on the $Q$-transform, which uses bisquare-windowed sinusoids to maximize detection efficiency. However, this wavelet basis is overcomplete and cannot be used for direct signal reconstruction and subtraction, which will likely be a necessary step for the LISA global fit (Littenberg et al., 2020; Littenberg & Cornish, 2023). Future studies should explore alternative bases more suitable for burst subtraction, such as orthogonal Meyer wavelets, as done in the cWB pipeline (Klimenko & Mitselmakher, 2004; Klimenko et al., 2008; Drago et al., 2020), or Morlet-Gabor wavelets with a manually enforced orthogonality condition, as done in BayesWave (Cornish & Littenberg, 2015; Cornish et al., 2021). Another limitation of our method is that we search for bursts in isolation and do not stack power coherently from multiple bursts. Our sensitivity could potentially be enhanced by filtering on templates which chain together many bursts, where the timing between bursts and their relative amplitudes represent fit parameters.