Introduction

The field of gravitational-wave astronomy has been rapidly expanding ever since the detection of the first binary black hole merger GW150914 (B. P. Abbott et al. 2016a). To date, nearly 90 gravitational-wave (GW) sources are cataloged by LIGO-Virgo-KAGRA (LVK) Collaboration, including dozens of binary black holes, two binary neutron stars, and three neutron star-black hole mergers (B. P. Abbott 2023a). Additional GW sources were independently cataloged (A. H. Nitz et al. 2021a;S. Olsen et al. 2022;A. H. Nitz et al. 2023;D. Wadekar et al. 2023;A. K. Mehta et al. 2023) using publicly available data (B. P. Abbott et al. 2021a;B. P. Abbott 2023b). The third observation (O3) run of Advanced LIGO (J. Aasi et al. 2015) and Advanced Virgo (F. Acernese et al. 2015) led to the detection of compact objects under 3 M e with GW190814 (R. Abbott et al. 2020), where the secondary compact object had a mass of ∼2.59 M e and a low spin (0.07). Several other events like GW190425, GW191219, GW200105, GW200115, and GW200210 identified during this run also had one of the component masses less than 3 M e . The ongoing fourth observation run detected GW230529 (A. G. Abac et al. 2024). This event's primary object had a mass ranging between 2.5 and 4.5 M e , making it an additional compact object, likely a black hole, existing within the "lower mass gap" (C. D. Bailyn et al. 1998;F. Özel et al. 2010; W. M. Farr et al. 2011). Although several studies have provided insights into the mass and spin distributions of compact sources detected through current GW detectors (J. Roulet et al. 2019;B. P. Abbott et al. 2021b), the possibility of discovering ultracompact objects with masses less than a solar-mass range remains an open question (B. P. Abbott et al. 2019;A. H. Nitz et al. 2021b; K. S. Phukon et al. 2021;R. Abbott et al. 2022;A. H. Nitz et al. 2022;LVK Collaboration 2023;A. L. Miller 2024).

Subsolar-mass (SSM) compact objects do not follow the standard stellar evolution pathway. These objects, if black holes, are expected to form through nonstellar evolution models and could be primordial black holes (PBHs; B. J. Carr et al. 2010). If they are neutron stars (V. Doroshenko et al. 2022), they might result from nonstandard supernova explosion models (B. Müller et al. 2024). Although the search for SSM black holes began quite early (T. Nakamura et al. 1997;C. Alcock et al. 2000;B. Abbott et al. 2005B. Abbott et al. , 2008)), no candidates have been found yet. Since then, numerous models have proposed various formation pathways for these sources. The most common formation mechanism of PBHs suggests their origin from the direct collapse of early, small-scale fluctuations (Y. B. Zel'dovich & I. D. Novikov 1967;S. Hawking 1971) due to certain features of the inflationary potential. Additionally, there are alternative formation channels where PBHs emerge from phase transitions (C. T. Byrnes et al. 2018) in the early Universe or through the collapse of topological defects like cosmic strings (S. W. Hawking 1989;A. Polnarev et al. 1991;C. Hong-bo et al. 1996).

If an SSM compact object, whether a black hole or neutron star, appears in a binary system, the emitted GW can be detected through ground-based interferometers. Several studies have investigated the search for SSM black holes (B. P. Abbott et al. 2018Abbott et al. , 2019;;A. H. Nitz et al. 2021bA. H. Nitz et al. , 2021c;;R. Abbott et al. 2022), but no significant detections have been made to date. Studies show that a small fraction of dark matter could be due to PBHs (B. J. Carr et al. 2010;B. Carr et al. 2021). While many cosmological investigations have ruled out their existence at extremely low masses (M. Sasaki et al. 2018), exploration continues in a mass range spanning several orders of magnitude. A confirmed detection within the LIGO-Virgo frequency band would provide critical insights into the formation mechanisms of PBHs and contribute to constraining the fraction of dark matter in the universe. On the other hand, a SSM neutron star in a binary system could provide a high tidal effect on its companion (F. Crescimbeni et al. 2024a). These effects would be measured with the future-generation GW detectors (F. Crescimbeni et al. 2024b).

The offline search for GWs from merging binaries uses the matched-filtering technique (B. S. Sathyaprakash et al. 1991; S. V. Dhurandhar et al. 1994aDhurandhar et al. , 1994b;;B. J. Owen et al. 1999;B. Allen et al. 2012;S. A. Usman et al. 2016;G. S. Davies et al. 2020). In this method, a bank of modeled signals, or templates, is correlated with well-calibrated interferometer data (X. Siemens et al. 2004;J. Abadie et al. 2010). However, this approach becomes computationally demanding, particularly for low-mass binaries, as the cost increases with the number of templates and the length of the signal model used as a matched-filter template. To mitigate the computational challenges, suboptimal choices are often made by limiting the search parameters. For instance, searches may only filter data above 45 Hz (B. P. Abbott et al. 2016b) or limit the duration of the templates to nearly 512 s (A. H. Nitz et al. 2021bNitz et al. , 2022)). While these restrictions help reduce computational costs, they also reduce the sensitive volume by approximately 24%, within which PBHs might be detected.

