Resolving Galactic binaries using a network of space-borne gravitational wave detectors

Gravitational wave (GW) astronomy is now well-established in the frequency band $\approx[10,1000]$ Hz using ground-based interferometric detectors. Starting with the discovery in 2015 of GW150914 [1], a binary black hole (BBH) merger signal, by the twin LIGO [2] detectors, the number of confirmed GW signals has grown to $90$ over three observing runs [3, 4, 5] of the LIGO and Virgo [6] detectors. While the majority of signals are from BBH mergers, there is also a double-neutron star [7] and two neutron star-black hole signals [8] in the haul so far.

The success of GW astronomy with ground-based detectors has further spurred the drive to open up the gravitational wave window at lower frequencies. Pulsar Timing Arrays (PTAs) are already being used for the $\approx[10^{-9},10^{-7}]$ Hz and a tantalizing signal, not yet confirmed to be from GWs, has emerged [9, 10, 11] in the most recent datasets collected by the NANOGrav [12], PPTA [13] and EPTA [14] collaborations. The $\approx[10^{-4},1]$ Hz band is the target for the space-borne LISA (Laser Interferometer Space Antenna) [15] detector, scheduled for launch in 2037. Targeting a similar frequency band, intensive work is in progress on the design of the Taiji [16] and Tianqin [17] missions, both of which have expected launch dates in the middle of the next decade.

While both LISA and Taiji will have heliocentric orbits, Tianqin is planned to be geocentric. In other respects, however, these detectors share the common feature of having three spacecrafts in a triangular formation that sense GW-induced distance distortion using inter-spacecraft laser links. The nominal distance between the spacecrafts will be $\approx 2.5\times 10^{6}$ km for LISA and $\approx 3\times 10^{6}$ km for Taiji, while Tianqin will have shorter arm lengths of $\approx\sqrt{3}\times 10^{5}\ \mathrm{km}$. The normals to the planes of both the LISA and Taiji constellations will have an approximate tilt of $60^{\circ}$ from the normal to the ecliptic and rotate once around the latter in a year. For Tianqin, the normal to the spacecraft plane will point towards a fixed sky location, namely that of the double white dwarf RX J0806.3+1527, at all times.

All planned space-based GW detectors will employ the technique of time-delay interferometry (TDI) [18], in which the laser frequency Doppler readouts from different arms are combined linearly after introducing time delays, to substantially reduce laser frequency noise. In this paper we consider the so-called $A$, $E$, and $T$ TDI combinations that have mutually independent noise.

The operational frequency bands of LISA, Taiji, and Tianqin are expected to be rich in astrophysical sources [19, 20]. These include a large population of compact object Galactic binaries (GBs) [21], comprised mostly of white dwarfs and reaching $O(10^{7})$ in number, tens to hundreds [22] of Extreme Mass Ratio Inspirals (EMRIs) – a system containing a massive black hole in the $10^{4}$ to $10^{7}\ M_{\odot}$ range orbited by a stellar-mass range compact object [23], and a handful of Massive black hole binaries (MBHB) containing comparable $\approx 10^{4}$ to $10^{6}\ M_{\odot}$ components [24]. Due to the generically broad antenna patterns of interferometric detectors and longer lifetimes of GW sources at lower frequencies, the data from these detectors will be simultaneously occupied by signals from a large number and variety of sources. The multisource resolution problem of disentangling these signals from each other presents a major data analysis challenge that has motivated a number of different approaches and algorithms.

To provide a benchmark for the development and comparison of data analysis methods, the LISA community has organized several mock LISA data challenges (MLDCs) [25, 26, 27, 28]. The most recent in this series are the LISA data challenges (LDCs) [29]. The first LDC (LDC1 or Radler), which is the challenge considered in this paper, consists of several subchallenges among which subchallenge 4 (LDC1-4) focuses purely on GB multisource resolution. LDC1-4 provides single realizations of TDI combinations, each containing Gaussian stationary instrumental noise added to the GW signals from $30\times 10^{6}$ GBs. A catalog of the parameters of the GBs is also provided, allowing the performance of a multisource resolution method to be quantified rigorously.

Investigations carried out with the mock data challenges [30, 28] indicate that, over an observation period of $\simeq 2$ yr with a single detector, it would be possible to confidently resolve $O(10^{4})$ individual GBs in the $\approx[0.1,15]$ mHz band [31, 32, 33, 34]. The remaining GBs will blend together to form a stochastic signal – the GB background – that will likely dominate over instrumental noise below $\approx 3$ mHz. Since the GB background will be the main factor limiting the sensitivity of searches for all other GW sources in this band, it is important to develop methods that can extract the maximum number of resolvable GBs without grossly overfitting the data.

