Detecting Gravitational Wave Bursts from Stellar-mass Binaries in the mHz Band

The dynamical formation channels of gravitational wave (GW) sources typically involve a stage when the compact object binary source interacts with the environment, which may excite its eccentricity, yielding efficient GW emission. For the wide eccentric compact object binaries, the GW emission happens mostly near the pericenter passage, creating a unique, burst-like signature in the waveform. This work examines the possibility of stellar-mass bursting sources in the mHz band for future LISA detections. Because of their long lifetime (∼107 yr) and promising detectability, the number of mHz bursting sources can be large in the local Universe. For example, based on our estimates, there will be ∼3–45 bursting binary black holes in the Milky Way, with ∼102–104 bursts detected during the LISA mission. Moreover, we find that the number of bursting sources strongly depends on their formation history. If certain regions undergo active formation of compact object binaries in the recent few million years, there will be a significantly higher bursting source fraction. Thus, the detection of mHz GW bursts not only serves as a clue for distinguishing different formation channels, but also helps us understand the star formation history in different regions of the Milky Way.


INTRODUCTION
Gravitational wave bursts are expected to be a natural consequence of wide, eccentric compact object binaries' GW emission.In particular, for highly eccentric binaries, the GW emission happens mostly near the pericenter passage, creating a unique, pulse-like signature that lasts much shorter than the orbital period (see, e.g., Kocsis et al. 2006a;O'Leary et al. 2009;Gould 2011;Kocsis & Levin 2012;Seto 2013).Most of the burst emission does not yield a merger immediately; therefore, if the GW sources have an orbital period that is smaller than the observational time, the bursting signal will appear as repeated bursts (RB).
As the source of GW bursts, eccentric compact object binaries are often expected to form via dynamical channels, during which the binary's interaction with the environment excites its eccentricity e, reduces the pericenter Corresponding author: Zeyuan Xuan zeyuan.xuan@physics.ucla.edudistance, and causes effective GW radiation.Generally, several dynamical mechanisms involve a highly eccentric stage during the GW source's evolution.For example, in a hierarchical triple system (a tight binary orbiting a third body on a much wider "outer orbit"), the inner binary can undergo large eccentricity oscillations due to gravitational perturbations from the tertiary.This socalled eccentric Kozai-Lidov (EKL) mechanism (Kozai 1962;Lidov 1962;Naoz 2016) may potentially contribute to the overall merger rate of stellar mass compact objects at significant levels (e.g., Wen 2003;Hoang et al. 2018;Hamers 2018;Stephan et al. 2019;Bub & Petrovich 2020;Deme et al. 2020;Wang et al. 2021).Additionally, wide compact object binaries in the galactic field may interact with the surrounding environment, which can excite the binary's eccentricity (e.g., Michaely & Perets 2019, 2020;Michaely & Naoz 2022), resulting in merger events that are potentially observable.Moreover, a variety of dynamical mechanisms, such as GW capture, binary-single, and binary-binary scattering interaction, can take place in dense star clusters, producing compact object binaries with non-negligible eccentricity (e.g., O'Leary et al. 2009;Thompson 2011;Aarseth 2012;Kocsis & Levin 2012;Breivik et al. 2016;D'Orazio & Samsing 2018;Zevin et al. 2019;Samsing et al. 2019;Martinez et al. 2020;Antonini & Gieles 2020;Kremer et al. 2020;Winter-Granić et al. 2023).Furthermore, dynamical interactions in a flattened black hole distribution, such as in a stellar disk or an active galactic nucleus accretion disk, may also lead to highly eccentric mergers (Tagawa et al. 2021;Samsing et al. 2022;Muñoz et al. 2022;Gautham Bhaskar et al. 2023).
The detection of repeated bursts can greatly enhance our understanding of the GW sources' dynamical formation.In particular, with the corresponding data analysis methods (e.g., Tai et al. 2014;Loutrel & Yunes 2017;Bécsy et al. 2020), the LIGO-VIRGO-KAGRA (LVK) detectors may detect the residual eccentricity of highly eccentric sources, thus helping with identifying the fraction of GW sources formed in a variety of dynamical channels (e.g., East et al. 2013;Samsing et al. 2014;Coughlin et al. 2015;Gondán et al. 2018a,b;Zevin et al. 2021).Furthermore, as one of the main targets of the future space-based gravitational wave detectors, the extreme mass ratio inspirals (EMRIs) are expected to have high eccentricity in the mHz band, emitting detectable GW bursts at a cosmological distance (e.g., Glampedakis 2005;Hopman & Alexander 2006;Rubbo et al. 2006;Amaro-Seoane et al. 2007;Barack 2009;Berry & Gair 2013;Chen & Han 2018;Fan et al. 2022;Oliver et al. 2023;Naoz et al. 2022;Naoz & Haiman 2023).With the unique signature of highly-eccentric orbit, we can even enhance the peculiar acceleration measurement of the GW source by a factor of ∼ 100 (Xuan et al. 2023) and probe the gravitational potential surrounding these bursting systems (Zhang et al. 2021;Romero-Shaw et al. 2023).
In this work, we focus on stellar-mass compact object binaries that are bursting in the mHz band.These sources will typically have a semi-major axis a ∼ > 0.1 au, up to ∼ 10 5 au, with the pericenter distance a(1 − e) ∼ 10 −3 au, which makes the peak frequency of GW emission in the milli-hertz band.Although the GWs from these sources may be weaker because of their light mass (compared with EMRIs) and larger orbital separation (compared with bursting sources in the LIGO band, see, e.g., Randall et al. 2022), they are expected to have a much longer lifetime and therefore a larger number of systems in the mHz band (e.g., Fang et al. 2019).Furthermore, since the mHz bursting sources have not undergone significant orbital shrinkage, their eccentricity is more likely to be extreme, and their evolution is less rapid.These features can be used to extract relevant information about compact object binaries' surrounding environment, as well as their formation history.
This paper is organized as follows.We first specify the definition of bursting GW sources in Section 2.1, then quantify their detectability (Section 2.2) and lifetime (Section 2.3).In Section 3, we constrain the population of bursting sources from the observational results of LIGO (Section 3.2), then carry out numerical simulations to predict the number of detectable bursts for LISA.We show the results separately, for the field (Section 3.3), globular clusters (Section 3.4), and the galactic nucleus (Section 3.5) of the Milky Way.In Section 4, the properties of stellar-mass bursting sources are summarized and the implications are discussed.
Unless otherwise specified, we set G = c = 1.