Observing long-duration GW signals poses significant challenges with existing search methods. The main difficulty arises from the necessity of using a very dense template bank ( ( )  10 7 ), which significantly increases the computational cost of the matched-filtering search. Furthermore, the search sensitivity can be compromised by nonstationary data, which may contain long-duration correlations that hinder current signal-vetoing techniques and statistical analyses. This nonstationarity can also impact the statistics and signal-vetoing methods used in current search pipelines. These issues are expected to become considerably more severe in the era of third-generation (3G) ground-based detectors-3G detectors such as the Cosmic Explorer (B. P. Abbott et al. 2017;D. Reitze et al. 2019) and the Einstein Telescope (S. Hild et al. 2009;M. Punturo et al. 2010;M. Maggiore et al. 2020;S. Di Pace et al. 2022;A. Grado 2023) are anticipated to detect binary mergers at rates 2 to 3 orders of magnitude higher than current detectors (M. Evans et al. 2021). These detectors, designed to operate from very low frequencies (starting from 2 Hz), will observe signals for several minutes or hours. Due to their high sensitivity in the lower-frequency band, the likelihood of detecting eccentric or precessing binaries will be higher, which will indirectly expand the template bankʼs parameter space both in dimensionality and parameter ranges, thereby increasing the computational cost of the search. Furthermore, since signals would remain in the sensitivity band for longer periods, the Earthʼs rotation will reduce search sensitivity by altering the detectorʼs response functions and affecting matched-filter statistics. Therefore, developing an efficient, cost-effective matched-filtering strategy for longduration GW signals is essential to advance the current stateof-the-art search techniques.

One approach to efficiently search for long-duration signals, such as those from SSM binaries, is implementing a hierarchical search strategy (R. Dhurkunde et al. 2022;K. Soni et al. 2022). In this method, a two-stage matched-filtering search is performed using multiple template banks of varying densities. In the first stage, the search is conducted over coarsely sampled data using a coarse bank to identify coincident triggers that could represent true GW events. These triggers are followed up in the second stage with a finer search, focusing on the neighborhood of the parameter space identified in the first stage. This two-stage approach effectively reduces the number of matched-filtering operations required for the search, significantly reducing computational time.

In this paper, we demonstrate how hierarchical search strategies can be effectively applied to search for long-duration GW signals from compact binary mergers in Advanced LIGO data. Using the SSM compact object search as a case study, we show that performing the search at lower frequencies in Advanced LIGO and tuning each stage of the existing hierarchical search strategy (K. Soni et al. 2022Soni et al. , 2024) not only reduce the computational cost of matched-filtering process but also enhance the search sensitivity. Specifically, we show that our approach increases the searchʼs volumetric sensitivity by approximately 10%-20% while reducing computation time by a factor of 2.5 compared to the standard two-detector PyCBC search. This methodology also accounts for the increase in background noise due to the larger template bank, as well as signal consistency tests and potential signal-to-noise ratio (SNR) losses resulting from approximations introduced during the hierarchical stages. Furthermore, we discuss how hierarchical search strategies could be generalized for computationally intensive compact binary coalescence (CBC) searches in 3G detectors.

Method

The matched-filtering search for long-duration signals as in the case of binaries containing SSM compact objects is expensive as it requires a very dense template bank. To optimize the search, often the length of the template is reduced to a manageable duration (∼512 s) so that search analysis can be performed. This could be enabled by performing matched filtering from 45 Hz rather than 15 Hz (B. P. Abbott 2023a). However, such adjustments affect the horizon distance of the binary and the expected SNR.

