# Hierarchical search method for gravitational waves from stellar-mass binary
black holes in noisy space-based detector data

Yao Fu  MOE Key Laboratory of Fundamental Physical Quantities Measurements,
PGMF, Department of Astronomy and School of Physics, Huazhong University of
Science and Technology, Wuhan 430074, China  Yan Wang  0000-0001-8990-5700
<EMAIL_ADDRESS>MOE Key Laboratory of Fundamental Physical Quantities
Measurements, PGMF, Department of Astronomy and School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China  Soumya D. Mohanty
0000-0002-4651-6438<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of Texas Rio Grande Valley, Brownsville, Texas 78520,
USA  Department of Physics, IIT Hyderabad, Kandai, Telangana-502284, India

###### Abstract

Future space-based laser interferometric detectors, such as LISA, will be able
to detect gravitational waves (GWs) generated during the inspiral phase of
stellar-mass binary black holes (SmBBHs). These detections contain a wealth of
important information concerning astrophysical formation channels and
fundamental physics constraints. However, the detection and characterization
of GWs from SmBBHs poses a formidable data analysis challenge, arising from
the large number of wave cycles that make the search extremely sensitive to
mismatches in signal and template parameters in a likelihood-based approach.
This makes the search for the maximum of the likelihood function over the
signal parameter space an extremely difficult task, with grid-based
deterministic global optimization methods becoming computationally infeasible.
We present a data analysis method that addresses this problem using both
algorithmic innovations and hardware acceleration driven by Graphics
Processing Units (GPUs). The method follows a hierarchical approach in which a
semi-coherent $\mathcal{F}$-statistic is computed with different numbers of
frequency domain partitions at different stages, with multiple particle swarm
optimization (PSO) runs used in each stage for global optimization. An
important step in the method is the judicious partitioning of the parameter
space at each stage to improve the convergence probability of PSO and avoid
premature convergence to noise-induced secondary maxima in the semi-coherent
$\mathcal{F}$-statistic. The hierarchy of stages confines the semi-coherent
searches to progressively smaller parameter ranges, with the final stage
performing a search for the global maximum of the fully-coherent
$\mathcal{F}$-statistic. We test our method on 2.5 years of a single LISA time
delay interferometry (TDI) combination and find that for an injected SmBBH
signal with a signal-to-noise ratio (SNR) between $\approx 11$ and $\approx
14$, the method can estimate (i) the chirp mass with a relative error of
$\lesssim 0.01\%$, (ii) the time of coalescence within $\approx 100$ sec,
(iii) the sky location within $\approx 0.2$ ${\rm deg}^{2}$, and (iv) orbital
eccentricity at a fiducial signal frequency of $10$ mHz with a relative error
of $\lesssim 1\%$.

##  I Introduction

Since the first direct detection of a gravitational wave (GW) signal from a
stellar-mass binary black hole (SmBBH) merger, GW150914 [1], the twin LIGO
observatories [2] in coordinated observations with the Virgo [3] and KAGRA [4]
detectors have recorded nearly $100$ events [5, 6, 7, 8] during the first
three observing runs (O1 to O3). Due to a lower frequency bound on the
sensitive bandwidth of ground-based detectors arising from seismic noise, all
observed signals have frequencies higher than $\approx 20$ Hz, which
corresponds to the late inspiral, and merger phases of stellar-mass binary
systems lasting for a $\sim 1$ min or less. Even with third-generation
instruments in the future, such as the Einstein Telescope [9] and Cosmic
Explorer [10], ground-based detection of GW signals will ultimately be limited
to $\gtrsim 1$ H due to gravity gradient noise [11]. Observations of GW
signals in the sub-Hz to millihertz range must necessarily use space-based
detectors, projected to be operational in the next decade starting with the
launch of the Laser Interferometer Space Antenna (LISA) [12]. At around the
same time, the proposed TianQin [13, 14] and Taiji [15, 16] (or one of them)
may also join LISA to provide a network of space-based GW detectors that will
have a significantly higher combined capability than any of the individual
ones [17, 18, 19, 20].

LISA will be sensitive to GW signals from a variety of sources in the
millihertz range [21, 22], such as galactic compact binaries, extreme-mass-
ratio inspirals (EMRIs), supermassive binary black hole mergers, SmBBHs,
intermediate-mass binary black holes, and possibly the stochastic background
from phase transitions and cosmic strings in the early Universe. Of particular
interest to this paper is the fact that, complementary to ground-based
detectors, LISA will be able to detect GWs in the early inspiral phase of
SmBBHs in the chirp mass range $10-100~{}M_{\odot}$. It is expected that we
will observe multiple SmBBHs [23, 24, 25, 26] throughout the mission lifetime
of LISA and when combined with other proposed space-based GW detectors, such
as TianQin [13, 14], Taiji [15, 16], and DECIGO [27, 28], the number of
detectable SmBBHs will be further increased [29].

Signals from SmBBHs will contain rich information, which will allow us to
study their astrophysical formation channels [30, 31, 32, 33] and constrain
fundamental physics theories [34, 35, 36, 37]. With a reasonable accuracy in
the predicted time of coalescence of an SmBBH system that has a lifetime to
merger of a few years, space-based detectors will be able to provide alerts to
both ground-based GW detectors as well as electromagnetic observatories,
enabling multi-wavelength and multi-messenger observations [23, 38, 39, 40] of
such events.

However, in order to realize the rich scientific prospects outlined above we
have to first overcome the major data analysis challenges [26] posed by the
detection and parameter estimation of SmBBHs signals. Unlike ground-based
detectors, SmBBH signals will persist in the sensitive frequency range of
space-based detectors for several years, spanning more than $10^{5}$ wave
cycles. This makes a matched filter-based search extremely sensitive to
mismatches in the parameters of the signal and the templates, requiring an
estimated ($10^{30}-10^{40}$) [26] templates to cover the entire signal
parameter space. In addition, each template, for 2.5 years of observations
with a sampling interval of 10 seconds, will require $\approx 10^{7}$ samples,
making the evaluation of each matched filter output computationally expensive.
In fact, the computational cost associated with SmBBH signals is comparable to
that of the more widely known problem of matching filter-based search for EMRI
signals [41, 42, 43] and arises from essentially the same causes, namely, low
amplitude signals lasting over many cycles requiring long integration times to
accumulate sufficient signal-to-noise ratio (SNR).

Several studies have been conducted to assess the SmBBH data analysis
challenge or propose exploratory methods to address it. For example, in [29],
the Fisher information matrix was used to evaluate the parameter estimation
accuracy of SmBBHs for LISA. As shown in [44], archival searches can uncover
many low SNR GW signals in LISA data. In general, methods that have been
proposed so far for SmBBH searches have strong underlying assumptions: In [44,
45], the signal parameter search space is significantly restricted, while the
performance of the methods in [46, 47, 48] have only been studied for
noiseless data.

In this paper, we present a data analysis method based on a hierarchical
approach enabled by Particle Swarm Optimization (PSO) [49, 50, 51] that
overcomes some of the above limitations and can search for SmBBH signals over
a wide parameter range in the presence of noise, thus significantly advancing
the state-of-the-art in this problem. PSO is a stochastic global optimization
method that has been applied successfully to a wide range of GW data analysis
problems [52, 53, 54, 55, 56], and proves to be highly effective in addressing
the present problem too. The hierarchical search is comprised of several
stages of semi-coherent searches, inspired by the ones used extensively in
continuous GW searches [57, 58, 59] with ground-based detectors, that are
terminated by a fully coherent one. An important feature of our method is the
use of parameter space partitioning at each stage to avoid premature
convergence to secondary noise-induced maxima of the semi-coherent search
statistic and improve the convergence of PSO to its global maximum. In
addition, we use multiple PSO runs at each stage and in each partition to
further improve its performance.

A hierarchical search has been investigated in [60] for the case of noiseless
data. Our method differs significantly in several aspects, such as the use of
parameter space partitioning, maximization over all extrinsic parameters in
the semi-coherent log-likelihood functions, and multiple PSO runs. We note
that the problem of secondary maxima induced by noise is one of the major
challenges in this search that does not appear in noiseless data.

The rest of the paper is organized as follows. Sec. II describes in detail the
search statistic. Sec. III provides an overview of the hierarchical search
method and its implementation in this paper. Sec. IV describes how we
accelerate our code with the frequency-domain response of the detector and
GPU-based parallelization. In Sec. V, we present our main results on the
performance of the method using simulated data realizations. Our conclusions
and discussion of future improvements are presented in Sec. VI. Several
technical details have been relegated to two appendices. Throughout this
paper, we work in natural units in which $G=c=1$.

##  II $\mathcal{F}$-statistic and Semi-Coherent Search

In this section, we provide a brief overview of the Generalized Likelihood
Ratio Test (GLRT) formalism [61], commonly called the fully-coherent
$\mathcal{F}$-statistic [62] in GW data analysis, along with our hierarchical
search scheme based on the semi- and fully-coherent $\mathcal{F}$-statistic.

###  II.1 Fully-coherent $\mathcal{F}$-statistic

The data from a space-based detector such as LISA consists of a set of time
series $s=\\{s^{I}(t)\\}$, $t\in[0,T]$, where $I$ denotes a specific time-
delay interferometry (TDI) combination [63, 64] used to suppress laser
frequency noise. Commonly used TDI combinations are the Michelson combinations
$I\in\\{X,Y,Z\\}$ and combinations called $\\{A,E,T\\}$ that have mutually
independent noise of instrumental origin [65].

The detection of GW signals in noise requires deciding between the null
hypothesis $H_{0}$, namely, the data $s$ does not contain a GW signal, and a
set of alternative hypotheses, denoted as $H_{1}$, corresponding to the
presence of an additive GW signal characterized by parameters $\Theta$. The
most commonly used decision rule in searches for parameterized GW signals with
theoretically known waveforms is the Generalized Likelihood Ratio Test (GLRT)
defined as [66]:

$\displaystyle\Lambda_{G}$ $\displaystyle=$
$\displaystyle\max_{\Theta}\Lambda(\Theta|s)\;,$ (1)
$\displaystyle\Lambda(\Theta|s)$ $\displaystyle=$
$\displaystyle\ln\left(\frac{p(s|H_{1};\Theta)}{p(s|H_{0})}\right)\,,$ (2)

where $\Lambda(\Theta|s)$ is called the log-likelihood ratio, and the joint
probability density functions describing the data under the hypotheses $H_{0}$
and $H_{1}$ (with signal parameters $\Theta$) are denoted by $p(s|H_{0})$ and
$p(s|H_{1};\Theta)$, respectively.

Assuming that the instrumental noise in each TDI combination is a realization
of a stationary Gaussian stochastic process and that the TDI combinations
chosen have mutually independent noise,

$\Lambda(\Theta|s)=\sum_{I}\left[\langle
s^{I},h^{I}(t;\Theta)\rangle_{I}-\frac{1}{2}\|h^{I}(t;\Theta)\|^{2}_{I}\right]\,.$
(3)

Here, $h^{I}(\Theta)$ represents the GW signal, $\langle a,b\rangle_{I}$ is
the noise-weighted inner product defined as [67],

$\langle
a,b\rangle_{I}=4\,\mathrm{Re}\int_{0}^{+\infty}\frac{\widetilde{a}(f)\cdot\widetilde{b}^{*}(f)}{S_{n}^{I}(f)}df\,,$
(4)

where $\widetilde{a}(f)$ and $\widetilde{b}(f)$ are the Fourier transforms of
continuous time functions $a(t)$ and $b(t)$, respectively, and $S_{n}^{I}(f)$
is the one-sided power spectral density (PSD) of the noise. The norm induced
by the inner product is denoted as $\|a\|^{2}_{I}=\langle a,a\rangle_{I}$. For
simplicity, in this work, we only employ the TDI X-combination for which the
corresponding one-sided PSD is provided in [63, 12]. While adding data from
other combinations is not fundamentally challenging, it increases the
computational costs proportionally and exceeds the computational resources
available for this study. As such, we drop the use of the TDI combination
index $I$ from here on.

In this paper, we employ the restricted post-Newtonian waveform family
described in Appendix A that is characterized by a $10$-dimensional parameter
space. The set of parameters is
$\Theta=\\{t_{c},M_{c},\delta_{\mu},\theta,\phi,e_{0},D_{L},\psi,\iota,\phi_{c}\\}$,
which denote, starting from $t_{c}$, the time of coalescence, the chirp mass,
the dimensionless mass difference, the sky location angles
($\theta=\pi/2-\beta$, $\phi=\lambda$, where $\lambda$ and $\beta$ are the
ecliptic longitude and ecliptic latitude, respectively), orbital eccentricity
at a fiducial signal frequency of $10$ mHz, luminosity distance, polarization
angle, the angle between the line of sight and the angular momentum direction,
and the phase at the coalescence. We consider only the $(2,2)$ mode of the GW
signal and leave the inclusion of higher-order modes arising from significant
orbital eccentricity, spin, tidal effects, and other influences, to future
work.

The GLRT involves maximizing the log-likelihood ratio defined in Eq. 3 over
the space of signal parameters. The value of the global maximum serves as the
detection statistic that is compared to a threshold, and if a detection occurs
due to the crossing of the threshold, its location in parameter space provides
the Maximum Likelihood Estimate (MLE) of the parameters of the true GW signal.
One of the main challenges across many GW data analysis problems, including
the SmBBH search, is the extremely high computational cost associated with
deterministic methods for solving the global optimization over $\Theta$
required in the GLRT.

It is straightforward to show that the response $h(t;\Theta)$ can be written
as [61, 62]

$h(t,\Theta)=\sum_{i=1}^{4}a_{i}A^{i}(t;\Theta_{I})\,,$ (5)

where $\Theta_{E}=\\{a_{i}\\}$, which are time-independent reparameterizations
of $\\{D_{L},\psi,\iota,\phi_{c}\\}\subset\Theta$, comprise the set of
extrinsic parameters, while
$\Theta_{I}=\\{t_{c},M_{c},\delta_{\mu},\theta,\phi,e_{0}\\}$ form the set of
intrinsic parameters. This allows $\Lambda_{G}$ to be expressed as

$\Lambda_{G}=\max_{\Theta_{I}}\left(\max_{\Theta_{E}}\Lambda(\Theta|s)\right)\,.$
(6)

One can analytically maximize $\Lambda(\Theta|s)$ over $\Theta_{E}$ to get the
fully-coherent $\mathcal{F}$-statistic [62, 61],

$\mathcal{F}(\Theta_{I}|s)=\max_{\Theta_{E}}\Lambda(\Theta|s)=\frac{1}{2}N^{T}(\Theta_{I})M^{-1}N(\Theta_{I})\,,$
(7)

where,

$\displaystyle N^{i}(\Theta_{I})$ $\displaystyle=$ $\displaystyle\langle
s|A^{i}(\Theta_{I})\rangle\,,$ (8) $\displaystyle M^{ij}$ $\displaystyle=$
$\displaystyle\langle A^{i}|A^{j}\rangle\,,$ (9) $\displaystyle a$
$\displaystyle=$ $\displaystyle N^{T}M^{-1}\,.$ (10)

The dimensional reduction above from $10$ to $6$ significantly reduces the
difficulty of the global optimization problem. Once the MLE estimates of
$\Theta_{I}$ are obtained, the corresponding MLE estimates of $\Theta_{E}$ are
obtained as shown in Appendix B.

As mentioned earlier, the global optimization of the fully-coherent
$\mathcal{F}$-statistic over the full astrophysical search range in
$\Theta_{I}$ space, called a fully-coherent search, is computationally
infeasible if carried out deterministically using a grid of points in this
space. Stochastic optimization methods are also likely to fail due to the
needle-in-the-haystack problem posed by the extremely small footprint of high
values around the global maximum compared to the required search range. A
practical solution is to use a hierarchical approach in which the search space
is first narrowed down to some promising regions in $\Theta_{I}$ space, making
follow up fully-coherent searches in each region computationally feasible. We
follow this approach using a semi-coherent version of the
$\mathcal{F}$-statistic as described below.

###  II.2 Semi-coherent $\mathcal{F}$-statistic

The semi-coherent step is inspired by the principal idea behind current
continuous wave searches using ground-based GW detectors [57, 58, 59]. This is
the partitioning of the time domain GW strain data from a detector into
smaller segments that are individually subjected to fully coherent searches.
Instead of finding just the global optimum in each search, a set of
significant locations in $\Theta_{I}$ are found. Different schemes have been
proposed and implemented [68, 69, 70, 71] for combining these locations into
candidates for follow up with fully-coherent searches using longer segment
lengths. It should be noted that there are alternative approaches [72]
possible in how the initial step of the hierarchy is implemented that do not
use the $\mathcal{F}$-statistic at all. To limit computational costs, the
threshold in the initial step of the above types of hierarchical approaches
must be set such that the number of false alarms that need to be followed up
becomes manageable. As such, semi-coherent searches generally incur a trade-
off between available computational resources and loss in detection
sensitivity.

The primary motivation for using time-segmentation in ground-based continuous
wave searches is the high cost of computing Discrete Fourier Transforms (DFTs)
due to the large sampling rate of the data. However, this is not an issue for
data from space-based detectors since the number of data samples from the
entire mission lifetime of LISA would be $\sim 10^{8}$, which is just about
the size of a ten hours segment from a single LIGO detector sampled at $2048$
Hz. Instead, for space-based detectors, the hierarchical approach can be
implemented using frequency domain partitioning and this is the approach we
follow here.

Let $\widetilde{\mathcal{F}}(\Theta_{I}|s_{l})$ denote the
$\mathcal{F}$-statistic function computed using data $s_{l}$ that is band-
limited to $[f_{l}^{(0)},f_{l}^{(1)}]$ and let

$\displaystyle\widehat{\mathcal{F}}(\Theta_{I})$ $\displaystyle=$
$\displaystyle\frac{1}{N_{\text{seg}}}\sum_{l=1}^{N_{\text{seg}}}\widetilde{\mathcal{F}}_{l}(\Theta_{I})\;,$
(11) $\displaystyle=$ $\displaystyle\frac{1}{N_{\text{seg}}}\sum_{l=1}^{N_{\rm
seg}}\max_{\Theta_{E}}\Lambda(\Theta|s_{l})$

define the semi-coherent $\mathcal{F}$-statistic, in which $N_{\rm seg}$ is
the number of segments. Since the inner product can be computed in the
frequency domain following Eq. 4, $\Lambda(\Theta_{I}|s_{l})$ is obtained by
simply redefining the inner product for segment $l$ as,

$\langle
s,h\rangle_{l}=4\,\text{Re}\int_{f_{l}^{(0)}}^{f_{l}^{(1)}}\frac{\widetilde{s}(f)\cdot\widetilde{h}^{*}(f)}{S_{n}(f)}df\,.$
(12)