Burst Definition
Heuristically, a burst is defined as a situation in which a significant amount of energy is emitted in a short amount of time compared to the orbital period.In particular, for eccentric binary sources, the pericenter time is usually defined with the orbital separation, velocity, and eccentricity, (r p , v p , e) as (e.g., O'Leary et al. 2009) where bin is the period of a binary with a mass m bin and semi-major axis a.Note that there is an order unity factor of (1 + e) −1/2 that we omit here.
As a proof of concept, consider a source with e > 0.9, we demonstrate below that, in this case, more than 88% of the GW energy emission is emitted in less than 1% of the orbital period during pericenter passage, regardless of its masses and semi-major axis.
The fraction of GW energy emitted during a burst can be estimated by integrating the power of GW emission, P , over the orbital true anomaly ψ.Particularly, P (ψ) is given by (Peters & Mathews 1963): in which m 1 , m 2 are the mass of the binary's components.
Throughout the paper, as a proof of concept, we adopt 1% as the fraction of burst time relative to the orbital period (note that other bursting fractions are straightforward to analyze.)In this case, the corresponding change in the orbital phase, δψ, can be solved by integrating the Keplerian motion of the binary near the pericenter passage, ψ ∼ 0, and require that T burst /T orb = 1%: 3 2 (1 + e cos ψ) 2 dψ , (3) where ψ is the time derivative of the orbital true anomaly: The fraction of energy emitted in the burst, compared with the total energy loss during one orbital period, is given by: where ⟨P ⟩ is the orbit-averaged GW energy emission power, given by Peters & Mathews (1963): We note that the fraction shown in Equation ( 5) is a function of the orbital eccentricity, as e will affect both the GW emission waveform (Eq.2) and the Keplerian motion, resulting in different δψ (Eq.3).Therefore, higher eccentricity values will yield a higher fraction of energy loss around the pericenter passage, making GW emission look more and more like bursts (see, e.g., Figure 1.The waveform in this case is calculated numerically using the x-model (see, e.g., eq.( 1)-( 13), (A1)-(A36) in Hinder et al. 2010)).
Adopting F(e) ∼ 90% and T burst /T orb = 1%, we can solve for the corresponding eccentricity threshold and find that a bursting source should have e > 0.9, regardless of its masses and semi-major axis.

Detectability of Bursting Sources in the Milli-hertz Band
The detection of GW bursts from highly eccentric sources is quite different from that of quasi-circular sources.In particular, when the binary's eccentricity is small (e.g., e ∼ 0.001 in Figure 1), the GW signal can be well approximated by a near-monochromatic, sinusoidal wave.Therefore, the GW templates of these sources are relatively straightforward to construct, which enables us to adopt the matched filtering method in the template fitting (Thorne et al. 1987;Finn 1992;Cutler & Flanagan 1994) and measure the source's parameters accurately.
Figure 1.The GW waveforms of a BBH system with the same orbital frequency but different eccentricities.We show a BBH system with m1 = m2 = 20 M⊙, orbital frequency f orb = 1.8 × 10 −5 Hz, luminosity distance D l = 8 kpc, and e = 0.001, 0.4, 0.9, respectively.As explained in the text, when the binary's orbital eccentricity increases, its GW energy emission will concentrate near each pericenter passage, turning the GW signal from a sinusoidal wave (e ∼ 0 case) into a "burst-like" waveform (e = 0.9 case).
On the other hand, when the source's eccentricity increases (e.g., e = 0.4, e = 0.9 in Figure 1), its GW emission will become stronger upon each pericenter passage, turning the signal into a burst-like waveform.For example, mHz bursting sources typically emit a bright GW signal during a period of 10 2 ∼ 10 3 s.However, for the rest of the orbital time (which is typically much longer than the signal time and can be days, months, or even years), the emission is suppressed.The transient nature of bursting sources makes them hard to detect, mostly because it adds to the computational expense of the template fitting and reduces their average signal-tonoise ratio.
Note that the transient nature of bursting sources can also lead to their long lifetime.Since the bursting sources spend a small amount of time emitting GW during each orbit (near the pericenter passage), the average power of orbital energy loss is small when compared with a circular binary with the same GW frequency as the eccentric bursting source's peak frequency.This means their orbit will shrink on a longer timescale.Therefore, bursting sources potentially have more significant numbers in future mHz GW detection.In the following sections, we will quantify how this feature enhances their detectability.
The signal-to-noise ratio (SNR) of a bursting source can be estimated analytically 1 .In particular, the energy loss of a bursting source mostly takes place near the pericenter passage; thus, consider a circular binary with the orbital radius the same as the bursting source's pericenter distance, r c = a(1 − e), then such a circular binary's GW emission power should be on the same order of magnitude as the "bursting power" of the bursting source.
This approximation is consistent with the nature of GW radiation, since the strain amplitude is proportional to the second-order time derivative of the binary's massquadrupole moment, while the peak GW frequency depends on the angular velocity at pericenter passage.Thus, a circular and eccentric binary with a similar radius/periapsis and velocity have a similar GW strain amplitude and frequency.
The bursting source has similar GW amplitude and frequency as the corresponding circular binary enclosed by the pericenter distance.The effective GW burst duration can be written as: Note that this expression also describes the period of the corresponding circular orbit.Moreover, highly eccentric sources emit GW emission in a wide range of frequencies, which peak at f peak ∼ f orb (1 − e) −3/2 (Peters & Mathews 1963).In the approximation that the GW burst emission from a highly-eccentric source has a similar frequency as the circular orbit with r c ∼ a(1 − e), the peak frequency is approximately 2 : in which f circ,GW is the GW frequency of the corresponding circular orbit, and T circ is the period of the circular orbit.
Similarly, the strain amplitude of a burst can be estimated by considering the circular orbit enclosed within the pericenter distance (Peters & Mathews 1963; Kocsis 1 Here we focus on giving a simple and direct estimation by comparing the system to a circular orbit with the radius of the pericenter distance.A more detailed analysis of the signal-tonoise ratio may be found in e.g.Kocsis et al. (2006b); O' Leary et al. (2009); Randall et al. (2022).
2 The peak frequency is often estimated as Leary et al. 2009).For consistency with other definitions in our treatment above, we adopt f circ,GW , which differs by an order of unity.
& Levin 2012), i.e., here D l is the luminosity distance of the binary, and the constant coefficient comes from the average of strain amplitude over the binary's orientation.Furthermore, m = (m 1 + m 2 )/2 is the average component mass, η s = 4m 1 m 2 /(m 1 +m 2 ) 2 = 4q/(1+q) 2 is unity for equal mass sources where q = m 1 /m 2 .
The analytical expression of signal-to-noise ratio (SNR), for a monochromatic source, is given by (see, e.g., Seto 2002): in which h A is the (time-domain) strain amplitude of the GW signal, √ f T is the number of observed cycles, and S n (f ) is the spectral noise density of LISA evaluated at GW frequency f (we adopt the LISA-N2A5 noise model, see, e.g., Klein et al. 2016a;Robson et al. 2019) 3 , f S n (f ) is the dimensionless spectral noise amplitude per logarithmic frequency bin.The numerator in Eq. (10) defines h c , the characteristic strain for a monochromatic source.
Note that T in Equation ( 10) stands for the time of observation when the signal is present, thus is different from the total observational time of LISA, T obs .For bursting sources, we sum over the pericenter passage time when they are emitting GW with significant amplitude (comparable with the corresponding circular source enclosed by the pericenter distance) to get T .For example, in a four-year LISA mission, if a bursting source bursts n times, then we have T = nT burst , while T obs = nT orb = 4 yr.
Plugging in Equation ( 7), (8), and (9) into Equation (10), then take into account the small amount of time when the source is bursting (T = nT burst = T obs /T orb × T burst if T obs /T orb = n ≥ 1, otherwise T = T burst if T obs ≤ T orb ), we can get an estimate of the bursting source's detectability: Note that for the first case (T obs ≥ T orb ) in Equation (11), the SNR can also be expressed as: ) which shows the dependency of (repeated) bursting sources' SNR on the eccentricity.
Equation ( 12) estimates the RB sources' SNR.Furthermore, it can be verified against the numerical results of eccentric binaries' detectability (e.g., Kocsis & Levin 2012;Hoang et al. 2019).For example, the SNR of a bursting BBH system with m 1 = m 2 = 10 M ⊙ , a = 1 au, e = 0.999 is estimated to be ∼ 200 in the Milky Way center, and ∼ 2 at the distance of the Andromeda.Therefore, Equation ( 12) implies that stellarmass mHz bursting sources can be detected in the Milky Way, and are marginally detectable in nearby galaxies.
For the second case (T obs ≤ T orb ) in Equation ( 11), the binary can only undergo one pericenter passage during the observation.Therefore, the signal-to-noise ratio represents a single GW burst's detectability.We show in Figure 2 that such SNR of a single burst can be expressed as a function of the binary's pericenter distance, and for stellar-mass binaries, the distance of single burst detection is mostly limited to the Milky Way.The SNR peaks at a burst frequency of about 6 mHz, the minimum of the LISA sensitivity curve, corresponding to a pericenter distance of ∼ 10 −3 au.
We emphasize that, the expression of SNR in Equation (11) are based on the average GW emission power, and thus may not describe the full detectability of bursting sources.For example, Figure 3 shows the comparison of time and frequency domain waveforms of a bursting BBH and a circular double white dwarf (DWD), calculated numerically using the x-model (Hinder et al. 2010).4For a comprehensive comparison, we adjust Figure 2. The signal-to-noise ratio of a single GW burst, as a function of the source's pericenter distance.We show a BBH system with m1 = m2 = 10 M⊙, placed at D l = 1, 8, 50, 765 kpc, respectively.The SNR of a single GW burst is plotted as a function of the binary's pericenter distance a(1 − e), and the corresponding GW burst frequency is plotted on the top x-axis.We note that repeated burst sources can have multiple bursts detected during the observation, and the overall SNR can be higher than the single burst case of this figure.See Figure 4 for the repeated burst case.
their parameters to have the same SNR ∼ 12.As shown in the figure, even if the bursting BBH system has the same SNR as the circular DWDs, its time-domain GW amplitude is thousands of times greater than the latter one.In other words, with the same SNR, the timedomain strain amplitude of a mHz "burst" is much greater than that of a continuous GW.Therefore, we may use this signature to enhance the bursting sources' detectability (e.g., East et al. 2013;Tai et al. 2014;Loutrel 2020;Wu et al. 2023).
We note, however, that the transient nature of bursting sources also adds to the difficulty of data analysis.In particular, the typical strategy for detecting GWs uses matched filtering (Finn 1992;Cutler & Flanagan 1994), which heavily relies on the accurate model of GW signal and is not well-suited for searching for discrete bursts localized in time and frequency (see, e.g., Tai et al. 2014;Loutrel & Yunes 2017;Loutrel 2020).In the context of LIGO data analysis, extensive efforts were made to identify transient events using time-frequency methods, such as power stacking (East et al. 2013), the TFCLUS-TER algorithm (Sylvestre 2002), wavelet decomposition (Klimenko & Mitselmakher 2004), and the Q-transform (Bassetti et al. 2005;Tai et al. 2014).These methods may lead to a lower SNR relative to the matched filtering but are robust to different kinds of transients and can combine the information from multiple bursts from a single source.For LISA data analysis, burst detection methods for stellar mass binaries are still underdeveloped, but some efforts have been made to detect highly eccentric EMRIs (see, e.g., Barack & Cutler 2004;Cornish & Larson 2003;Hopman et al. 2007;Porter & Sesana 2010;Mikóczi et al. 2012) and burst from scattering of black holes (Kocsis & Levin 2012).
Moreover, the detection of bursts with the LIGO-VIRGO-KAGRA network has the advantage that the source may be more precisely localized in the sky, given the ability to measure the GW arrival time difference at different detector sites.In contrast, for persistent inspiraling sources in the LISA band, their localization depends on the modulation of the signal caused by the annual motion of the LISA satellite constellation along its orbit.While "time delay interferometry" (see, e.g., Tinto & Larson 2004;Tinto et al. 2021) may be utilized to identify unmodeled astrophysical GW sources with LISA, corresponding mock data analysis methods are currently unavailable for cases with only a small number of bursts in the observational period.Thus, it is currently not known how accurately the source may be localized in the sky, if at all in this case.Additionally, distance information may also be unavailable, as it degenerates with the total mass of the source and the separation at the closest approach.Thus, confusion of several distinct binaries with repeating single GW bursts may be a major uncertainty for these types of sources.It remains unclear how many repeated bursting sources (with what SNR) are needed for secure detection.For conservative purposes, in this paper, we adopt the average SNR as the criteria for evaluating the bursting sources' detectability.This approach represents the power of the GW signal and provides an optimistic estimate of bursting sources' detectability in the LISA band with current data analysis methods.