The horizon distance (K. S. Thorne 1987; B. Allen et al. 2012) for an inspiraling binary is given by

where f min and f max are the minimum and maximum frequencies of the LIGO's sensitivity range. S n ( f ) is the power spectral density of the noise in the detector and ρ is the expected matched-filter SNR for an inspiraling binary with chirp mass  observed in the detector's frame. For a particular source of chirp mass of a few that have a fixed SNR in the detector's frame, changing the operating frequency band can affect the detectability of a signal. This means that the fractional SNR loss, as also shown in R. Magee et al. (2018), due to a change in the operating frequency band would be

, 2

SNR loss max min max max with respect to the standard operational band of Advanced LIGO for matched filtering for generic CBC search (B. P. Abbott 2023a).

If we perform a search from 35 Hz, assuming the other source parameters do not change, the percentage loss in SNR is about 3.1% for a frequency band of 35-2048 Hz. This loss is relatively small compared to the SSM searches conducted in the 45-2048 Hz band (R. Abbott et al. 2022), which experience an SNR loss of about 8%-9%. This comparison can also be seen in Figure 1. Ideally, the lower-frequency limit can be further lowered to match that used for a generic CBC search. However, this step will increase the computation demand and require very long-duration templates in the bank. Therefore, in our work, we select a lower-frequency cutoff of 35 Hz. By making this choice, we expect the loss in astrophysical volumetric sensitivity to the inspiral stage by ∼10% of binary coalescence, which is lower than 24% as observed in R. Abbott et al. (2022).

Generating templates at lower frequencies enhances the sensitivity of the search, but it can also lead to a higher density of the template bank. This increase in density may raise the computational cost of matched filtering in traditional offline PyCBC or flat search (S. A. Usman et al. 2016;G. S. Davies et al. 2020). However, this added computational burden can be mitigated by adopting hierarchical search strategies.

Review of Hierarchical Search

The two-detector flat search performs matched filtering over a discretely sampled data segment using a dense bank of templates and generates triggers with SNRs (ρ) (S. A. Usman et al. 2016). In contrast, a hierarchical search strategy offers a more efficient approach by enabling a multistage matched-filtering process, where the number of templates is progressively reduced at each stage. As described in K. Soni et al. (2022), this search divides a flat search into two stages: coarse and fine. During the coarse stage, the data is matched filtered using a coarse template bank, which consists of sparsely placed templates. The sparseness of these templates is determined by the minimal match (B. J. Owen 1996) at which the bank is constructed, typically set lower (below 0.97) than the value used for constructing a bank for the flat search.

Performing a coarse search reduces the number of matchedfiltering operations. This reduction is further enhanced when the data are sampled at a lower frequency (512 Hz). However, this approach comes with the trade-off of potentially lower matched-filter SNR values for the resulting triggers. To address this, the SNR thresholds are lowered (ρ = 3.5) compared to those used in a flat search (ρ = 4). Given these triggers could be generated by non-Gaussian features or glitches (B. P. Abbott et al. 2021c) present in the data, the SNRs are further downweighted using chi-square vetoes (B. Allen 2005; A. H. Nitz 2018). Only those triggers that pass these vetoes are then subjected to a coincidence test (S. A. Usman et al. 2016), during which a ranking statistic is computed to assess their significance (A. H. Nitz et al. 2017).

The coincident triggers obtained from the coarse search, with ranking statistics above a certain threshold (approximately 7, as used in K. Soni et al. 2022), are followed up in the second stage for a finer search. In this stage, a focused search is conducted within the vicinity or neighborhood (nbhd) of the coarse templateʼs parameter space. This nbhd is a region around a coarse template where the minimal match between templates within the nbhd (∼10-100) and the coarse template ranges from 0.75 to 0.99.

To avoid the computational burden of calculating the nbhd for each coarse template on the fly, a precomputed nbhd bank is used. This bank contains all the nbhd regions and the corresponding templates for each coarse template. During the second-stage search, a union of all the nbhds corresponding to the coarse triggers in each data segment is used for matched filtering. To further improve the SNR, the data sampling rate is increased to 2408 Hz, which is higher than that used in the coarse search. Triggers with SNRs above 4 that pass all chisquare tests are then subjected to a final coincidence search, compiling a list of foreground candidates.