Given that, until recently, LISA was the only anticipated mission, much of the literature on data analysis for space-based detectors has focused on a single detector. Some notable exceptions are the works that have analyzed a mission concept called the Big Bang Observer (BBO) [35], proposed as a successor to LISA. BBO is envisioned to contain $4$ LISA-like detectors but with shorter arm lengths of $\approx 5\times 10^{4}$ km. Three of the detectors will be spaced $120^{\circ}$ apart in heliocentric orbits while the fourth will be coincident in location with one of them. In addition to multiple detectors, enhancements such as higher laser power are assumed in order to reach the overall sensitivity expected to reveal the inflationary GW stochastic background.

While BBO remains a conceptual mission, the scenario of multiple space-based detectors does not appear to be a pipe dream any longer given the advent of Taiji and Tianqin, which will likely overlap in operation with LISA. As such, the question of what new science capabilities will be unlocked by the existence of a network of space-based detectors has become a timely one. Several recent works have approached this question using the computationally inexpensive Fisher information formalism [36, 37, 38, 39, 40, 41, 42] in the context of single sources. It has been shown that considerable benefits can be derived by combining the data from multiple space-based detectors. However, the analysis of single sources using the Fisher information approach, which only provides an asymptotic lower bound on parameter estimation errors, does not provide a realistic estimate of performance for the full multisource resolution problem. For the latter, there is no substitute for the analysis of realistic mock data with an actual, not ideal, data analysis pipeline.

In this paper, we report an analysis of mock data from a network of two space-based detectors and provide realistic estimates of the detection and estimation performance on the GB multisource resolution problem. The mock data is generated using the same set of GBs as in LDC1-4 under the assumption, which is adequate for a first study and simplifies mock data generation, that Taiji is identical in design to LISA, and that one detector trails the Earth while the other leads it by $20^{\circ}$. The mock data is analyzed using the latest version of the GB resolution pipeline, GBSIEVER (Galactic Binary Separation by Iterative Extraction and Validation using Extended Range), introduced in [33], that has been generalized to handle TDI data from multiple detectors.

As in the earlier paper [33], referred to as P1 from here on, it is convenient to define the following types of GB sources when discussing our results. (i) True: sources in the LDC1-4 catalog. (ii) Reported: the final list of estimated sources returned by GBSIEVER. (iii) Identified: The initial list of sources produced by the single source search step in GBSIEVER before various cuts are imposed to construct the set of reported sources. (iv) Confirmed: The reported sources that match true sources as determined by a prescribed metric for association. The fraction of confirmed sources in the set of reported ones is called the detection rate.

Our key result is that, at comparable detection rates, a LISA-Taiji network is capable of resolving $\approx 75\%$ more confirmed binaries than either of the detectors alone. In addition, we find that the residual left after subtracting the reported signals from the data tends to have fewer cases of missed strong sources compared to the single detector case. The power spectral density (PSD) of the residual reaches the theoretical lower bound on confusion noise, obtained under the assumption of a perfect but infeasible multisource resolution method for a single detector [43], without significant overfitting. Thus, not only can a detector network probe deeper into the GB population for resolvable sources but the lowered GB background will enhance the detectability of non-GB sources. While parameter estimation errors are reduced with the LISA-Taiji network for sources that are common with the LISA-only analysis, the scaling of error variances with signal strength that is expected from a Fisher information based analysis is not observed for all the parameters.

The rest of the paper is organized as follows. The data used in this paper are described in Sec. II. The baseline single-source detection and parameter estimation method for a detector network is described in Sec. III. A self-contained but brief overview of GBSIEVER is described in Sec. IV. This is followed by the results of our analysis in Sec. V. Our conclusions and prospects for future work are presented in Sec. VI.

Space-based GW detectors will use TDI to strongly suppress laser frequency noise. TDI is implemented by linearly combining readouts of frequency shifts between the incoming and outgoing light along each arm after introducing known time delays [18]. There are a variety of combinations and levels of approximations, called TDI generations, that have been proposed. We use the so-called $A$ and $E$ combinations, which have mutually independent instrumental noise, corresponding to the first generation of TDIs used in LDC1 [44]. (The $T$ combination is dropped because the GW signal in it is highly attenuated.)

To obtain the TDI GW signals, we start with the two polarizations in the transverse traceless (TT) gauge of a plane GW incident at the Solar System barycenter (SSB) origin. For a GB, the polarizations are well-modeled in the SSB frame as linear chirps,