For continuous GW searches in ground-based detectors, where the signal is
nearly monochromatic, the SNR and the number of wave cycles are nearly evenly
distributed across time and, hence, the semi-coherent search can use uniform
partitioning of the data in the time domain. In contrast, this is not true for
SmBBH signals since the GW power is neither uniformly distributed in time nor
in frequency since these are evolving sources. One possible approach is to
determine the frequency-domain segment boundaries based on the SNR in each
segment. However, as mentioned in [60], this partitioning scheme introduces a
dependence of the segment boundaries on the template parameters, necessitating
a recalculation of the segment limits for each evaluation of $\widehat{F}$.
This can increase program complexity and computation time. Therefore, we opt
to compute $f_{l}^{(1)}$ and $f_{l}^{(0)}$ by fixing the number of segments
for all templates but keeping $f_{1}^{(0)}$ and $f_{N_{\rm seg}}^{(1)}$
dependent on some of the intrinsic parameters. In our implementation, these
limits are set to be the lowest and highest instantaneous frequencies
occurring during the observation period in the template associated with
$\Theta_{I}$.

Unlike the fully coherent case, the calculation of the semi-coherent
$\widehat{\mathcal{F}}(\Theta_{I})$ allows different values of the extrinsic
parameters $\Theta_{I}$ for each segment. This reduces the sensitivity of the
detection statistic to the intrinsic parameters, and the SNR in each segment
is now lowered. However, the profile of the $\mathcal{F}$-statistic is
significantly broadened around the global maximum of the fully-coherent
$\mathcal{F}$-statistic, which reduces the difficulty of the search
considerably. Fig. 1 illustrates, in both the absence and presence of noise,
the profiles of the semi-coherent $\mathcal{F}$-statistic with different
$N_{\rm seg}$ as a function of the chirp mass $M_{c}$. Here, other parameters
are set to their true values. Note that $N_{\rm seg}=1$ is equivalent to the
case of the fully-coherent $\mathcal{F}$-statistic.

Figure 1:  Comparison of the semi-coherent $\mathcal{F}$-statistics
($N_{\text{seg}}>1$) with the fully-coherent $\mathcal{F}$-statistic
($N_{\text{seg}}=1$) under different conditions. The true value of the chirp
mass, $M_{c}=70.3032M_{\odot}$, is marked with a red dashed line, while the
other parameters are fixed at their true values, and the SNR is 11.3. In the
left panel, one can see that, under noiseless condition, the peaks for the
semi-coherent $\mathcal{F}$-statistics are significantly broadened compared to
the fully-coherent $\mathcal{F}$-statistic, which can mitigate the difficulty
of the search for the global maximum. However, under noisy conditions, as
shown in the right panel, as $N_{\text{seg}}$ increases, the peaks become less
distinct relative to the fluctuations caused by noise.

Fig. 1(a) shows how, in the absence of noise, the semi-coherent method expands
the width of the peak and absorbs the sharp secondary peaks near the true
value of $M_{c}$ into the main peak. This makes it easier to search for the
peak. However, as shown in Fig. 1(b), the peak also becomes less prominent
compared to the fluctuations away from the peak caused by noise. Thus,
adopting an appropriate $N_{\rm seg}$ according to the type of the GW signal
is a key issue in a semi-coherent search.

##  III Hierarchical search pipeline

In this section, we provide a brief overview of Particle Swarm Optimization
(PSO) and discuss the specific implementation used in this work. This is
followed by a description of the search pipeline that uses PSO in multiple
stages of a hierarchical search pipeline.

###  III.1 PSO

PSO is an iterative stochastic search method for the global maximum of a
function $f(x)$, $x\in\mathbb{R}^{N}$, called the fitness function, over a
specified compact subset $D\subset\mathbb{R}^{N}$ called the search space. [In
our case, the fitness function is the semi-coherent $\mathcal{F}$-statistic
(c.f., Eq. 11) and the search space is the space of intrinsic parameters
$\Theta_{I}$.] In each iteration, the fitness is sampled at multiple locations
and the locations are updated for the next iteration based on the sampled
fitness values. The sampled locations, which are typically fixed in number,
are called particles in PSO, with the rules used for updating their positions,
called the dynamical equations, stated below. The collection of particles is
called a swarm.

Let $x_{i,j}[k]$ denote component $j$ of the position vector $x_{i}$ of
particle $i\leq N_{\rm part}$ at iteration $k$. Attached to each particle is a
displacement vector $v_{i}[k]$ (called velocity in PSO terminology) with
component $j$ denoted as $v_{i,j}[k]$. The position and velocity are updated
as,

$x_{i,j}[k+1]=x_{i,j}[k]+\min(v_{i,j}[k+1],v_{\rm max})\,,$ (13)
$\displaystyle v_{i,j}[k+1]$
$\displaystyle=w[k]v_{i,j}[k]+c_{1}r_{1}(p_{i,j}[k]-x_{i,j}[k])$ (14)
$\displaystyle+c_{2}r_{2}(l_{i,j}[k]-x_{i,j}[k])\,,$

where $v_{\rm max}$ limits the maximum step size in any dimension, $w[k]$
denotes a deterministically decaying number called the inertia weight,
$c_{1,2}$ are called acceleration constants, $r_{1,2}\sim U([0,1])$ are
independent random variables with uniform distribution over $[0,1]$,
$p_{i}[k]$ (personal best) is the location with the best fitness in the
history of particle $i$, and $l_{i}[k]$ (local best) is the best location
found in the history of its local neighborhood. (When the local neighborhood
includes all the other particles in the swarm, $l_{i}[k]$ becomes the global
best, $g[k]$.) An example of a neighborhood definition is the ring topology
[73] in which particle indices are arranged on a circle and a subset of
consecutive indices constitute a neighborhood.

Each of the three terms in Eq. 14 has a special role. The term with the
inertial weight, called the inertia term, causes the particle to move past its
current location irrespective of its fitness value. This allows the particle
to escape local maxima and contributes to the exploration of the search space.
The second term, called the cognitive term, that depends on the personal best
can be interpreted as a random force that tries to attract the particle back
to a good location that it has already found. The last term, called the social
term, that depends on the local best is also a random force that seeks to pull
the particle to a good location found by its neighbors. The cognitive and
social terms promote exploration around previously found good solutions, a
stage of the search called exploitation, but the randomness in these terms
prevents too quick a convergence. Thus, the dynamical equations create a
trade-off between the exploration and exploitation behavior of the swarm, with
the former dominating in the beginning. The use of a local instead of a global
best extends the exploration phase by slowing down information flow within the
swarm. The inertia weight $w[k]$ is typically allowed to decay linearly with
iterations, providing yet another control mechanism for the shift in swarm
behavior from exploration to exploitation.

Most PSO implementations work best when the search space $D$ is a hypercube,
with component $j$ of $x$ limited to $[a_{j},b_{j}]$. The iterations are
commonly initiated with random locations $x_{i,j}[0]\sim U([a_{j},b_{j}])$ and
velocities $v_{i,j}\sim U[a_{j}-x_{i,j},b_{j}-x_{i,j}]$. While these initial
conditions guarantee that particles will be localized within $D$ in the first
position update, it is possible for them to escape $D$ in subsequent
iterations. Therefore, boundary conditions need to be specified to handle such
particles. Among the many boundary conditions proposed in the literature [74],
we choose the let-them-fly condition in which nothing is done to the particles
except for setting their fitness to $-\infty$ (in a maximization problem) when
they escape $D$. Similarly, among the diverse termination conditions proposed
for PSO, we pick the simplest one in which termination happens once a
specified number of iterations is completed. The final solution for the
maximization problem is the global best location and the corresponding fitness
value.

Like most other stochastic heuristics for global optimization, there is no
guarantee of convergence (even asymptotically) to the global maximum in PSO
but it can be tuned to provide an acceptable probability, $P_{\rm success}$,
of successful convergence to a given region containing the global maximum.
This probability can be boosted exponentially by running $N_{\rm runs}$
independent runs of PSO, with statistically independent pseudorandom streams
for $r_{1,2}$ and initialization, and picking the best solution across all the
runs. The resulting probability of success becomes $1-(1-P_{\rm
success})^{N_{\rm runs}}$, which increases rapidly with $N_{\rm runs}$ even
for moderate $P_{\rm success}$.

Much of the overview provided above is common to several applications of PSO
in GW data analysis. In fact, we have adopted the same settings for most of
the parameters above as given in [51] except for $N_{\rm part}=128$, which is
significantly higher than the usual choice of $\approx 40$ particles, and
different $N_{\rm runs}$ at different stages of the hierarchical search. The
large number of particles used here follows from the specific implementation
of GPU acceleration used in our code, which requires $N_{\rm part}$ to be a
power of $2$, and the fact that lower $N_{\rm part}$ values of $32$ and $64$
did not yield a high probability of successful convergence. As described
below, the number of PSO iterations used in the termination criterion varies
for the different stages of the hierarchy. The remaining PSO parameters are
set to widely used values [73]: $c_{1}=c_{2}=2.0$ and $v_{\rm max}=0.5$.

###  III.2 Hierarchical search

Despite the significant broadening of the peak in the fitness function, as
shown in Fig. 1, it is still much smaller in width than the full search range
used in our analysis. In addition, the semi-coherent fitness function is very
rugged with multiple local maxima induced by noise in the data. Therefore, a
direct search for the peak of the semi-coherent $\mathcal{F}$-statistic using
a single PSO run has a low probability of success in the presence of noise. To
cope with this problem, we divide the search space into smaller regions along
certain dimensions, search over the regions independently, and propagate their
results in a hierarchical manner through multiple stages in order to
progressively narrow down the search around the global maximum of the semi-
coherent $\mathcal{F}$-Statistic. The details of this hierarchical pipeline
are provided below and illustrated schematically in Fig. 2.