To assess the significance of potential GW events, the false alarm rate (FAR) is estimated based on the background (S. A. Usman et al. 2016). Unlike the flat search, which estimates the background by applying millions of time shifts to triggers from a single detector, the hierarchical search employs a hybrid approach in its second stage (K. Soni et al. 2024). At first, a few time shifts are applied to generate coincident background triggers. Then, an exponential fit is applied to the cumulative distribution of these background triggers, and the fitted curve is used to calculate the FAR for the foreground triggers obtained in the second stage. This method of estimating the background is particularly effective for longduration signals, as the expected background distribution tends to follow the tail of a Poisson distribution (S. A. Usman et al. 2016).

Template Bank

We construct two aligned-spin banks-a coarse bank at a minimal match of 0.92 and a fine bank of 0.97-using a geometric placement algorithm (I. W. Harry et al. 2014) for the hierarchical search. Both banks are designed to cover parameters where m 1 ranges between 0.2 and 10 M e and m 2 between 0.2 and 1.0 M e in the detector frame. The dimensionless spins are restricted to 0.1 for component masses smaller than 0.5 and 0.9 for larger masses. The orbital precession and the existence of multipoles are neglected in the search space. These bank parameter ranges are consistent with the bank used for the LVK search (R. Abbott et al. 2022; LVK Collaboration 2023). From here, we refer to this bank as the flat bank .

In contrast to the flat bank, where templates commence at a frequency of 45 Hz, the templates in our banks start at a frequency of 35 Hz. This choice reduces the fractional loss in the matched-filter SNR, as shown in Section 2. As a result, even though our coarse bank is constructed at a lower minimal match, it is approximately 1.6 times the size of the flat bank. These distinctions are summarized in Table 1.

The hierarchical search requires the construction of the nbhd bank for the second-stage search. Therefore, we use the generated fine bank to construct the nbhd bank using the method described in Section IV of K. Soni et al. (2022).

Figure 2 shows the parameter space covered by the coarse templates in the chirp mass and effective spin plane. This plot shows that the distribution of templates within a nbhd is not uniform across the parameter space, primarily due to boundary effects. These effects occur because the match between neighboring templates gradually decreases as the mismatch in the τ 3 mass parameter increases relative to τ 0 , thereby significantly extending the nbhd region along this coordinate.

As shown in Figure 2, the number of templates in an nbhd typically ranges from a few tens to hundreds. This can significantly affect the final background in the second stage. To improve background estimation, more noise coincidences from the coarse search need to be followed up. However, since the number of templates in the nbhds is relatively small, the search cost is not expected to be significantly higher than that of the flat search.

Application to SSM Search

We perform an SSM search on publicly available data sets using a two-stage hierarchical approach, as outlined in Section 2. The data consist of approximately 5 days of coincident observations from the O3 run of the two LIGO detectors LIGO Hanford and LIGO Livingston, covering the period from 2019 April 1 to 8.

To begin, we first conduct a coarse search over the data sampled at 512 Hz using the coarse bank described in Section 2.2. The templates are generated at 35 Hz using the TaylorF2 (R. Sturani 2010) waveform model with phase corrections up to 3.5 post-Newtonian order. The lengths of these templates range between 10 2 and 10 3 at 35 Hz. To prevent the templates from wrapping around during the fast Fourier transform operation in the matched-filtering step, we ensure that the data segment length exceeds that of the longest template in the bank. Consequently, we set the data segment length to 2048 seconds for our analysis.

We identify triggers with a matched-filter SNR and reweighted SNR (B. Allen 2005; S. A. Usman et al. 2016) of 3.5 or greater. This threshold is chosen to increase the likelihood of detecting true signals that might be missed due to lower data sampling and the use of a coarse bank. To reduce the impact of short-duration glitches in the data, the triggers are further weighted using a chisquare and sine-Gaussian vetoes. The surviving triggers then undergo a coincidence test, where they are shifted in time by 5000 seconds, and a ranking statistic (Λ) (A. H. Nitz et al. 2017; G. S. Davies et al. 2020) is computed.