where $\mathcal{A}$ is the amplitude of the wave, $\iota$ is the inclination angle between the GB orbital angular momentum and the line of sight from the SSB origin to the GB, $\phi_{0}$ is the initial phase, $f$ is the frequency at the start of observations, and $\dot{f}$ is the secular frequency drift. We use the formalism in [31] to generate the TDI GW signals for GB sources.

The TDI time series $\overline{y}^{I}_{D}$ for combination $I$ and detector $D$ is given by

where $\overline{x}\in\mathbb{R}^{N}$ denotes a row vector, $\overline{s}^{I}_{D}(\theta)$ denotes a single GB signal corresponding to source parameters $\theta$, $\overline{n}^{I}_{D}$ is a realization of the instrumental noise, and $N_{s}$ is the number of GBs. For a GB, $\theta$ consists of $\{\mathcal{A},\phi_{0},\iota,\psi,\lambda,\beta,f,\dot{f}\}$, where $\mathcal{A}$, $\phi_{0}$, $\iota$, $f$, and $\dot{f}$ were defined following Eq. 1 to Eq. 3), $\psi$ is the polarization angle defining the orientation of the binary orbit projected on the sky, and the longitude and the latitude of the source in the SSB frame are denoted by $\lambda$ and $\beta$, respectively. In the rest of the paper, $\theta$ will serve as a stand in for a GB source itself where convenient. Following LDC1-4 we set the sampling frequency for the uniformly sampled time series $\overline{y}^{I}_{D}$, for all $I$ and $D$, to be $f_{s}=1/15$ Hz, and the number of samples to $N=4194304$ corresponding to an observation period $T_{\rm{obs}}\approx 2$ yr.

For the generation of mock TDI data, the LISA community has adopted a standard code called LISACode [45] for generating instrumental noise and a code called FastGB for GB signals. (The latter is based on the formalism in [46] and is included in the software distribution associated with LDC data [47].) The simulated LDC1-4 data contains the TDI signals from $30$ million GBs, with GW signal frequencies in $\approx[0.1,15]$ mHz. In the astrophysically realistic model used for the GB population [48, 49], the density of GBs drops with increasing frequency with most of the sources above $4$ mHz being individually resolvable. The instrumental noise realization is drawn from a Gaussian, stationary, stochastic process. The PSD of the noise is denoted by $S_{D}^{I}(\nu)$ at Fourier frequency $\nu$. In LISACode, the PSD of the instrumental noise in the A and E combinations are identical and derived from the design sensitivity of LISA.

We generate simulated Taiji TDI data by shifting the centroid of the orbit $40^{\circ}$ ahead of LISA and using the same software suite and GB catalog as LDC1-4. Statistical independence of the instrumental noise in the two detectors is ensured by using different seeds for the corresponding pseudo-random sequences in LISACode. We keep the arm length of Taiji the same as that of LISA. Besides allowing LISACode to be used with minimal changes at the code level, this leads to a simpler and more direct comparison, as befits a first study, of the performance gained in switching from a single detector to a network. To acknowledge this simplification, the simulated Taiji data is called Taiji-mod in the remainder of the paper. Extending our analysis to a realistic Taiji would be quite straightforward when its own data simulation framework is made available in the public domain.

The GW signal in each TDI combination has amplitude and frequency modulations arising from the time-varying antenna response due to the rotation of the spacecraft constellation and the Doppler effect due to its orbital motion, respectively. Consequently, the spectrum of the TDI GB signal is broadened over a frequency range that is $\approx 10^{4}\times f$ times larger than the separation of Fourier frequencies for a $\approx 2$ yr observation period. The broadening increases the overlap of spectra from multiple GBs, with the number of overlaps increasing at lower frequencies and signal strengths. This precludes the differentiation of the signals using a simple Fourier transform and makes multisource resolution an extremely challenging data analysis problem.

GBSIEVER implements an iterative scheme in which the parameters of only a single source are estimated at a time and the corresponding signal is subtracted from the data. The parameter estimation step uses maximum likelihood estimation (MLE) where the log-likelihood function is constructed under the assumption that the data contains a single source added to Gaussian, stationary noise. The maximization of the log-likelihood in the case of a GB can be carried out analytically over a subset of the parameters, leaving behind a function, commonly called the $\mathcal{F}$-statistic in the GW literature [50], that needs to be numerically maximized over the remaining ones. In GBSIEVER, the latter maximization is carried out using particle swarm optimization (PSO) [51].