Lifetime of Bursting Sources in the Milli-Hertz Band
Most of the stellar-mass bursting sources start their evolution from a wide configuration.These systems undergo some dynamical interaction, such as EKL, fly-by, or a strong encounter, resulting in a highly eccentric configuration.At this stage, the system undergoes repeated GW bursts.Following this RB stage, the binary's orbit shrinks, and if GW emission dominates the evolution, it becomes circularized, yielding a merged system (e.g., see Figure 8).Particularly, an RB source's lifetime, τ RB , can be estimated using the merger timescale for binaries with extreme eccentricity (Peters 1964): Equation ( 13) describes the isolated evolution of bursting binaries.However, some bursting sources that undergo dynamical interaction can have eccentricity oscillation throughout their evolution (e.g., EKL merger systems Hoang et al. 2019;Randall & Xianyu 2019;Deme et al. 2020;Emami & Loeb 2020;Chandramouli & Yunes 2022).In this case, we adopt the maximum eccentricity during the oscillation, e max , as the value of e in Equation ( 13), which yields a lower limit of the timescale for a source to burst.
We are interested in comparing the length of a bursting system's RB stage with its latter stage of evolution (inspiral with moderate eccentricity).Thus, for simplicity, we assume that the binary circularizes and shrinks to a radius r c ∼ a(1 − e) after the RB evolution.Such systems are well described in EKL systems or captured systems (e.g., see Kocsis & Levin 2012;Naoz 2016).Further, this approximation is consistent with the concept that the pericenter of a highly eccentric GW source remains nearly the same (in the absence of dynamical evolution, see for details Peters & Mathews 1963;Peters 1964).
Thus, under this notion, when the bursting binary's eccentricity drops to a moderate value, its remaining inspiral time, τ inspiral , can be estimated as (Peters 1964): We note that Equation ( 14) also serves as a general estimate of millihertz circular BBHs' merger timescale (if we replace f burst in the equation with f circ,GW ).Therefore, the equation implies that a stellar mass circular binary will stay in the mHz band, in particular f circ,GW ≳ 1mHz, for ∼ 10 3 − 10 5 yr (e.g., Chen et al. 2020).
As can be seen from Equation ( 13) and ( 14), the timescale for a highly-eccentric source's RB stage is longer than the following inspiral stage with moderate eccentricity: which reflects the fact that in the RB stage, the eccentric binary will spend most of its time at large separation, yielding a lower GW emission and orbital energy loss compared to a circular orbit with a separation a(1 − e).
For example, the EKL mergers of BBHs in the galactic nucleus are expected to have e ≳ 0.999, induced by the supermassive black (SMBH) via the EK mechanism (Naoz 2016;Hoang et al. 2018;Hoang et al. 2019).Thus, Equation ( 15) means the RB stage of EKL mergers will be hundreds of times longer than their inspiral time with moderate eccentricity and the same a(1 − e).
Figure 4 shows the lifetime (Eq.15) and detectability of RB sources with a more detailed calculation for the SNR than in Equation ( 11 We also over-plot, in Figure 4, the merger timescales for isolated binaries in each panel (thin solid lines in the figure).For comparison purposes, in the left panel, we further show the direct merger limit for BBHs (black solid line at the bottom).In the right panel, we show the Roche limit for DWD in the black dash-dotted line (which is the limit that DWDs will undergo mass transfer).
Overplotted in Figure 4 is a representative example.Specifically, starting from a given initial condition, the RB source will evolve in the a−e parameter space following the changes in its semi-major axis and eccentricity due to the GW emission (Peters 1964 Using these equations, we plot the evolutionary track (the solid lines with arrows) on each panel in the absence of further dynamical evolution.Consider, for example, the evolutionary track in the right panel.This DWD system began in the undetectable regime; however, as it evolves, it becomes detectable via RB, with SNRs corresponding to the color on the map.When the system crosses the lifetime line, it will merge within that time, thus undergoing repeated bursts for that remaining time.In particular, the DWD example crosses the 10 6 yr mark with a ∼ 0.03 au and e ∼ 0.99, which means it will survive for an additional 10 6 years as a detectable repeated burst source.
As can be seen in Figure 4, the lifetime of a detectable RB source in the Milky Way can be as much as ∼ 10 7 yr, since the edge of the detectable region Here we show a BBH system (Left Panel) with m1 = m2 = 10 M⊙, and a DWD system (Right Panel) with m1 = m2 = 0.75 M⊙, placed at D l = 8 kpc and observed for 4 yr with LISA.The color represents the compact binary's signal-to-noise ratio as a function of its semi-major axis a and eccentricity 1 − e.As a proof of concept, we depict an example system's evolution (purple line with star and arrows) for both parameter spaces and plot the corresponding orbital time on the top x-axis.The solid lines with color represent the merger timescale τmerger for a given (a, e) configuration of the binary.In other words, when the evolution track of the binary gets into the colored region and gets across a given merger timescale, we can estimate the remaining time for it to emit detectable GW bursts before the merger.We note that the signal-to-noise ratio is suppressed in the region below the dashed lines (tmerger < T obs ) because the binary will merge during the observation, thus having less observational time than the total LISA mission time.Therefore, when calculating the SNR for these binaries following Equation.(B11), we replaced the T obs with tmerger.
We note that in Figure 4, another population of BBHs potentially exists, which is characterized by moderate eccentricity (e ∼ 0) and long lifetime (t merger ∼ 10 7 yr).However, these systems exhibit a much lower GW frequency (∼ 0.1mHz) compared to the bursting sources considered in this paper (≳ 1mHz).For example, a BBH system with e ∼ 0, a ∼ 0.03 au would be detectable in the Milky Way (∼ 10kpc) for 10 7 yr, with an orbital period of T orb ∼ 10 4 s and GW frequency f circ,GW ∼ 10 −4 Hz.Note that because we expect a large population of circular DWDs in the sub-millihertz band (see, e.g., Nissanke et al. 2012;Lamberts et al. 2018;Xuan et al. 2021), the identification of the near-circular, low-frequency BBHs is beyond the scope of this study.Thus, in this paper, we focus on the dynamically formed, highly-eccentric bursting sources, with f GW ≳ 1mHz.
For comparison purposes, we show in Figure 5 the maximum detectable distance of RB sources, assuming a 4 yr LISA mission and a detection threshold of SNR = 5.As shown in the figure, the detection of bursting BBHs can be promising within the distance of nearby galaxies, while the detection of bursting DWDs is limited to the Milky Way.This result is consistent with the estimation in Section 2.2 (see, Equations ( 11), ( 12) and Figure 2).
One important caveat is that the estimation of SNR in Figure 4 and 5 is based on the averaged power of the GW signal, where for a binary with T orb larger than T obs , the number of bursts we expect to detect during the observation is smaller than one.In this case, the average SNR in Figure 4 is smaller than the SNR of detected bursts because we detect only a fraction of such bursting sources, only having a probability of less than one for detection.However, for these non-repeated bursting sources, a single burst falling within the observational window will be brighter than shown in the Figure according to Eq. ( 12), and hence can still be bright enough to detect with LISA.Moreover, in the region below the dashed black line of Figure 4, the merger timescale of a GW source is shorter than the observational duration, which means the fast-merging sources in this region will not be observed for the whole LISA mission, thus having lower SNR.Therefore, we can further include these facts and finally get: S n (f burst ) .
(17) in which τ insp is the binary's remaining inspiral time.
Equation ( 17) gives the full detectability of (both repeated and non-repeated) burst sources, providing that at least one burst takes place during the observation.The suppression of SNR caused by τ insp is taken into account in the numerical results of this paper.

General Considerations
As shown in Section.2, the RB stage of a stellar-mass GW source is well detectable within ∼ < 10 kpc, and lasts longer than the inspiral stage with moderate eccentricity.Therefore, the population of bursting sources can have a significant contribution to the detection of dynamically formed GW sources, especially in the Milky Way.
As a proof of concept, we heuristically estimate the expected populations of bursting BBHs.Due to the large uncertainties of star formation and detailed dynamical evolution, we focus on the expected order of magnitude of the bursts.We consider three regimes that are expected to have eccentric sources, namely: Galactic Field, Globular Clusters (GCs), and Galactic Nucleus (GN), and show the results in Figure 6.For comparison purposes, we also present the number estimation of GW bursts from the EMRIs in the universe.
For other kinds of bursting sources, such as DWDs, their number, parameter distribution, and evolution is less understood, mostly due to the uncertainty in the formation scenario and the complexity of tidal effects in extreme eccentricity cases.Therefore, we will leave these populations for future work.
Before providing the estimation for bursting sources in the three regions mentioned above, we first describe a general estimate agnostic to the astrophysical formation channel of the source.

Steady-State Approximation and Constraints from LVK
Motivated by LVK observations, we first estimate the number of RB sources under the assumption that LVK represents a steady-state population.In particular, assume there is a continuous birth and death of compact object binaries in the Milky Way, keeping the total number of systems unchanged.In this case, the replenishing rate of a GW source is balanced by its merger rate.In other words, Γ rep = Γ merger .
For GW sources formed through dynamical interaction, the expected number of systems in the RB stage, N RB , is found by multiplying the average lifetime of the RB stage (Eq.( 13)) with the formation rate of the source: which is also applicable to other evolutionary stages, such as the inspiral with moderate eccentricity.
The LVK observations indicate that BBHs merge at a rate of 15-38 Gpc -3 yr −1 (Abbott et al. 2021).Therefore, we can estimate the BBHs merger rate per galaxy by dividing the galaxy number density.Adopting the galaxy number density of 0.02 Mpc −3 (Conselice et al. 2005), the merger rate of BBHs in an averaged galaxy is estimated as: Γ merger,avg ∼ 10 −6 yr −1 . (20) Thus, plugging Equation (20) into Equation ( 19) places a constraint on the number of bursting BBHs we can detect locally.For example, bursting GW sources can spend up to a few 10 7 yr in the detectable regime (Figure 4).Therefore, the upper bound of RB BBHs' number, from the steady-state approximation and LIGO's constraint, is on the order of N RB ∼ 10 7 yr 10 −6 yr −1 ∼ 10 sources in the Milky Way.Below we use this approach to place constraints on different formation channels.We emphasize that this constraint only serves as a rough estimation of bursting BBHs' number.It can be well-exceeded if taking into account the non-equilibrium  ties.Specifically, we show the number of detectable bursting sources from each channel (x-axis) and their contribution to the number of bursts detected per year (y-axis).For details of the simulation, see Appendix A. Note that the number of detectable bursting sources is roughly proportional to the observational time.The solid points represent the population with convincing evidence of existence, while the hollow points are for the potentially existing population of bursting sources.We emphasize that if future detection does not show evidence of the population marked in hollow points, it serves as a strong constraint on the corresponding formation channels.
formation history of compact binaries.For example, different regions in a galaxy may undergo starbursts while others have a quiescence phase, making the number of bursting sources fluctuate in a wide range.(See Section 3.5 for a detailed analysis.)Further, assuming that the observed rate of mergers all form at mHz or smaller frequencies5 due to their long lifetime, bursting binaries could be the dominant source in the local mHz BBHs population (see Equation (19)).For example, consider the lifetime of circular (non-bursting) BBHs in the mHz band, τ inspiral ∼ 10 3 − 10 5 yr (see Equation ( 14)).Provided that both bursting and non-bursting sources are detectable within the Milky Way, the shorter lifetime of circular sources results in their small number in the local population, N circ,BBHs ≲ 10 5 yr 10 −6 yr −1 ∼ 0.1 in the mHz band, while bursting sources represent the majority of mHz sources (N RB ∼ 10).
Below, we consider three characteristic BBHs formation channels and provide heuristic estimations for their bursting properties.

BBHs in the Field
Isolated stellar binaries in the galactic field have been proposed as a possible channel explaining the LVK observations (e.g., de Mink & Belczynski 2015;Dominik et al. 2015;Belczynski et al. 2016;Eldridge et al. 2017;Giacobbo et al. 2018;Olejak et al. 2020).However, in the context of bursting systems, we focus on wide bina-ries in the field, which are less likely to be seen with LVK without external perturbations.Providing that external perturbations such as flybys or galactic tides will affect their evolution, these wide-field binaries can be driven to extreme eccentricity, emit bursting GW signals, and become mergers in the end.
We adopt the model proposed in Michaely & Naoz (2022), as well as the steady-state approximation, to calculate the eccentricity and semi-major axis distribution for the BBH mergers in the disk of the Milky Way.In particular, following Michaely & Perets (2019), the merger rate from the wide-binaries field channel can be within the range of 5 +5 −3 Gpc −3 yr −1 in the local universe (up to 50Gpc −3 yr −1 if taking elliptical galaxies into account, see Michaely et al. in prep).Based on this rate estimate, assuming a optimistic LISA observation timescale T LISA,obs = 10 yr and detection threshold SNR= 5, we find that the minimum (maximum) number of detectable bursting BBHs in the MW, induced by fly-by interaction, is 0.7 (3.3), and they are expected to emit 12 (61) GW bursts per year.In these bursting sources, we identified the lower (upper) bound of the repeated burst sources number as 0.4 (1.8), and 0.3 (1.5) for non-repeated burst sources.However, because of the long orbital period of non-repeated burst sources, the chance that we will observe GW bursts during their pericenter passage is small.Thus, the number expectation of non-repeated bursts is negligible (∼ 0.1).This is consistent with the estimate in Kocsis et al. (2006a).See Appendix A.1 for more details.
Here, we adopt the eccentricity distribution of BBHs in Martinez et al. (2020, see their figure 4) and the spatial distribution of GCs in the Milky Way from Arakelyan et al. (2018).Since we are interested in the RB stage where BBHs are in a wide configuration with low GW frequency (∼ mHz), we evolve the systems shown in Martinez et al. (2020), which is in a higher frequency band, back to the former RB stage.The bursting time of these sources is calculated by counting the time difference between the point when they become detectable and the point when their eccentricity drops below 0.9 (see the example tracks in Figure 4).In the simulation, we use steady-state approximation to calculate the expectation of bursting sources number (For detailed information, see Appendix A.2).
We adopt the BBH merger rate in GCs to be: R 0 = 7.2 +21.5 −5.5 Gpc −3 yr −1 (e.g., Rodriguez et al. 2016;Kremer et al. 2020;Antonini & Gieles 2020).Assuming T LISA,obs = 10 yr and SNR=5, we find that the minimum (maximum) number of detectable bursting sources in the Milky Way GCs is 1.2 (20.1), which corresponds to a number of GW bursts 15.2 (253) per year.Among these sources, the minimum (maximum) number of repeated burst sources is 0.34 (5.7), while the minimum (maximum) number of non-repeated burst sources is 0.85 (14.4).

BBHs in the Galactic Nucleus
The Milky Way's galactic nucleus offers a natural place for the formation of bursting BBHs (see, e.g., Kocsis & Levin 2012;Hoang et al. 2019;Stephan et al. 2019;Arca Sedda et al. 2023;Zhang & Chen 2024).In particular, BBHs orbiting around the supermassive black hole in the galactic nucleus will undergo eccentricity excitation via the EKL mechanism, resulting in the bursting signatures on their GW signal.
The orbital evolution of BBHs in the GN is quite different from the isolated binaries (e.g., Equation ( 16)) since they are strongly affected by the gravitational perturbation from the SMBH tertiary.Thus, to get the bursting properties of these sources, we carried our detailed simulations of hierarchical triple systems, including the secular equations up to the octupole level of approximation (Naoz et al. 2013a), general relativity precession (e.g., Naoz et al. 2013b), andGW emission (Peters 1964;Zwick et al. 2020).These calculations follow a similar approach to Hoang et al. (2018).
We first present the results under the steady-state approximation (see details in Appendix A.3).This approximation can well-describe the main population of stars in the GN (old population, aged 2 ∼ 8 Gyr Chen et al. 2023).Assuming a 10 yr LISA observation with the signal-to-noise ratio threshold 5, we got the expectation of EKL induced, bursting BBHs number ∼ 1 in the mHz band, with the number of detectable bursts ∼ 100 per year.All the bursting sources in the simulation are repeated burst sources, which is consistent with the fact that wide binaries will evaporate quickly in the active dynamical environment of GN.The result is calibrated using the m − σ relation (Merritt 1999;Kormendy & Ho 2013) and observational results of the Milky Way center.
However, unlike the previous channels (i.e., in the field and GCs), we go beyond the steady state approximation because observations suggest a recent (2 ∼ 8 Myr) star formation occurred in the GN (Paumard et al. 2006;Lu et al. 2008;Lu et al. 2009;Do et al. 2013;Chen et al. 2023).Specifically, this recent star formation may have formed a young nuclear star cluster (YNC) within 0.5 pc from the central SMBH, with a total mass of ∼ (1.4 − 3.7) × 10 4 M ⊙ and a top-heavy mass distribution (Lu et al. 2009).Therefore, it can hold ∼ 100 − 400 BHs as a result of the stellar evolution.Furthermore, since stars often reside in binaries, and high-mass stars reside in higher multiples (Pribulla & Rucinski 2006;Tokovinin et al. 2008;Raghavan et al. 2010;Sana & Evans 2011;Sana et al. 2012;Moe & Di Stefano 2017), we expect these black holes to form ∼ 100 BBH systems in the YNC.This expectation is supported by myriad observational and theoretical arguments (e.g., Ott et al. 1999;Martins et al. 2006;Rafelski et al. 2007;Pfuhl et al. 2014;Alexander & Pfuhl 2014;Naoz et al. 2018;Gautam et al. 2019;Chu et al. 2023).For the newly- born BBHs, we carried out Monte-Carlo simulations to get the evolution of RB sources' number as a function of the age of YNC (see Figure 7), as well as the number of bursts in the Appendix (see Figure 9).As is shown by the simulation, although the YNC has a relatively small mass (M ∼ 10 4 M ⊙ ), there can be 2 ∼ 4 RB sources with 60 ∼ 150 bursts detected per year.These numbers are larger than the expectation of bursting systems from the old population of stars under the steady-state approximation, even if the latter one has a much larger total mass of stars (∼ 10 7 M ⊙ ).
Combining the simulation results of bursting BBHs in the old population of stars with those in the YNC, we get the number estimation of detectable bursting BBHs in the inner 1000 au ∼ 0.5 pc of the galactic center as 2 ∼ 6.These sources will contribute to 160 ∼ 254 bursts per year.
Furthermore, observations suggest that inwards of S0-2's orbit ( ∼ < 1000 au) there is a hidden mass of ∼ 3000 M ⊙ (e.g., Do et al. 2019;GRAVITY Collaboration et al. 2020).This is supported by theoretical arguments of the stability of S0-2, (e.g., Naoz et al. 2020;Zhang et al. 2023;Will et al. 2023).We, thus, explored the possibility that these unidentified objects are all stellar-mass BHs binaries formed within the recent 10 7 yr, and find that there can be ∼ 3 − 20 detectable RB sources with ∼ 10 − 300 bursts per year under this assumption (see the grey data point in Figure 6).We note that the existence of BBHs in the inner 1000 au is highly uncertain, but any potential population in this region can have a significantly high fraction of bursting systems.Therefore, GW burst detection in the future will serve as a strong constraint on the number of BBHs in the inner 1000 au of the GN.For more details, see Appendix A.3.
These results highlight the close relationship between the number of bursting BBHs and their formation history.In the young population of BBHs, RB systems will remain detectable up to ∼ 10 7 yr after being driven to extreme eccentricity by the SMBH, which is longer than the YNC's age.However, despite the large total mass of the old population, they are likely to reach a steady state because of the low formation rate of BBHs.Therefore the young population will appear to have a much higher RB source fraction than the average value, as depicted in Figure 7.
In other words, RB sources can serve as a tracer of active star formation and dynamical evolution of compact object binaries in the recent ∼ 10 7 yr history of the central Milky Way.We speculate that excess detection of bursting sources in certain regions of the Milky Way may indicate a recent episode of compact object formation.
In particular, Naoz & Haiman (2023) estimate that there can be as much as 100 ∼ 2000 EMRIs detected per year if there is a significant fraction of galaxies holding an SMBH binary at their galactic center.Their simulation shows that most of the detectable EMRIs via SMBH binary channel are in a highly eccentric configuration, with the orbital period in the range of ∼ 1−10 yr.Therefore, we estimate the number of detectable bursting sources and bursts per year by multiplying the EM-RIs number with their orbital frequency (∼ 0.1−1 yr −1 ).It turns out that under the assumption of Naoz & Haiman (2023), there can be ∼ 1000 − 20000 bursting EMRIs detected during the LISA mission, which contribute to ∼ 100 − 20000 GW bursts per year.
We emphasize that, similar to the case of bursting BBHs in the inner 1000 au of the Milky Way (see Section 3.5), the existence of SMBH binaries in the universe is poorly constrained.Therefore, in Figure 6, we use a hollow diamond to highlight the uncertainty of EMRI bursting sources formed via the SMBH binary channel.Similar to the arguments in Naoz & Haiman (2023), whether or not we will detect a large number of bursting sources as expected by this channel can strongly constrain the fraction of SMBH binaries in the universe.

DISCUSSION
Many dynamically formed GW sources are expected to undergo an evolutionary stage at which the compact object binary is driven to extreme eccentricity.Such a system will emit GW bursts effectively upon each pericenter passage (see, for example, Figure 1), creating a pulse-like pattern in the signal before finally losing enough orbital energy and becoming a merger.Therefore, GW bursts from highly eccentric compact object binaries can be a promising tracer of their dynamical formation.
In Section 2, we estimate the detectability and lifetime of these sources analytically (see Eqs. ( 13)-( 17)) and compare the results to numerical calculations (see Figure 4).Particularly, we show that stellar-mass bursting sources should be detectable in the Milky Way (∼ 10kpc), with a much longer lifetime than other mHz sources (up to 10 7 yr).Moreover, a bursting source yields a larger strain amplitude than a low eccentricity source with the same average signal-to-noise ratio (see Figure 3).Considering this feature, we can potentially enhance these sources' detectability and parameter extraction in future data analysis.
Notably, BBH systems in many different regimes, such as in the field, globular clusters, and galactic nuclei, can naturally form bursting GW sources in the mHz band.For example, in Section 3, we show that bursting sources can dominate the population of mHz BBHs in the Milky Way since their longer lifetime (Equation ( 15)) results in a larger number expectation (Equation ( 19)).Adopting the constraints from LIGO's observation, we estimate ∼ 10 bursting BBHs detectable for the future LISA mission under the steady-state approximation (see Section 3.2).
In Figure 6, we present the estimated number of bursts as a function of bursting sources for different formation environments.In particular, we consider stellar-mass BHs in the galactic field, the globular clusters, and the galactic nucleus.We find that these channels can contribute to the number of detectable stellar-mass bursting BBHs in a range of 3 ∼ 45, with 10 2 ∼ 10 4 mHz GW bursts observed during the LISA mission.
We highlight that the number of detectable bursting sources can well exceed the steady-state value since the formation history of bursting sources significantly affects their number in the mHz band.In particular, the Milky Way is expected to be a star-forming galaxy.Therefore, in Section 3.5, we go beyond the steady-state approximations and simulate the time evolution of bursting sources in the galactic center Young Nuclear Star Cluster, which is expected to form ∼ 2−8 Myr ago.The simulation result supports that long-living bursting sources can remain detectable for ∼ 10 7 yr.This region, with an active star formation in the recent few million years, will have a much higher fraction of bursting sources (2 ∼ 4 bursting sources out of ∼ 100 BBHs) than those old ones (∼ 1 bursting sources out of ∼ 1500 BBHs).In future observation, the distribution of bursting sources, as well as their actual number, can serve as a valuable tool to probe the GW sources' creation time in different regions of the Milky Way.
To conclude, unlike other mHz sources with moderate eccentricity, the detection of stellar-mass bursting systems is mostly limited in the Milky Way and nearby galaxies (see, e.g., Figure 5).However, these sources' bursting nature also results in a long lifetime and a potentially large number in the entire mHz population.Therefore, it is very likely that we can find many bursting sources at a close distance (e.g., a few dozen bursting BBHs in the Milky Way), with hundreds to thousands of GW bursts detected in the LISA observation.By studying the properties of these sources, we can better understand the early stages of dynamical formation and map the compact objects' formation history in our galaxy.
The authors thank Cheyanne Shariat, Sanaea Rose, Brenna Mockler, Xian Chen, and Xiao Fang for their valuable discussions and the anonymous referee for their useful comments.ZX acknowledges partial support from the Bhaumik Institute for Theoretical Physics summer fellowship.ZX, SN, and EM acknowledge the partial support from NASA ATP 80NSSC20K0505 and from NSF-AST 2206428 grant and thank Howard and Astrid Preston for their generous support.For the BBHs in the galactic field, flyby gravitational interactions with other neighbors may excite their high eccentricity, driving the binary into a merger (e.g., Michaely & Perets 2019;Michaely & Naoz 2022).In this work, we adopt the model in Michaely & Naoz (2022) to calculate the properties of bursting BBHs born from this channel.

APPENDIX
In particular, we choose the configuration of a Milky-Way type galaxy, with the stellar density profile and wide BBHs fraction following Eq.( 23)∼( 26) in Michaely & Perets (2019).For simplicity, we further assume that all the BBHs have the mass of 10 − 10 M ⊙ , with the log-uniform distribution in the semi-major axis (from 100 to 50000 au).Second, we run a Monte-Carlo simulation, randomly choose the position and semi-major axis of BBHs following the distribution mentioned above, and calculate the fly-by merger rate as a function of semi-major axis (see Eq. ( 16)∼( 23) in Michaely & Naoz 2022).
After undergoing a fly-by interaction, the BBHs with a given semi-major axis need to have the eccentricity exceed a threshold, e crit , to make themselves merge before the next fly-by (see Eq. ( 2) in Michaely & Naoz 2022).Therefore, we randomly generate the eccentricity of fly-by induced merger in the range of e crit to 1, following the thermal distribution F (e) = 2e.
Using the methodology mentioned above, we can specify each merger system's semi-major axis and eccentricity in the simulation.This information allows us to calculate how many bursts they emit, and how long they stay detectable.For example, we can evolve a system following Eq.( 16) and trace its signal-to-noise ratio using Eq. ( 17).Once the system becomes detectable, the merger timescale and number of bursts can be calculated via integration (see the example track in Figure .4 for demonstration).The number of bursting systems in the fly-by channel is calculated using the steady-state approximation.In other words, we multiply the corresponding merger rate with the average detectable time for each semi-major axis, then do the summation over different semi-major axis to get the expectation of sources' number (see Eq. ( 19)).
Assuming a 10 yr LISA observation with the signal-to-noise ratio threshold 5, we got the expectation of fly-by induced, bursting field BBHs number is N field ∼ 3.3 in the mHz band, with the number of detectable bursts N burst,field ∼ 61 per year and merger rate Γ field ∼ 5 × 10 −7 yr −1 in the Milky Way (or ∼ 10 Gpc −3 yr −1 in the universe).Among these sources, we identified ∼ 1.8 repeated burst source and ∼ 1.5 non-repeated burst source (i.e., T orb > 10 yr).
According to Michaely & Perets (2019), the BBH merger rate from the field channel can be within the range of 5 +5 −3 Gpc −3 yr −1 in the local universe.Moreover, the number of bursts and bursting sources is proportional to the merger rate of BBHs.Therefore, to get a realistic estimation of the bursting sources' properties, the numbers we got from the simulation should be multiplied by a factor of 0.5 +0.5 −0.2 , which is reflected in the results in Section 3.3.

A.2. BBHs in Globular Clusters
As is shown in Section 3.4, we use the eccentricity distribution of BBHs at a given frequency (see Fig. 4 in Martinez et al. 2020) and the spatial distribution of globular clusters in the Milky Way (Arakelyan et al. 2018) to generate the initial condition of bursting BBHs after undergoing dynamical interaction in GCs.For conservation purposes, we exclude the single-single capture and few body capture channels since the GW sources from these two channels are mostly formed above the mHz band, thus may not be suitable targets for LISA (Kocsis 2020).
Adopting the same method as is shown in Section A.1, we carried out Monte-Carlo Simulations and generated the parameters of bursting BBHs following the distribution mentioned above (Samsing & D'Orazio 2019;D'Orazio & Samsing 2018;Martinez et al. 2020), then evolved them to the merger.The number of bursting sources and bursts is calculated under the steady-state approximation, taking into account the burst sources' properties from different channels of merger in GCs.Particularly, for a given sub-channel in GCs, we get the averaged bursting lifetime and the total number of bursts from the simulation, then multiply it with the merger rate (see Eq. ( 19)).
Assuming a 10 yr LISA observation with the signal-to-noise ratio threshold 5 and GCs merger rate Γ GC ∼ 1 × 10 −6 yr −1 in the Milky Way ( ∼ 20 Gpc −3 yr −1 in the universe), we find the number of detectable bursting BBHs in the GCs of Milky Way to be N GC ∼ 14, with the number of detectable bursts N burst,GC ∼ 176 per year.Among these sources, we identified ∼ 4 repeated burst source and ∼ 10 non-repeated burst source.Similar to the case of field binaries in Section 3.3, the number of bursting sources is proportional to the merger rate (R 0 = 7.2 +21.5 −5.5 Gpc −3 yr −1 (Rodriguez et al. 2016;Kremer et al. 2020;Antonini & Gieles 2020)).Therefore, the result in the simulation should be multiplied by a factor of 0.36 +1.12 −0.26 , which is reflected in the results of Section 3.4.

A.3. BBHs in the Galactic Nucleus
Stellar-mass BBHs surrounding the supermassive black hole in the galactic nucleus will naturally be in the configuration of a hierarchical triple system, thus undergoing eccentricity oscillation via the EKL mechanism.Based on the observational results, we carried out Monte-Carlo simulations of these BBHs' evolution.The simulations of hierarchical triple systems include the secular equations up to the octupole level of approximation (Naoz et al. 2013a), general relativity precession (Naoz et al. 2013b), and GW emission (Peters 1964).In particular, we take into account three different populations of BBHs in the galactic nucleus: 1. BBHs from the Main Population of Stars in the Galactic Nucleus.This population is expected to have a distance to SMBH within 5 pc, age ∼ 2 − 8 Gyr, total mass ∼ 1.8 × 10 7 M ⊙ (see, e.g., Pfuhl et al. 2011;Launhardt et al. 2002).In particular, we randomly generate BBHs with log uniform distribution in mass and semi-major axis, ranging from 6 to 100 M ⊙ and 0.1 to 50 au, respectively.For the spatial distribution of these BBHs, we adopt the isotropic distribution, with the radial density profile in the galactic center following Hoang et al. (2018).The number of systems is calibrated using the m − σ relation.For conservation purposes, we rule out the systems that are too close to the SMBH to be classified as hierarchical triple systems (Naoz 2016).
After generating the initial condition for BBHs, we evolve these systems numerically, counting their bursting time and number of emitted bursts.The integration is stopped once the system becomes a GW merger or reaches the evaporation timescale (see Eq.( 3) in Hoang et al. 2018).Because of the old age of the main population, we use the steady-state approximation to work out the expectation of bursting sources' number, i.e.: in which Γ rep,all is the replenishing rate of BBHs in the GN, f RB is the fraction of BBHs in the GN that become a repeated burst source, and τ RB is the average lifetime time of bursting sources when they are detectable.Under the steady-state approximation, the total replenishing rate, Γ rep,all , equals the total empty rate at which BBHs merge or evaporate.Therefore, in the simulation, we determine the Γ rep,all by multiplying the inverse of all the BBHs' average lifetime with their total number (∼ 1500) in the inner 0.5 pc.
From the result of the simulation, Γ rep,all ∼ 3 × 10 −6 yr −1 .We note that this quantity equals the combined rate of GW merger plus evaporation, thus is higher than the GW merger rate constrained by LIGO.Moreover, the fraction of BBHs that turns into a bursting source is f RB ∼ 3%, and the average bursting time is τ RB ∼ 1.4 × 10 7 yr.Therefore, we estimate the number of main population bursting sources as ∼ 1 in the inner 0.5 pc of the Milky Way.A similar approach can be applied to calculate the number of bursts, and the result from the simulation is ∼ 100 yr −1 .
2. BBHs in the Young Nuclear Star Cluster According to the observation (Paumard et al. 2006;Lu et al. 2009Lu et al. , 2013)), the YNC has a distance to SMBH within ∼ 0.5 pc, age ∼ 2 − 8 Myr, and total mass ∼ 1.4 − 3.7 × 10 4 M ⊙ , with top-heavy mass distribution.Because of their small mass and young age, we need to go beyond the steadystate approximation and simulate the time evolution of bursting GW sources as a function of time.
For simplicity, we assume that this YNC is born in a starburst at t = 0, and all the massive stars are in the binary system.After evolving for a few Myr, the massive stars have all died and become black holes, and we get the corresponding blackhole mass following the results in Woosley et al. (2002); Belczynski (2020).According to the simulation result, there will be ∼ 100 − 400 BBHs newly formed in the YNC.Because of the existence of mass gap (see, e.g., Bond et al. 1984;Fryer et al. 2001), most of these black holes have mass ∼ 10 M ⊙ .We distribute the BBHs isotropically around the SMBH, with the same density profile as the main population, and trace their evolution as hierarchical triple systems.
As is discussed in Section 3.5, although the YNC has a relatively small total mass, its bursting BBHs number can be as much as 4, which is much higher than the main population with old age and a slow replenishing rate.We show the simulation results in Figure .7. For completeness, we also show the expected number of bursts detected during the LISA mission as a function of the YNC's age (see Figure 9).
3. BBHs in the Inner 1000 au of the Galactic Nucleus The observation of stellar dynamics in the galactic center suggests that there can be up to ∼ 3000 M ⊙ unknown mass within the inner 1000 au of the SMBH (see, e.g., Ghez et al. 2008;GRAVITY Collaboration et al. 2020;Will et al. 2023).Therefore, we explored the possibility that these unknown objects are made up of stellar-mass BBHs (i.e., ∼ 150 BBHs).Since we are only aiming at a heuristic estimation, the simulation assumes that all these BBHs have the mass of 10 − 10 M ⊙ , with the spatial density profile the same as for bursting BBHs in the YNC, and was born at t = 0.
Adopting the same approach as modeling the time evolution of bursting sources in the YNC, we get the properties of this inner 1000 au population (see the hollow grey triangle in Figure 6).In particular, the existence of BBHs within the inner 1000 au, if any, will give a significantly higher fraction of bursting sources.When the cluster is a few Myr old, up to ∼ 13% of all the BBHs can emit GW bursts, contributing to ∼ 20 bursting sources and 200 bursts per year.On the other hand, if there is no such bursting source in future observation, we can put stringent constraints on the existence of BBHs in the center of our galaxy.

B. NUMERICAL APPROACH TO CALCULATE THE SNR OF HIGHLY ECCENTRIC BINARIES
The time domain waveform of eccentric binaries, h(a, e, t) , can be decomposed into different harmonics, with the frequency f n = nf orb (Peters & Mathews 1963;Kocsis & Levin 2012) For highly eccentric binaries, the transient nature makes its frequency domain waveform split into millions of harmonics (peaks, see Figure 3).They are separated by an interval of f orb , and each peak has the width of ∆f ∼ n ḟorb T obs (caused by the orbital frequency shift of GW sources).Therefore, we can transform the expression of SNR into the summation of harmonics: SNR 2 = 8h 2 0 (a) n g(n, e) S n (nf orb )n 2 T obs . (B10) As is shown in Equation (B10), the huge number (millions) of harmonics adds up to the numerical difficulty of calculating the SNR.However, we can further simplify this expression using the knowledge of the envelope of the frequency spectrum.In other words, instead of directly summing over millions of peaks, we can calculate the averaged amplitude in a wider frequency bin, and work out the density of peaks below that envelope (area enclosed by the power spectrum), then get a useful approximation of the SNR.
In particular, let's consider the integration in the frequency domain, take a value of f and the integration interval df around it.For each integration interval, the number of peaks contained in this frequency bin is ∆n ∼ df /f orb , while the amplitude of these peaks can be represented by the average value g(n avg , e), with n avg = f /f orb .
Thus, the summation in Equation (B10) can be turned into the integration of a smoothed function: T obs df f orb = 8h 2 0 (a)f orb T obs g(n avg , e) S n (f )f 2 df (B11) Equation (B11) serves as a fast numerical estimation of highly eccentric binaries' SNR.The accuracy of integration can be changed by adjusting df .In the limit of df = f orb , we recover the strict expression of Equation (B10).

Figure 3 .
Figure3.Comparison between the GW signal of a bursting BBH and a near-circular DWD, both with the same signal-to-noise ratio SNR ∼ 12.We show a BBH system (blue line) with m1 = m2 = 20 M⊙, a = 4.43 au, e = 0.999 and a DWD system (green line) with m1 = m2 = 0.5 M⊙, a = 0.001 au, e = 0.01, both of them are placed at D l = 8 kpc, and observed for T obs = 4 yr.These parameters are chosen such that these sources will have SNR = 12.Upper Panel shows the characteristic strain (Eq.B8) compared with the LISA noise curve ( f Sn(f )) (red line), Middle Panel and Bottom Panel show the timedomain waveforms.This plot highlights that bursting signals in the milli-hertz band can have significant amplitude in the time domain (thousands of times greater than a circular source with the same SNR), even if their averaged SNR is limited.
), as depicted in Appendix B. Particularly, we consider equal mass 10 M ⊙ BBH systems (left panel), and equal mass 0.75 M ⊙ DWD systems (right panel) at a luminosity distance D l = 8 kpc from the detector.The Figure shows the GW signal-to-noise ratio, as a color map of the binary's semi-major axis a and e.The detectable regions are labeled on the map.