Figure 2:  Flowchart of the hierarchical stages. Each stage is shown by a
stack of colored boxes (the boxes are colored with alternating colors for
visual clarity). Each stack has two axes: the horizontal and vertical axes are
labeled by parameters for which the search range is not subdivided or divided,
respectively. Each colored box schematically represents one partition of the
search range for the latter set of parameters. In stage–1, starting at the
left, the search range of $t_{c}$ is uniformly divided into 24 contiguous
intervals. In stage–2, based on the results of stage–1, the search ranges of
both $t_{c}$ and $M_{c}$ are narrowed down, while the search ranges of other
parameters remain unchanged. In stage–3, the search ranges of all parameters
except $\delta_{\mu}$ are narrowed down. Each colored box is followed by a box
showing the number of PSO runs used in that restricted search range. Each box
labeled MAX represents the operation of selecting the best solution out of the
PSO runs connected to it. Finally, in stage–3, the two search ranges differ in
the range of $\theta$, with the search ranges over all the other parameters
kept identical.

Stage–1: The search range of the parameter $t_{c}$ is divided into $24$
equally wide and contiguous intervals. In each interval, $24$ independent PSO
runs are carried out to search for the global maximum of the semi-coherent
$\mathcal{F}$-statistic using $N_{\rm seg}=16$ segments. The best solution
among the $24$ runs in each region is then recorded, producing a total of $24$
candidate locations for the global maximum in the full parameter space. The
number of iterations in this step is $8000$ for each PSO run.

Stage–$2$: In this stage, the $24$ candidate locations obtained in stage 1
form the centers for new search regions that are restricted in both $t_{c}$
and $M_{c}$. Based on our experience with different settings, we set the size
of each restricted region to $0.01$ yr along $t_{c}$ and $0.6$ $M_{\odot}$
along $M_{c}$, while keeping the search ranges for all other parameters the
same as in stage $1$. Each of the $24$ regions is now searched with $8$ PSO
runs and the single best solution from all these $24\times 8=192$ runs is
propagated to the next stage. The number of PSO iterations and
$N_{\text{seg}}$ for semi-coherent $\mathcal{F}$-statistic remain the same as
in stage $1$. The single solution provides a candidate location in a search
region that is restricted in $t_{c}$ and $M_{c}$ but unrestricted in all other
parameters.

Stage–3: In this stage, the single solution obtained in the previous stage
forms the center of narrower search ranges with $0.004$ yr for $t_{c}$, $0.3$
$M_{\odot}$ for $M_{c}$, $0.4$ rad for $\theta$, $0.4$ rad for $\phi$, and
$0.006$ for $e_{0}$. As for $\delta_{\mu}$, due to the relative insensitivity
of the semi-coherent $\mathcal{F}$-statistic to this parameter, we continue to
maintain its original search range. Besides restricting the search ranges, we
also create a search region centered on $\pi-\theta$, due to a degeneracy of
the $\mathcal{F}$-statistic in $\theta$. Each of the two search spaces is now
subjected to $48$ PSO runs with 15,000 iterations per run. At this point,
$N_{\text{seg}}$ is set to 1, which means that this stage uses the fully-
coherent $\mathcal{F}$-statistic. The best solution across the two search
regions provides the final estimate of the intrinsic parameters that are then
used to obtain the corresponding estimates of the extrinsic parameters
following the expressions given in Appendix B.

While we restrict the search ranges for only $t_{c}$ and $M_{c}$ in the first
two stages, it is conceivable to have a more general scheme in which the
search ranges of additional parameters are restricted. However, our experience
shows that this is redundant in the first two stages because the semi-coherent
$\mathcal{F}$-statistic is most sensitive to these two parameters and
relatively insensitive to the others. Moreover, subdividing the ranges of the
parameters increases the computational cost as the number of such parameters
goes up. Thus, subdividing the ranges of only $t_{c}$ and $M_{c}$ offers a
compromise between effectively locating the interesting search region for the
latter stages and the overall computational cost of the search.

It is worth noting that we do not reduce the number of search ranges when
going from stage–1 to stage–2. This is because, in the stage-1 search, there
is a high probability that secondary maxima induced by noise in search regions
far from the one containing the true parameter have larger values. Therefore,
confining the stage-2 search to a subset of the stage-1 search ranges incurs
the risk of prematurely converging to secondary maxima. It is preferable to
further refine the solution in each search range, as done in stage–2, before
deciding to pick the best solution.

##  IV Efficient computation of the semi-coherent $\mathcal{F}$ statistic

The hierarchical stages described above involve a large number of PSO runs
over independent search regions, with each run requiring a large number of
iterations. To mitigate the computational costs involved, we use an efficient
calculation of the frequency domain TDI response to a plane GW and implement
GPU acceleration of parallelizable steps.

The computationally most expensive step in PSO is the evaluation of the semi-
coherent $\mathcal{F}$-statistic, $\widehat{F}$, for each particle.
Calculating $\widehat{F}$ involves computing the noise-weighted inner product
between the templates and the data, which requires Fourier transforms and
consumes a significant amount of time. According to [75], for sources where
the timescale of the gravitational wave radiation reaction is much smaller
than the timescale of the detector modulation, it is possible to directly
obtain their frequency-domain responses.

A TDI combination is constructed as a sum of time-delayed single-arm
responses. The general form of the GW-induced single-arm response is given by:

$x(t)=F(t)h(t+d(t))\,.$ (15)

Here, $F(t)$ is the time-domain modulation associated with the rotating
antenna pattern of the arm and precessional effects in the source orbit
evolution. Note that the former depends on the orbit of the detector. For
LISA, the orbit used in this study is an analytic approximation, with accuracy
extended to second-order eccentricity [76]. While this provides convenience
and efficiency, it is worth noting that our work is not dependent on it.
Numerical orbits with higher precision can also be used, and the evaluation of
the detector antenna pattern at specific times requires numerical
interpolation using appropriate orders of Chebyshev polynomials, as
demonstrated in [77]. $d(t)$ represents the time delay due to the orbital
motion of the detector. A GW chirp signal, such as the one from an SmBBH
inspiral, has a well-defined instantaneous frequency and its Fourier transform
can be expressed as [78, 79, 80]:

$\widetilde{h}(f)=A(f)e^{i\Phi(f)}\,,$ (16)

where $A(f)$ and $\Phi(f)$ can be obtained using the stationary phase
approximation. Given the much slower timescale of the modulation and time
delay variations compared to the instantaneous period of the GW signal, the
Fourier transform of $x(t)$ is well-approximated by [75]:

$\widetilde{x}(f)=T(f)\widetilde{h}(f)\,,$ (17)

where

$T(f)=\sum_{k\geq
0}\frac{(-i)^{k}}{k!}(T_{Ak})^{k}F_{T_{f,\epsilon}}\left[\frac{d^{k}}{d\tau^{k}}\left(\frac{F(\tau-d(\tau))}{1+\dot{d}(\tau)}e^{2i\pi
fd(\tau)(1-\dot{d}(\tau))}\right)\right](t_{f})\,.$ (18)

Here, ($a_{N,k}^{\epsilon}$ are defined in Appendix 2 of [75])

$T_{A1}=\frac{7}{6}\frac{1}{2\pi f}\,,$ (19)
$t_{f}=\frac{1}{2\pi}\frac{d\Phi(f)}{df}\,,$ (20)
$T_{f}=\frac{1}{4\pi}\left|\frac{d^{2}\Phi}{df^{2}}\right|\,,$ (21)
$\epsilon=-\text{sgn}\left(\frac{d^{2}\Phi}{df^{2}}\right)=-1\,,$ (22)
$F_{T_{f,\epsilon}}^{N}(F)(t)\equiv\frac{1}{2}\sum_{k=0}^{N}a_{N,k}^{\epsilon}\left(F(t+kT_{f})+F(t-kT_{f})\right)\,.$
(23)

In this work, we set both $k$ and $N$ to be zero, and neglecting the higher-
order terms of $d(t)$ in the frequency-domain response, the matching with that
obtained by Fourier-transform of the time-domain signal exceeds 0.999, which
meets our requirements. Using the frequency-domain response to compute
$\widehat{F}$ not only directly simplifies the computation — no Fourier
transforms are required at this point — but also provides the conditions for
deep parallelization of the program.

In the hierarchical method presented in Sec. III, there are several nested
parallelizable steps. First, each PSO run can be executed independently.
Secondly, within each PSO run, the computation of the semi-coherent
$\mathcal{F}$-statistic for each particle can be carried out in parallel.
Finally, in the frequency-domain response, the computations for each frequency
can be parallelized. We implement these layers of parallelization using a
nested set of hardware elements, starting with nodes of a cluster at the top
level, CPU cores at the intermediate level, and GPUs at the deepest level. At
present, our parallelization scheme allows the code implementing the
hierarchical search described in Sec. III to perform a complete analysis of a
single data realization in $\approx 2.5$ days using $8$ NVIDIA A100 GPUs
attached to a single compute node with $48$ CPU cores.

##  V Results

In this section, we discuss the results obtained by applying our pipeline to
different sets of simulated data.

###  V.1 Noiseless data

First, we validate the hierarchical search method on noiseless data. The
values of the parameters and corresponding search ranges for the simulated GW
signals are given in Table 1. Since there is no noise in the data, we can
adjust some of the configuration parameters to increase the execution speed.
For example, we set $N_{\text{seg}}=128$ in both stage–1 and stage–2, which
significantly reduces the search complexity by broadening the global maximum
of the semi-coherent $\mathcal{F}$-statistic while not incurring any risk of
noise-induced secondary maxima. At the same time, we reduce the number of PSO
iterations to 2000 for the first two stages and 5000 for the stage–3. Finally,
the number of PSO runs per search region was reduced substantially. These
adjustments reduced the runtime of the program by a factor of $10$ compared to
the fiducial settings given in Sec. III.