In the next step, we perform a finer search in the second stage on coincident triggers with Λ 7. During this stage, matched filtering is conducted again on data sampled at 2048 Hz, using a union of nbhds around the identified coarse templates. The data sampling rate is increased by a factor of 4 to improve the matched-filter SNR. Triggers from this stage are collected if their SNRs and reweighted SNRs exceed a threshold of 4. These triggers are then reweighted using chisquare and sine-Gaussian vetoes before undergoing a coincidence test, over the same time-shift interval as in the first stage. Finally, a list of foreground and background triggers is compiled, and the FAR for the foreground triggers is calculated following the procedure described in K. Soni et al. (2024).

The primary differences between the SSM search using the hierarchical method and the search adopted by the LVK (R. 45 Hz, our approach begins at a lower frequency of 35 Hz. This choice enhances the sensitivity of our search to potential GW signals, particularly those that might be present at lower frequencies. By adjusting the matched-filtering start frequency, we aim to improve the detection capabilities for signals that might be overlooked in the flat search. Note that our analysis utilizes data from the two LIGO detectors, whereas the LVK analysis incorporates data from the Virgo detector as well. However, this distinction does not impact the comparison of results in Section 3.1, as both analyses consider two-detector configurations for consistency.

Results

The hierarchical search yielded a list of GW candidates, many of which were statistically insignificant due to their FAR values exceeding 1 per year. These candidates along with the ones identified through a flat search using the same data set are summarized in Table 2. While a few candidates are common to both search pipelines, none are statistically significant. As shown in Figure 3, the foreground events overlap with the background distributions for both searches, indicating that the candidates are primarily noise coincidences.

Figure 3 also highlights the disparity in backgrounds between hierarchical and flat searches, with the hierarchical search producing a larger background. This difference is due to the usage of different numbers of templates within their respective search pipelines. As shown in Table 1, the coarse bank is approximately 1.6 times and the fine bank is 4.8 times denser the flat bank owing to the template generation at 35 Hz. If the flat search had employed the fine bank, its background distribution would likely resemble that of the hierarchical search, especially in the tail. Although the number of templates used in the second-stage search is low (∼10-1000 per data segment), the background generated in the second stage is expected to increase leading to more instances of noise coincidence. However, this also improves the chances of detecting GW sources that might be missed in a flat search.

Sensitivity and Search Efficiency

The sensitivity of a search is determined by how many signals it can detect at a particular significance level within a given observation time (T). This can be quantified by estimating the observable volume-time (VT) product (V. Tiwari 2018). For a constant merger rate of the population of binaries, the average VT sensitivity product is given by

where N det is the number of detected sources in the search and N inj is the total number of injected sources. V 0 is the volume defined as

where dV dz c is the differential comoving volume in an expanding universe with redshift z.

To test the sensitivity of our search method, we conducted a comparative analysis through an injection campaign on a simulated binary population. In this population, we assumed that one of the compact objects has a mass below a solar mass, while the other ranges from 1 to 10 solar masses. We created three distinct sets of injections, each defined by different spin conditions: high spin, low spin, and a mixed case where only one compact object has low spin. The distributions and ranges of the component masses and spins for these three scenarios are detailed in Table 3.

For the given parameter space, we generated GW signals using the waveform approximants listed in Table 3, starting at a frequency of 35 Hz. Each signal was injected into the data with a minimum interval of 100 seconds between injections, assuming an isotropic distribution for the sky locations of the sources. Following this procedure, approximately 6500 injections were made across the three sets, and a search was performed using both the traditional flat method and our hierarchical approach.

Figure 4 shows the VT ratio computed for the two search methods. The hierarchical search outperforms the flat search, for each identified candidate. The FARs of the detected events in the flat search were determined using the time-shift method, whereas those for the hierarchical search were determined by the method described in K. Soni et al. (2024). showing an improvement in the VT ratio by approximately 1.1-1.2. This improvement is primarily due to the SNR enhancement, when the search is performed at a lower frequency of 35 Hz instead of 45 Hz. When the search is performed at 35 Hz, we expect the gain in the SNR to be approximately 6% and an astrophysical volumetric gain of ∼20%. However, the injection study results indicate that the detection volume improves by only 10%-20%. This discrepancy likely arises from an increase in the noise background when the lower bound of the matched-filtering frequency is reduced.

We compared the matched-filtering cost by comparing the number of CPU core hours required by each detector in use in the two searches. As shown in Table 4, the number of CPU core hours required by flat search is more than coarse and fine searches owing to the different number of templates used in each search. The numbers show that the overall cost of a hierarchical search is approximately 2.5 times less than a flat search. This is a huge advantage of hierarchical search given that the search sensitivity is also more than that of the flat search.

Conclusion and Discussion

The search for long-duration GW signals from compact binary mergers, such as SSM binaries, low-mass BNSs, precessing binaries, and binaries with moderate-mass ratios (m 1 /m 2 , m 1 m 2 ) and noncircular orbits, is challenging with the current matched-filtering search methods. This is primarily because of the requirement for large, densely populated template banks. To mitigate the computational burden, suboptimal choices are often made, which inevitably limit the sensitivity of such searches. With the advent of 3G detectors, these challenges are expected to become more pronounced as GW signals will be observed over longer durations, ranging from tens to hundreds of minutes. This extended observation window will substantially increase the computational demands due to the rapid expansion of the parameter space. Consequently, the development of efficient hierarchical search strategies is critical, not only to enhance current detection capabilities but also to ensure readiness for the vastly more computationally expensive searches required by next-generation detectors.

In this paper, we demonstrated for the first time how the hierarchical search strategy can be used to search long-duration signals, such as those from binaries containing an SSM compact object, without restricting the parameter space. In Section 2, we presented a preliminary calculation indicating that the SNR improves by approximately 6% when the matched filtering is conducted starting at a frequency of 35 Hz. This means that an approximately 20% increase in the sensitive

Table 3 Overview of Three Injection Sets Focusing on One of the Component Masses in SSM Ranges 3). The VT ratio, binned over inverse false alarm rates (IFARs), demonstrates an improvement in search sensitivity across all chirp mass (M chirp ) bins. Our findings in this paper highlight that near-optimal sensitivity can be achieved using a hierarchical search strategy for very long-duration signals, even compared to a direct flat search. We specifically showed that by optimizing different aspects of the search process in two stages, such as adjusting the frequency of operation, data sampling rates, and the density of the template banks, we achieved computational savings of up to a factor of 2.5 while simultaneously enhancing the sensitivity of the SSM search.