Figure 4 .
Figure 4.The Dependence of the (RB) compact binary's detectability on its semi-major axis and eccentricity.Here we show a BBH system (Left Panel) with m1 = m2 = 10 M⊙, and a DWD system (Right Panel) with m1 = m2 = 0.75 M⊙, placed at D l = 8 kpc and observed for 4 yr with LISA.The color represents the compact binary's signal-to-noise ratio as a function of its semi-major axis a and eccentricity 1 − e.As a proof of concept, we depict an example system's evolution (purple line with star and arrows) for both parameter spaces and plot the corresponding orbital time on the top x-axis.The solid lines with color represent the merger timescale τmerger for a given (a, e) configuration of the binary.In other words, when the evolution track of the binary gets into the colored region and gets across a given merger timescale, we can estimate the remaining time for it to emit detectable GW bursts before the merger.We note that the signal-to-noise ratio is suppressed in the region below the dashed lines (tmerger < T obs ) because the binary will merge during the observation, thus having less observational time than the total LISA mission time.Therefore, when calculating the SNR for these binaries following Equation.(B11), we replaced the T obs with tmerger.

Figure 5 .
Figure 5.The maximum detectable distance of the (RB) compact binary.Here we consider the same systems as in Figure 4, but choose different values of luminosity distance D l , and show the equal signal-to-noise ratio (SNR = 5) contours for a 4 yr LISA observation.The left (right) panel shows the BBHs (DWDs) system.The detectability of bursting sources increases at the top-left part of the plot.