The principal difference between P1 and the present paper is the shift in the mathematical formalism implemented in GBSIEVER, which is a straightforward generalization from a single detector to, as described in this section, a network of detectors.

The principal algorithmic components of GBSIEVER and their corresponding settings are reviewed in this section. The description is self-contained but brief since further details are available in P1. We do not provide an extensive review of PSO here, which is used for the global optimization task in Eq. 12, since it is a widely known algorithm. A pedagogical introduction to PSO and its application to MLE can be found in [52]. The specific variant of PSO used in GBSIEVER is discussed in P1 and its parameter settings are described in detail in Sec. III of [53].

The presentation of results in this section is organized as follows. Sec. V.1 contains results on source resolution. An investigation of the residual is presented in Sec. V.2. Sec. V.3 contains results on GB parameter estimation.

We investigated the GB multisource resolution problem in the context of a network of comparable space-based GW detectors, namely, LISA and Taiji, using a previously developed iterative source subtraction pipeline, GBSIEVER, that was extended to analyze TDI combinations from multiple detectors. Our results provide a realistic assessment of the benefits expected from a network of space-based detectors that goes beyond Fisher information studies of parameter estimation errors for single GBs.

In the present work, we looked at detection and estimation performance of the LISA-Taiji network relative to LISA alone. There is a significant increase, by $\approx 75\%$, in the number of confirmed sources that can be extracted while maintaining approximately the same detection rate. This increase was obtained without changing any of the cuts (SNR and $R_{\rm ee}$) imposed on the identified sources or the cross-validation setup for primary and secondary runs.

One of the key results we have reported is the comparison of the residual after subtraction of reported sources found in the network analysis with the ideal GB background of a single detector. We find that a network of space-based detectors will allow the GB background to be reduced substantially below the lowest level possible with a single detector. The PSD of the residual lies entirely below the ideal GB background obtained with a detection threshold of ${\rm SNR}_{\rm bg}^{\rm min}=7.0$. It coincides with the ${\rm SNR}_{\rm bg}^{\rm min}=5.0$ background over much of the frequency range in which the latter dominates over instrumental noise. Since the residual from subtracting GBs forms the noise background for non-GB searches, one expects that their performance will be significantly improved in the $\approx[0.5,4]$ mHz band. Further studies need to be carried out to quantify this expectation.

Our study of parameter estimation errors brings forth two salient points that highlight the difference between a single-source and multisource problem. (i) There is no significant change in the distribution of the differences in parameter values of confirmed and true sources between a single-detector and a network of comparable ones. (ii) Fisher information based error analysis of single sources may not be adequate in the multi-source resolution problem since the expected scaling of errors with SNR is observed only for the subset of intrinsic parameters but not the extrinsic ones. While there is a clear explanation for (i), we do not understand the origin of (ii) yet.

We have assumed simplified models for both the LISA and Taiji detectors, treating them as identical in design except for an angular separation in their heliocentric orbit. In future work, we will use more realistic models (the orbits and individual TDI combinations) for LISA and Taiji data to investigate the performance improvements for LISA-Taiji network. We will also take into account the uncertainty in prior knowledge of the PSDs of instrumental or GB background noise. This can be achieved, in principle, by a modest expansion of the search space for PSO by including the ratios of the noise PSDs in the same frequency band as free parameters to be optimized along with those of a single GB. We have further assumed that the LISA and Taiji data are collected over the same observation period. Since the GBs are persistent sources, and since the $\mathcal{F}$-statistic treats data from different detectors and TDI combinations additively, much of our analysis would go through, barring code modifications, even if the observation periods were non-overlapping. The inclusion of strongly evolving GBs could be an interesting exception here, if the observation periods are far apart, that will require a more careful approach.

Zhang, Zhao and Liu are supported by the National Key Research and Development Program of China grants No. 2021YFC2203000 and No. 2020YFC2201400, the National Natural Science Foundation of China through Grant No. 12047501, the 111 Project under grant No. B20063, and Lanzhou City’s scientific research funding subsidy to Lanzhou University. Group activities are financially supported by the Morningside Center of Mathematics and a part of the MPG-CAS collaboration in low frequency gravitational wave physics. Partial support from the Xiandao B project in gravitational wave detection is also acknowledged. We are grateful to Prof. Shing-Tung Yau for his long term and unconditional support and Prof. Yun-Kau Lau for assembling our research group. We gratefully acknowledge the use of high performance computers at the Supercomputing Center of Lanzhou University and State Key Laboratory of Scientific and Engineering Computing, CAS.