Table 1: The parameters used in the search of noiseless data and their search ranges. The signal amplitude is chosen such that the SNR relative to the instrumental noise has the value shown here. $D_{L}$, $\iota$, $\psi$, and $\phi_{c}$ are extrinsic parameters, whose estimated values are obtained from the relationship given in the Appendix B. Hence, they do not have an associated search range. Parameter | True value | Search range
---|---|---
$t_{c}/{\rm yr}$ | 3.5000 | [3.3, 3.8]
$M_{c}/M_{\odot}$ | 70.3032 | [20, 80]
$\delta_{\mu}$ | 0.1585 | [0, 0.9]
$\theta/{\rm rad}$ | 0.7854 | [0, $\pi$]
$\phi/{\rm rad}$ | 2.1216 | [0, $2\pi$]
$e_{0}(f_{\rm gw}=0.01\,\text{Hz})$ | 0.0100 | [0, 0.05]
$D_{L}/{\rm Mpc}$ | 920.9084 | /
$\iota/{\rm rad}$ | 0.5329 | /
$\psi/{\rm rad}$ | 1.2501 | /
$\phi_{c}/{\rm rad}$ | 0.8419 | /
${\rm SNR}$ | 11.3137 | /

Table 2 presents a comparison of the true values with the estimated ones for
noiseless data. It can be seen that, except for $\delta_{\mu}$, the intrinsic
parameters are recovered with very small relative errors. Part of the error in
$\delta_{\mu}$ can be attributed to degeneracies between $\delta_{\mu}$ and
several other parameters in the $\mathcal{F}$-statistic. To test this
hypothesis, we reran the analysis fixing $\delta_{\mu}$ to its true value. The
corresponding results displayed in Table 2 show that errors in almost all the
parameters are reduced with some of the reductions being significant. This
shows that allowing $\delta_{\mu}$ to be a free parameter causes errors in
other parameters to be compensated by an error in $\delta_{\mu}$. Table 2 also
contains the estimated values of the extrinsic parameters calculated by the
procedure in Appendix B. It can be seen that the relative errors in these
parameters are orders of magnitude larger than those in the intrinsic
parameters. We believe that this is due to some degree of degeneracy within
the extrinsic parameters. In the future, using more accurate waveforms that
include higher-order modes should improve the accuracy of estimation for the
extrinsic parameters.

Thus, the analysis of noiseless data validates our code and allows us to
proceed to the analysis of realistic noisy data.

Table 2: Results of the search in noiseless data. It can be seen that the intrinsic parameters, except $\delta_{\mu}$, are close to their true values. However, there is a considerable error in the extrinsic parameters. This is due to the problem of degeneracy in the GW waveform. Even in a noise-free scenario, obtaining a highly accurate value for these parameters remains challenging. The results of the search when $\delta_{\mu}$ is fixed at the true value are also given here, and it can be seen that the error is basically further reduced. Parameter | True value | Estimated value | Relative error | Estimated value ($\delta_{\mu}$ fixed) | Relative error
---|---|---|---|---|---
$t_{c}/\rm yr$ | 3.5000000 | 3.5000029 | 0.0000836% | 3.4999994 | 0.0000172%
$M_{c}/M_{\odot}$ | 70.30316 | 70.30662 | 0.00492% | 70.30319 | 0.0000436%
$\delta_{\mu}$ | 0.1585 | 0.4784 | 201% | / | /
$\theta/\rm rad$ | 0.78540 | 0.78554 | 0.0174% | 0.78535 | 0.00621%
$\phi/\rm rad$ | 2.12160 | 2.12162 | 0.00128% | 2.1214 | 0.00788%
$e_{0}(f_{\rm gw}=0.01\text{Hz})$ | 0.01000000 | 0.00998489 | 0.151% | 0.00999905 | 0.00947%
$D_{L}/\rm Mpc$ | 920.9084 | 880.9283 | 4.34% | 895.3753 | 2.77%
$\iota/\rm rad$ | 0.5329 | 0.6088 | 14.2% | 0.5826 | 9.33%
$\psi/\rm rad$ | 1.2501 | 1.2846 | 2.76% | 1.2831 | 2.64%
$\phi_{c}/\rm rad$ | 0.8419 | 0.2348 | 72.1% | 0.7808 | 7.26%

###  V.2 Noisy data

Within the analysis of noisy data, we considered three separate cases, labeled
as cases A, B, and C. Table 3 provides the values of the true signal
parameters in the three cases. The search ranges are identical for all three
cases. In the following, we discuss the results from these cases in sequence.

Table 3: The parameters in Case A, Case B, and Case C and their respective search ranges. Parameter | Case A | Case B | Case C | Search range
---|---|---|---|---
$t_{c}/\rm yr$ | 3.5000 | 3.5000 | 3.5000 | [3.3, 3.8]
$M_{c}/M_{\odot}$ | 70.3032 | 54.3572 | 28.1923 | [20, 80]
$\delta_{\mu}$ | 0.1585 | 0.0400 | 0.2485 | [0, 0.9]
$\theta/\rm rad$ | 0.7854 | 0.7854 | 0.7854 | [0, $\pi$]
$\phi/\rm rad$ | 2.1216 | 2.1216 | 2.1216 | [0, $2\pi$]
$e_{0}(f_{\rm gw}=0.01\text{Hz})$ | 0.0100 | 0.0100 | 0.0100 | [0, 0.05]
$D_{L}/\rm Mpc$ | 920.9084 | 603.5043 | 123.9608 | /
$\iota/\rm rad$ | 0.5329 | 0.5329 | 0.5329 | /
$\psi/\rm rad$ | 1.2501 | 1.2501 | 1.2501 | /
$\phi_{c}/\rm rad$ | 0.8419 | 0.8419 | 0.8419 | /
$\rm SNR$ | 11.3137 | 11.1723 | 14.1421 | /

The functioning of the hierarchical search is illustrated in Figures 3 to 5
pertaining to Case A. Each figure shows the progression in the localization of
the estimated parameter around the true value as one proceeds from stage–1 to
stage–3. We see that there is a significant improvement for all parameters
except $\delta_{\mu}$. The figures for stage–1 show clearly the effect of
secondary maxima induced by noise. For the $t_{c}$ parameter, the secondary
maxima are often larger than the semi-coherent $\mathcal{F}$-statistic values
in the search range containing the true value of this parameter. It is worth
noting that the latter may be just some values found by PSO, not necessarily
local maxima at this range. They could be close to secondary maxima, but they
could simply result from noise. As discussed earlier in Sec. III, confining
the stage–2 search to a subset of stage–1 search regions based on the semi-
coherent $\mathcal{F}$-statistic values could, therefore, lead to premature
convergence. The results for stage–2 show the importance of restricting the
search range in $M_{c}$ along with $t_{c}$. We see that, this leads to a
significant improvement in the semi-coherent $\mathcal{F}$-statistic value and
that it can often exceed the value at the true location. This is a positive
indicator that PSO was able to successfully converge to the global maximum.
Finally, we see how the search range is narrowed down by several orders of
magnitude at the beginning of stage–3, thus allowing the fully-coherent
$\mathcal{F}$-statistic to be computed successfully.

Figure 3:  Part of the PSO results about $t_{c}$ (top row) and $M_{c}$ (bottom
row). In each panel, the red circle marks the true values of the parameters
and their fitness, dots with different colors come from different search
ranges of $t_{c}$. The black dots represent the PSO results from the range of
$t_{c}$ that encloses its true value. From left to right are stage–1, stage–2,
and stage–3. Figure 4:  As in Fig. 3, the PSO results about $\delta_{\mu}$ and
$\theta$. Figure 5:  As in Fig. 3, the PSO results about $\phi$ and $e_{0}$.

Tables 4 to 6 present the estimated parameters and their relative errors for
the three cases. It can be seen that the method was able to estimate the
intrinsic parameters with very high accuracy. For example, in Case A, the
relative error in $M_{c}$ is less than $0.01\%$ and $t_{c}$ is determined with
an accuracy of $\approx 100$ sec. The sky localization is accurate to within
$0.2~{}{\rm deg}^{2}$, and the relative errors for eccentricity and chirp mass
are less than $1\%$ in all cases.

Table 4: Case A search results. The leftmost column lists the true parameters of the injected signal. The second column lists the estimated values found from stage–3 and the third column presents the relative error between the true and estimated values of each parameter. Parameter | True value | Estimated value | Relative error
---|---|---|---
$t_{c}/\rm yr$ | 3.500000 | 3.499997 | 0.0000899%
$M_{c}/M_{\odot}$ | 70.3032 | 70.3101 | 0.00986%
$\delta_{\mu}$ | 0.1585 | 0.5797 | 265%
$\theta/\rm rad$ | 0.7854 | 0.7909 | 0.696%
$\phi/\rm rad$ | 2.1216 | 2.1291 | 0.354%
$e_{0}(f_{\rm gw}=0.01\text{Hz})$ | 0.010000 | 0.009953 | 0.473%
$D_{L}/\rm Mpc$ | 920.9084 | 397.4851 | 56.8%
$\iota/\rm rad$ | 0.5329 | 1.2719 | 139%
$\psi/\rm rad$ | 1.2501 | -0.5077 | 141%
$\phi_{c}/\rm rad$ | 0.8419 | 0.1554 | 81.5%
Table 5: Case B search results. The leftmost column lists the true parameters of the injected signal. The second column lists the estimated values found from stage–3 and the third column presents the relative error between the true and estimated values of each parameter.  Parameter | True value | Estimated value | Relative error
---|---|---|---
$t_{c}/\rm yr$ | 3.5000000 | 3.5000142 | 0.000405%
$M_{c}/M_{\odot}$ | 54.3572 | 54.3617 | 0.00838%
$\delta_{\mu}$ | 0.0400 | 0.5721 | 1330%
$\theta/\rm rad$ | 0.7854 | 0.7898 | 0.558%
$\phi/\rm rad$ | 2.1216 | 2.1192 | 0.112%
$e_{0}(f_{\rm gw}=0.01\text{Hz})$ | 0.01000 | 0.009970 | 0.300%
$D_{L}/\rm Mpc$ | 603.5043 | 513.3089 | 14.9%
$\iota/\rm rad$ | 0.5329 | 0.7259 | 36.2%
$\psi/\rm rad$ | 1.2501 | 0.4386 | 64.9%
$\phi_{c}/\rm rad$ | 0.8419 | 0.4058 | 51.8%
Table 6: Case C search results. The leftmost column lists the true parameters of the injected signal. The second column lists the estimated values found from stage–3 and the third column presents the relative error between the true and estimated values of each parameter. Parameter | True value | Estimated value | Relative error
---|---|---|---
$t_{c}/\rm yr$ | 3.500000 | 3.5000039 | 0.0001126%
$M_{c}/M_{\odot}$ | 28.1923 | 28.1932 | 0.00315%
$\delta_{\mu}$ | 0.07692 | 0.41195 | 4355%
$\theta/\rm rad$ | 0.7854 | 0.7864 | 0.129%
$\phi/\rm rad$ | 2.1216 | 2.1228 | 0.0575%
$e_{0}(f_{\rm gw}=0.01\text{Hz})$ | 0.0100 | 0.009961 | 0.388%
$D_{L}/\rm Mpc$ | 123.9608 | 74.1714 | 40.2%
$\iota/\rm rad$ | 0.5329 | 1.1470 | 115%
$\psi/\rm rad$ | 1.2501 | -0.3975 | 132%
$\phi_{c}/\rm rad$ | 0.8419 | 0.7437 | 11.7%