Figure 6 .
Figure6.The number of bursting sources and bursts detected per year, assuming a 10-yr LISA mission and detection threshold SNR = 5.The figure depicts the estimation results of bursting BBHs for different channels, with the error bars reflecting the theoretical uncertainties.Specifically, we show the number of detectable bursting sources from each channel (x-axis) and their contribution to the number of bursts detected per year (y-axis).For details of the simulation, see Appendix A. Note that the number of detectable bursting sources is roughly proportional to the observational time.The solid points represent the population with convincing evidence of existence, while the hollow points are for the potentially existing population of bursting sources.We emphasize that if future detection does not show evidence of the population marked in hollow points, it serves as a strong constraint on the corresponding formation channels.

Figure 7 .
Figure 7. Number of detectable bursting sources in the galactic nucleus YNC, as a function of their age.Here we show the expectation of observable bursting systems in the Milky Way center young nuclear cluster, as a function of the cluster's age.The deep blue line stands for a LISA observational period T obs = 10 yr while the light blue line represents T obs = 4 yr.The YNC's age is constrained by observation to be 2 ∼ 8 Myr.

Figure 8 .
Figure 8.An example of the repeated burst phase and inspiral phase during the evolution of an EKL merger.We show a BBH system with m1 = 57.7 M⊙, m2 = 51.9M⊙, a1 = 1.97 au, e1 = 0.38, placed near a SMBH with M = 4×10 6 M⊙, a2 = 4706 au, e2 = 0.97, i = 88 • , D l = 8 kpc .Upper Panel shows the signal-to-noise ratio of this source, for the observational time 4 yr.Middle Panel and Bottom Panel show its eccentricity and semi-major axis evolution as functions of time.As is shown in the figure, this example system undergoes strong EKL oscillation because of the SMBH.It spends a significant fraction of time in the repeated burst phase, bursting detectable GW signal via highly eccentric orbit, before finally reaching the inspiral phase with moderate eccentricity.

Figure 9 .
Figure 9. Number of detectable bursts during LISA mission, as a function of the YNC's age.Here we show the simulation result of observable bursts in the Milky Way center young nuclear cluster, as a function of the cluster's age.The deep blue line assumes the LISA observational period τ obs = 10 yr while the light blue line represents τ obs = 4 yr.The large fluctuation of the light blue line is caused by the limited number of systems in the simulation. ): (1 − e 2 ) [1 + (121/304)e 2 ] .
Robson et al. 2019;nal-to-noise ratio, in the frequency domain, is defined as (see, e.g.,Robson et al. 2019; Chen et al.in which h c is the characteristic strain:h 2 c (a, e, f ) = 4f 2 | h(a, e, f )| 2 , (B8)and h(a, e, f ) is the Fourier transform of the GW signal:| h(a, e, f )| 2 = i is the i-th Bessel function evaluated at ne.