Looking ahead, 3G detectors are expected to introduce several significant challenges for CBC searches due to their enhanced low-frequency sensitivity, as discussed in Section 1. With the ability to detect GW signals from a broader range of CBC sources, including eccentric and precessing binaries over extended durations, the search space will expand dramatically. This expansion in search space will, in turn, increase computational demands to potentially unmanageable levels. Therefore, a hierarchical approach may be proposed to address these challenges effectively.

In the 3G era, the hierarchical search could be structured into stages, with the first stage focused on efficiently identifying potential GW candidates. This efficient identification could be achieved by focusing on certain features in the signalʼs morphology and maximizing the likelihood of detecting it across those source parameters. In contrast, the second stage aims to improve the SNR and address any losses from the first stage. Since the primary goal of the first stage is to locate regions where signals are likely to be present, general optimizations such as reducing the template bank size, adjusting data sampling, and defining the operating frequency range for matched filtering can be applied, as demonstrated in this work. The density of the template bank could be reduced by coarsening the bank and adjusting the frequency at which templates are generated. However, this step should focus on regions of the parameter space where a higher SNR is expected due to features in the binaryʼs orbit. For example, in eccentric binaries, the effects of eccentricity are most prominent at lower frequencies, while at higher frequencies, the binaryʼs orbit is expected to circularize. Therefore, the template bank, matchedfiltering frequency band, and data sampling rate could be adjusted to prioritize higher frequencies, reducing the search cost. To enhance detection probability, an nbhd search could be performed in the second stage with templates incorporating eccentricities at lower frequencies. A similar strategy could be applied to binaries with moderate precession. For these binaries, the first stage can be adjusted to search for signals with variable starting frequencies for matched filtering, excluding less significant merger-ringdown phases, which could be relaxed in the second stage. Since the Earthʼs rotation is expected to impact the antenna response functions, potentially reducing search sensitivity, approximate response functions that include the effect of the sourceʼs orientation with respect to the detector could be introduced only in the second stage of the hierarchical search, thereby improving overall sensitivity. The second stage could also be optimized to recover any SNR losses from the first stage due to the various optimizations applied earlier.