Case B and Case C differ from Case A primarily in the injected values of the
chirp mass. This design allows us to validate our method at different points
within its search range. The parameters $t_{c}$, $\theta$, $\phi$, and
$e_{0}$, as well as the extrinsic parameters (except $D_{L}$), remain the same
as in Case A. The parameter $D_{L}$ is adjusted to achieve the desired SNRs of
the injected signals. As shown in Tables 5 to 6, the relative errors in the
intrinsic parameters are consistent with those in Case A, despite the noise
realizations being independent draws based on the LISA instrumental noise PSD.

Moreover, the relative errors in $D_{L}$ exceed $10\%$ across all three cases
with relatively low SNRs, posing a challenge for distance measurement of these
sources using GW alone. However, this may be mitigated by the high estimation
accuracy achieved in $\theta$ and $\phi$, making it possible for
electromagnetic follow-up to pinpoint their host galaxies.

##  VI Conclusion

We have presented a hierarchical search method for SmBBH signals in noisy LISA
data that shows promising performance for an SNR range of $\approx 11$ to
$\approx 14$ for $2.5$ yr of a single TDI combination. The method uses a
hierarchy consisting of partitioned parameter search ranges and semi-coherent
$\mathcal{F}$-statistic with different numbers of segments. The search for the
global maximum of the semi-coherent $\mathcal{F}$-statistic in each search
range is performed using PSO. By employing frequency-domain responses and GPU
parallel computing to accelerate the code, we were able to efficiently search
for single GW signals within a reasonable timeframe of $2.5$ days. We
validated the reliability of this pipeline using simulated data containing
noise based on instrumental design power spectral density of LISA. Our results
demonstrate that it is possible to address the data analysis challenges posed
by SmBBH signals and produce highly accurate estimates of intrinsic
parameters. In particular, the time of coalescence can be localized within
$\approx 100$ sec and the sky location can be estimated within $\approx 0.2$
${\rm deg}^{2}$, which is more than adequate to enable multi-messenger and
multi-wavelength astronomy for SmBBH sources.

LISA can also measure the orbital eccentricity during inspirals of the SmBBH
with a relative error of 1$\%$, providing the ability to distinguish between
their origin from binary star evolution or the dynamical interactions in dense
star clusters. However, our waveform model does not include higher-order modes
that would be present in a system with high orbital eccentricity. In addition,
our waveform model does not include the spins of the individual black holes,
which may be important in some cases. We will address these limitations,
including the use of multiple TDI combinations, in future studies.

Finally, the signals we simulated all merged shortly outside the observation
time window, whereas in reality, some sources may merge within the window,
while others may merge long after the window. This implies that the search
range for the merger time should be expanded, and the problem of multiple
sources overlapping should also be addressed.

###### Acknowledgements.

Y.W. gratefully acknowledges support from the National Key Research and
Development Program of China (No. 2022YFC2205201), the National Natural
Science Foundation of China (NSFC) under Grant No. 11973024, Major Science and
Technology Program of Xinjiang Uygur Autonomous Region (No. 2022A03013-4), and
Guangdong Major Project of Basic and Applied Basic Research (Grant No.
2019B030302001). S.D.M. is supported by U.S. National Science Foundation (NSF)
Grant No. PHY-2207935. We acknowledge the High Performance Computing Platform
at Huazhong University of Science and Technology for providing computational
resources.

##  Appendix A Gravitational Waveforms

In this paper, we employ the restricted post-Newtonian waveforms in the
frequency domain, where “restricted” means that we neglect higher-order
amplitude modulation and focus only on the post-Newtonian effects in the
phase:

$h(f)=Af^{-7/6}e^{-i\Psi(f)}$ (24)

where (for $f>0)$:

$A=-\sqrt{5/24}(M_{c}^{5/6})/(\pi^{2/3}D_{L})\,,$ (25) $\displaystyle\Psi(f)$
$\displaystyle=2\pi ft_{c}-\phi_{c}-\frac{\pi}{4}+\frac{3}{128}(M_{c}\pi
f)^{-5/3}\sum_{k=0}^{6}\alpha_{k}x^{k/2}$ (26)
$\displaystyle-\phi_{p}(t(f))-\phi_{D}(t(f))+\phi_{e}(f)\,.$

The details of the parameters, including the coefficients $\alpha_{k}$, can be
found in [81, 82, 83, 84]. In addition, for simplicity, we have omitted
certain parameters, such as the spin of the individual black holes, from the
waveforms in the current paper, resulting in a total of ten parameters
represented by $\Theta$.

##  Appendix B Relationships between extrinsic parameters and intrinsic
Parameters

In Sec. II.1, the fully-coherent $\mathcal{F}$-statistic is obtained by
maximizing the log-likelihood ratio function over $\Theta_{E}=\\{a_{i}\\}$,
which are time-independent reparameterizations of
$\\{D_{L},\psi,\iota,\phi_{c}\\}$ as follows [85]:

$\displaystyle a_{1}=\frac{1}{D_{L}}\left(\cos
2\psi\frac{1+\cos^{2}\iota}{2}\cos\phi_{c}-\sin
2\psi\cos\iota\sin\phi_{c}\right)\,,$ (27) $\displaystyle
a_{2}=-\frac{1}{D_{L}}\left(\sin
2\psi\frac{1+\cos^{2}\iota}{2}\cos\phi_{c}+\cos
2\psi\cos\iota\sin\phi_{c}\right)\,,$ $\displaystyle
a_{3}=-\frac{1}{D_{L}}\left(\cos
2\psi\frac{1+\cos^{2}\iota}{2}\sin\phi_{c}+\sin
2\psi\cos\iota\cos\phi_{c}\right)\,,$ $\displaystyle
a_{4}=\frac{1}{D_{L}}\left(\sin
2\psi\frac{1+\cos^{2}\iota}{2}\sin\phi_{c}-\cos
2\psi\cos\iota\cos\phi_{c}\right)\,.$

Once the maximum likelihood estimate of the intrinsic parameters
($\Theta_{I}$) is obtained, $a_{i}$ can be determined by Eq. 10. Subsequently,
the corresponding maximum likelihood estimate of
$\\{D_{L},\psi,\iota,\phi_{c}\\}$ can be obtained by the following
relationship,

$\displaystyle D_{L}$ $\displaystyle=\frac{2}{A},$ (28) $\displaystyle\iota$
$\displaystyle=\arccos\left(\frac{-A_{\times}}{A_{+}+\sqrt{A_{+}^{2}-A_{\times}^{2}}}\right)\,,$
$\displaystyle\psi$
$\displaystyle=\frac{1}{2}\arctan\left(\frac{A_{+}a_{4}-A_{\times}a_{1}}{-(A_{\times}a_{2}+A_{+}a_{3})}\right)\,,$
$\displaystyle\phi_{c}$
$\displaystyle=\arctan\left(\frac{c(A_{+}a_{4}-A_{\times}a_{1})}{-c(A_{\times}a_{3}+A_{+}a_{2})}\right)\,,$

where:

$\displaystyle A_{+}$
$\displaystyle=\sqrt{(a_{1}+a_{4})^{2}+(a_{2}-a_{3})^{2}}+\sqrt{(a_{1}-a_{4})^{2}+(a_{2}+a_{3})^{2}}\,,$
(29) $\displaystyle A_{\times}$
$\displaystyle=\sqrt{(a_{1}+a_{4})^{2}+(a_{2}-a_{3})^{2}}-\sqrt{(a_{1}-a_{4})^{2}+(a_{2}+a_{3})^{2}}\,,$
$\displaystyle A$ $\displaystyle=A_{+}+\sqrt{A_{+}^{2}-A_{\times}^{2}},$
$\displaystyle c$ $\displaystyle=\text{sgn}(\sin(2\psi))\,.$

## References

  * Abbott  _et al._ [2016a] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, _et al._ , Observation of gravitational waves from a binary black hole merger, Physical review letters 116, 061102 (2016a).
  * Aasi  _et al._ [2015] J. Aasi, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, _et al._ , Advanced ligo, Classical and quantum gravity 32, 074001 (2015).
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, _et al._ , Advanced Virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Somiya [2012] K. Somiya, Detector configuration of KAGRA-the Japanese cryogenic gravitational-wave detector, Classical and Quantum Gravity 29, 124007 (2012), arXiv:1111.7185 [gr-qc] .
  * Abbott  _et al._ [2016b] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, _et al._ , Binary black hole mergers in the first advanced ligo observing run, Physical Review X 6, 041015 (2016b).
  * Abbott  _et al._ [2019] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, _et al._ , Gwtc-1: a gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs, Physical Review X 9, 031040 (2019).
  * Abbott  _et al._ [2021] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, _et al._ , Gwtc-2: compact binary coalescences observed by ligo and virgo during the first half of the third observing run, Physical Review X 11, 021053 (2021).
  * Abbott  _et al._ [2023] R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, _et al._ , Gwtc-3: compact binary coalescences observed by ligo and virgo during the second part of the third observing run, Physical Review X 13, 041039 (2023), arXiv:2111.03606 [gr-qc] .
  * Punturo  _et al._ [2010] M. Punturo, M. Abernathy, and et al., The Einstein Telescope: a third-generation gravitational wave observatory, Classical and Quantum Gravity 27, 194002 (2010).
  * Srivastava  _et al._ [2022] V. Srivastava, D. Davis, K. Kuns, P. Landry, S. Ballmer, M. Evans, E. D. Hall, J. Read, and B. S. Sathyaprakash, Science-driven Tunable Design of Cosmic Explorer Detectors, Astrophys. J. 931, 22 (2022), arXiv:2201.10668 [gr-qc] .
  * Hild  _et al._ [2011] S. Hild, M. Abernathy, F. Acernese, P. Amaro-Seoane, N. Andersson, and others., Sensitivity studies for third-generation gravitational wave observatories, Classical and Quantum Gravity 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, _et al._ , Laser Interferometer Space Antenna, arXiv e-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Luo  _et al._ [2016] J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, _et al._ , TianQin: a space-borne gravitational wave detector, Classical and Quantum Gravity 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Mei  _et al._ [2021] J. Mei, Y.-Z. Bai, J. Bao, E. Barausse, L. Cai, and others., The TianQin project: Current progress on science and technology, Progress of Theoretical and Experimental Physics 2021, 05A107 (2021), arXiv:2008.10332 [gr-qc] .
  * Hu and Wu [2017] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, National Science Review 4, 685 (2017), https://academic.oup.com/nsr/article-pdf/4/5/685/31566708/nwx116.pdf .
  * Ruan  _et al._ [2020] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Taiji program: Gravitational-wave sources, International Journal of Modern Physics A 35, 2050075 (2020).
  * Zhang  _et al._ [2022] X.-H. Zhang, S.-D. Zhao, S. D. Mohanty, and Y.-X. Liu, Resolving Galactic binaries using a network of space-borne gravitational wave detectors, Phys. Rev. D 106, 102004 (2022), arXiv:2206.12083 [gr-qc] .
  * Gong  _et al._ [2021] Y. Gong, J. Luo, and B. Wang, Concepts and status of Chinese space gravitational wave detection projects, Nature Astronomy 5, 881 (2021), arXiv:2109.07442 [astro-ph.IM] .
  * Torres-Orjuela  _et al._ [2023] A. Torres-Orjuela, S.-J. Huang, Z.-C. Liang, S. Liu, H.-T. Wang, _et al._ , Detection of astrophysical gravitational wave sources by TianQin and LISA, arXiv e-prints , arXiv:2307.16628 (2023), arXiv:2307.16628 [gr-qc] .
  * Jin  _et al._ [2023] S.-J. Jin, Y.-Z. Zhang, J.-Y. Song, J.-F. Zhang, and X. Zhang, The Taiji-TianQin-LISA network: Precisely measuring the Hubble constant using both bright and dark sirens, arXiv e-prints , arXiv:2305.19714 (2023), arXiv:2305.19714 [astro-ph.CO] .
  * Amaro-Seoane  _et al._ [2023] P. Amaro-Seoane, J. Andrews, M. Arca Sedda, A. Askar, Q. Baghi, _et al._ , Astrophysics with the Laser Interferometer Space Antenna, Living Reviews in Relativity 26, 2 (2023), arXiv:2203.06016 [gr-qc] .
  * Arun  _et al._ [2022] K. G. Arun, E. Belgacem, R. Benkel, L. Bernard, E. Berti, _et al._ , New horizons for fundamental physics with LISA, Living Reviews in Relativity 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Sesana [2016] A. Sesana, Prospects for multiband gravitational-wave astronomy after gw150914, Physical Review Letters 116, 231102 (2016).
  * Kyutoku and Seto [2016] K. Kyutoku and N. Seto, Concise estimate of the expected number of detections for stellar-mass binary black holes by elisa, Monthly Notices of the Royal Astronomical Society 462, 2177 (2016).
  * Wong  _et al._ [2018] K. W. Wong, E. D. Kovetz, C. Cutler, and E. Berti, Expanding the lisa horizon from the ground, Physical review letters 121, 251102 (2018).
  * Moore  _et al._ [2019] C. J. Moore, D. Gerosa, and A. Klein, Are stellar-mass black-hole binaries too quiet for lisa?, Monthly Notices of the Royal Astronomical Society: Letters 488, L94 (2019).
  * Sato  _et al._ [2017] S. Sato, S. Kawamura, M. Ando, T. Nakamura, K. Tsubono, A. Araya, I. Funaki, K. Ioka, N. Kanda, S. Moriwaki, _et al._ , The status of DECIGO, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 840 (2017) p. 012010.
  * Kawamura  _et al._ [2021] S. Kawamura, M. Ando, N. Seto, S. Sato, M. Musha, I. Kawano, J. Yokoyama, T. Tanaka, K. Ioka, T. Akutsu, _et al._ , Current status of space gravitational wave antenna DECIGO and B-DECIGO, Progress of Theoretical and Experimental Physics 2021, 05A105 (2021), arXiv:2006.13545 [gr-qc] .
  * Liu  _et al._ [2020] S. Liu, Y.-M. Hu, J.-d. Zhang, J. Mei, _et al._ , Science with the tianqin observatory: Preliminary results on stellar-mass binary black holes, Physical Review D 101, 103027 (2020).
  * Nishizawa  _et al._ [2016] A. Nishizawa, E. Berti, A. Klein, and A. Sesana, elisa eccentricity measurements as tracers of binary black hole formation, Physical Review D 94, 064020 (2016).
  * Breivik  _et al._ [2016] K. Breivik, C. L. Rodriguez, S. L. Larson, V. Kalogera, and F. A. Rasio, Distinguishing between formation channels for binary black holes with lisa, The Astrophysical Journal Letters 830, L18 (2016).
  * Samsing and D’Orazio [2018] J. Samsing and D. J. D’Orazio, Black hole mergers from globular clusters observable by lisa i: eccentric sources originating from relativistic n-body dynamics, Monthly Notices of the Royal Astronomical Society 481, 5445 (2018).
  * Gerosa  _et al._ [2019] D. Gerosa, S. Ma, K. W. Wong, E. Berti, R. O’Shaughnessy, Y. Chen, and K. Belczynski, Multiband gravitational-wave event rates and stellar physics, Physical Review D 99, 103004 (2019).
  * Chamberlain and Yunes [2017] K. Chamberlain and N. Yunes, Theoretical physics implications of gravitational wave observation with future detectors, Physical Review D 96, 084039 (2017).
  * Gnocchi  _et al._ [2019] G. Gnocchi, A. Maselli, T. Abdelsalhin, N. Giacobbo, and M. Mapelli, Bounding alternative theories of gravity with multiband gw observations, Physical Review D 100, 064024 (2019).
  * Tso  _et al._ [2019] R. Tso, D. Gerosa, and Y. Chen, Optimizing ligo with lisa forewarnings to improve black-hole spectroscopy, Physical Review D 99, 124043 (2019).
  * Toubiana  _et al._ [2020a] A. Toubiana, S. Marsat, S. Babak, E. Barausse, and J. Baker, Tests of general relativity with stellar-mass black hole binaries observed by lisa, Physical Review D 101, 104038 (2020a).
  * McGee  _et al._ [2020] S. McGee, A. Sesana, and A. Vecchio, Linking gravitational waves and x-ray phenomena with joint lisa and athena observations, Nature Astronomy 4, 26 (2020).
  * Caputo  _et al._ [2020] A. Caputo, L. Sberna, A. Toubiana, S. Babak, E. Barausse, _et al._ , Gravitational-wave Detection and Parameter Estimation for Accreting Black-hole Binaries and Their Electromagnetic Counterpart, Astrophys. J. 892, 90 (2020), arXiv:2001.03620 [astro-ph.HE] .
  * Vitale [2016] S. Vitale, Multiband gravitational-wave astronomy: parameter estimation and tests of general relativity with space-and ground-based detectors, Physical Review Letters 117, 051102 (2016).
  * Cornish [2011] N. J. Cornish, Detection strategies for extreme mass ratio inspirals, Classical and Quantum Gravity 28, 094016 (2011).
  * Chua  _et al._ [2017] A. J. Chua, C. J. Moore, and J. R. Gair, Augmented kludge waveforms for detecting extreme-mass-ratio inspirals, Physical Review D 96, 044005 (2017).
  * Babak  _et al._ [2017] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, _et al._ , Science with the space-based interferometer lisa. v. extreme mass-ratio inspirals, Physical Review D 95, 103012 (2017).
  * Ewing  _et al._ [2021] B. Ewing, S. Sachdev, S. Borhanian, and B. Sathyaprakash, Archival searches for stellar-mass binary black holes in lisa data, Physical Review D 103, 023025 (2021).
  * Wang  _et al._ [2024] H. Wang, I. Harry, A. Nitz, and Y.-M. Hu, Space-based gravitational wave observatories will be able to use eccentricity to unveil stellar-mass binary black hole formation, Phys. Rev. D 109, 063029 (2024), arXiv:2304.10340 [astro-ph.HE] .
  * Buscicchio  _et al._ [2021] R. Buscicchio, A. Klein, E. Roebber, C. J. Moore, D. Gerosa, E. Finch, and A. Vecchio, Bayesian parameter estimation of stellar-mass black-hole binaries with lisa, Physical Review D 104, 044065 (2021).
  * Toubiana  _et al._ [2020b] A. Toubiana, S. Marsat, S. Babak, J. Baker, and T. Dal Canton, Parameter estimation of stellar-mass black hole binaries with lisa, Physical Review D 102, 124037 (2020b).
  * Lyu  _et al._ [2023] X. Lyu, E.-K. Li, and Y.-M. Hu, Parameter Estimation of Stellar Mass Binary Black Holes under the Network of TianQin and LISA, arXiv e-prints , arXiv:2307.12244 (2023), arXiv:2307.12244 [gr-qc] .
  * Kennedy and Eberhart [1995] J. Kennedy and R. Eberhart, Particle swarm optimization, in  _Proceedings of ICNN’95 - International Conference on Neural Networks_ , Vol. 4 (IEEE, 1995) pp. 1942–1948.
  * Shi [1998] Y. Shi, A modified particle swarm optimizer, in  _Proc of IEEE Icec Conference_ (1998).
  * Mohanty [2018] S. Mohanty,  _Swarm intelligence methods for statistical regression_ (CRC Press, 2018).
  * Weerathunga and Mohanty [2017] T. S. Weerathunga and S. D. Mohanty, Performance of particle swarm optimization on the fully-coherent all-sky search for gravitational waves from compact binary coalescences, Phys. Rev. D 95, 124030 (2017), arXiv:1703.01521 [gr-qc] .
  * Normandin and Mohanty [2020] M. E. Normandin and S. D. Mohanty, Towards a real-time fully-coherent all-sky search for gravitational waves from compact binary coalescences using particle swarm optimization, Phys. Rev. D 101, 082001 (2020), arXiv:2002.02150 [astro-ph.IM] .
  * Zhang  _et al._ [2021] X.-H. Zhang, S. D. Mohanty, X.-B. Zou, and Y.-X. Liu, Resolving Galactic binaries in LISA data using particle swarm optimization and cross-validation, Phys. Rev. D 104, 024023 (2021), arXiv:2103.09391 [gr-qc] .
  * Wang  _et al._ [2015] Y. Wang, S. D. Mohanty, and F. A. Jenet, Coherent Network Analysis for Continuous Gravitational Wave Signals in a Pulsar Timing Array: Pulsar Phases as Extrinsic Parameters, Astrophys. J. 815, 125 (2015), arXiv:1506.01526 [astro-ph.IM] .
  * Wang and Mohanty [2017] Y. Wang and S. D. Mohanty, Pulsar Timing Array Based Search for Supermassive Black Hole Binaries in the Square Kilometer Array Era, Phys. Rev. Lett. 118, 151104 (2017), arXiv:1611.09440 [astro-ph.IM] .
  * Riles [2023] K. Riles, Searches for continuous-wave gravitational radiation, Living Reviews in Relativity 26, 3 (2023), arXiv:2206.06447 [astro-ph.HE] .
  * Steltner  _et al._ [2021] B. Steltner, M. A. Papa, H.-B. Eggenstein, B. Allen, V. Dergachev, R. Prix, B. Machenschalk, S. Walsh, S. J. Zhu, O. Behnke, _et al._ , Einstein@ home all-sky search for continuous gravitational waves in ligo o2 public data, The Astrophysical Journal 909, 79 (2021).
  * Steltner  _et al._ [2023] B. Steltner, M. Papa, H.-B. Eggenstein, R. Prix, M. Bensch, B. Allen, and B. Machenschalk, Deep einstein@ home all-sky search for continuous gravitational waves in ligo o3 public data, The Astrophysical Journal 952, 55 (2023).
  * Bandopadhyay and Moore [2023] D. Bandopadhyay and C. J. Moore, LISA stellar-mass black hole searches with semicoherent and particle-swarm methods, Phys. Rev. D 108, 084014 (2023), arXiv:2305.18048 [gr-qc] .
  * Jaranowski  _et al._ [1998] P. Jaranowski, A. Krolak, and B. F. Schutz, Data analysis of gravitational-wave signals from spinning neutron stars: The signal and its detection, Physical Review D 58, 063001 (1998).
  * Crowder and Cornish [2007] J. Crowder and N. J. Cornish, Solution to the galactic foreground problem for lisa, Physical Review D 75, 043008 (2007).
  * Armstrong  _et al._ [1999] J. W. Armstrong, F. B. Estabrook, and M. Tinto, Time-Delay Interferometry for Space-based Gravitational Wave Searches, Astrophys. J. 527, 814 (1999).
  * Tinto and Dhurandhar [2005] M. Tinto and S. V. Dhurandhar, Time-Delay Interferometry, Living Reviews in Relativity 8, 4 (2005).
  * Prince  _et al._ [2002] T. A. Prince, M. Tinto, S. L. Larson, and J. W. Armstrong, Lisa optimal sensitivity, Phys. Rev. D 66, 122002 (2002).
  * Kay [1993] S. M. Kay,  _Fundamentals of statistical signal processing: estimation theory_ (Prentice-Hall, Inc., USA, 1993).
  * Finn [1992] L. S. Finn, Detection, measurement, and gravitational radiation, Phys. Rev. D 46, 5236 (1992), arXiv:gr-qc/9209010 [gr-qc] .
  * Brady and Creighton [2000] P. R. Brady and T. Creighton, Searching for periodic sources with LIGO. II. Hierarchical searches, Phys. Rev. D 61, 082001 (2000), arXiv:gr-qc/9812014 [gr-qc] .
  * Cutler  _et al._ [2005] C. Cutler, I. Gholami, and B. Krishnan, Improved stack-slide searches for gravitational-wave pulsars, Phys. Rev. D 72, 042004 (2005), arXiv:gr-qc/0505082 [gr-qc] .
  * Krishnan  _et al._ [2004] B. Krishnan, A. M. Sintes, M. A. Papa, B. F. Schutz, S. Frasca, and C. Palomba, Hough transform search for continuous gravitational waves, Phys. Rev. D 70, 082001 (2004), arXiv:gr-qc/0407001 [gr-qc] .
  * Antonucci  _et al._ [2008] F. Antonucci, P. Astone, S. D’Antonio, S. Frasca, and C. Palomba, Detection of periodic gravitational wave sources by hough transform in the f versus plane, Classical and Quantum Gravity 25, 184015 (2008).
  * Dhurandhar  _et al._ [2008] S. Dhurandhar, B. Krishnan, H. Mukhopadhyay, and J. T. Whelan, Cross-correlation search for periodic gravitational waves, Phys. Rev. D 77, 082001 (2008), arXiv:0712.1578 [gr-qc] .
  * Bratton and Kennedy [2007] D. Bratton and J. Kennedy, Defining a standard for particle swarm optimization, in  _2007 IEEE Swarm Intelligence Symposium_ (2007) pp. 120–127.
  * Engelbercht [2005] A. Engelbercht, Fundamentals of computational swarm intelligence, University of Pretoria, South Africa (2005).
  * Marsat and Baker [2018] S. Marsat and J. G. Baker, Fourier-domain modulations and delays of gravitational-wave signals, arXiv e-prints , arXiv:1806.10734 (2018), arXiv:1806.10734 [gr-qc] .
  * Cornish and Rubbo [2003] N. J. Cornish and L. J. Rubbo, LISA response function, Phys. Rev. D 67, 022001 (2003).
  * Zhou  _et al._ [2021] M.-Y. Zhou, X.-C. Hu, B. Ye, S. Hu, D.-D. Zhu, X. Zhang, W. Su, and Y. Wang, Orbital effects on time delay interferometry for TianQin, Phys. Rev. D 103, 103026 (2021), arXiv:2102.10291 [astro-ph.IM] .
  * Vallisneri [2005] M. Vallisneri, Synthetic LISA: Simulating time delay interferometry in a model LISA, Phys. Rev. D 71, 022001 (2005), arXiv:gr-qc/0407102 [gr-qc] .
  * Rubbo  _et al._ [2004] L. J. Rubbo, N. J. Cornish, and O. Poujade, Forward modeling of space-borne gravitational wave detectors, Phys. Rev. D 69, 082003 (2004).
  * Estabrook and Wahlquist [1975] F. B. Estabrook and H. D. Wahlquist, Response of Doppler spacecraft tracking to gravitational radiation., General Relativity and Gravitation 6, 439 (1975).
  * Feng  _et al._ [2019] W.-F. Feng, H.-T. Wang, X.-C. Hu, Y.-M. Hu, and Y. Wang, Preliminary study on parameter estimation accuracy of supermassive black hole binary inspirals for TianQin, Phys. Rev. D 99, 123002 (2019), arXiv:1901.02159 [astro-ph.IM] .
  * Cutler [1998] C. Cutler, Angular resolution of the lisa gravitational wave detector, Phys. Rev. D 57, 7089 (1998).
  * Królak  _et al._ [1995] A. Królak, K. D. Kokkotas, and G. Schäfer, Estimation of the post-newtonian parameters in the gravitational-wave emission of a coalescing binary, Phys. Rev. D 52, 2089 (1995).
  * Buonanno  _et al._ [2009] A. Buonanno, B. R. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, Comparison of post-newtonian templates for compact binary inspiral signals in gravitational-wave detectors, Phys. Rev. D 80, 084043 (2009).
  * Crowder and Cornish [2007] J. Crowder and N. J. Cornish, Solution to the galactic foreground problem for LISA, Phys. Rev. D 75, 043008 (2007), arXiv:astro-ph/0611546 [astro-ph] .