Follow-up procedure for gravitational wave searches from isolated neutron stars using the time-domain $\mathcal{F}$-statistic method

TL;DR

This paper introduces a follow-up procedure using the time-domain $ ext{F}$-statistic for detecting gravitational waves from rotating neutron stars, improving analysis in all-sky searches with LIGO and Virgo data.

Contribution

The paper presents a new follow-up method based on matched-filtering $ ext{F}$-statistic, integrated into the existing pipeline for gravitational wave searches from neutron stars.

Findings

The follow-up procedure shows good agreement with theoretical predictions.

Tests with artificial signals demonstrate the method's effectiveness.

The method will be incorporated into routine all-sky search pipelines.

Abstract

Among promising sources of gravitational waves are long-lived nearly periodic signals produced by rotating, asymmetric neutron stars. Depending on the astrophysical scenario, the sources of asymmetry may have thermal, viscous, elastic and/or magnetic origin. In this work we introduce a follow-up procedure for an all-sky search for gravitational wave signals from rotating neutron stars. The procedure denoted as Followup implements matched-filtering $\mathcal{F}$-statistic method. We describe data analysis methods and algorithms used in the procedure. We present tests of the Followup for artificial signals added to white, Gaussian noise. The tests show a good agreement with the theoretical predictions. The Followup will become part of the Time-Domain $\mathcal{F}$-statistic pipeline that is routinely used for all-sky searches of LIGO and Virgo detector data.

Equations

$$h_0=4\times 10^{-25}\left(\frac{ϵ}{10^{-6}}\right)\left(\frac{I_3}{10^{45}{\text{ g cm}}^2}\right){\left(\frac{f}{100\text{ Hz}}\right)}^2\left(\frac{100\text{ pc}}{d}\right),$$

$$s(t)=\sum _{i=1}^4{a}_i{h}_i(t),$$

$$\begin{array}[]{c}h_{1}(t)=a(t)\cos\phi(t),\quad h_{2}(t)=b(t)\cos\phi(t),\\[8.61108pt] h_{3}(t)=a(t)\sin\phi(t),\quad h_{4}(t)=b(t)\sin\phi(t),\end{array}$$

$$\varphi (t)=\omega t+\dot{\omega }{t}^2+\frac{\mathbf{n}\cdot {\mathbf{r}}_d(t)}{c}(\omega +2\dot{\omega }t),$$

$$⟨g(t)⟩=\frac{1}{{\sigma }^2}\sum _{t=1}^{N}g(t),$$

$$\mathcal{F}=\frac{|F_a{|}^2}{⟨a^2(t)⟩}+\frac{|F_b{|}^2}{⟨b^2(t)⟩},$$

$$F_a:=⟨x(t)a(t)\exp [-\mathrm{i}\varphi (t)]⟩,$$

$$F_b:=⟨x(t)b(t)\exp [-\mathrm{i}\varphi (t)]⟩.$$

$${\rho }^2=⟨s^2(t)⟩.$$

$${\Gamma }_{{\theta }_i{\theta }_j}=⟨\frac{\partial s(t)}{\partial {\theta }_i}\frac{\partial s(t)}{\partial {\theta }_j}⟩,i,j=1,\dots ,M.$$

$$\mu$$

$${\Sigma }^2$$

$${\rho }_c=\sqrt{2({\mathcal{F}}_c-2)}.$$

$$\mathsf{a}:=\left(\begin{array}{c}{a}_1\\ ⋮\\ a_n\end{array}\right),\mathsf{h}(t;\mathbf{\xi }):=\left(\begin{array}{c}{h}_1(t;\mathbf{\xi })\\ ⋮\\ h_n(t;\mathbf{\xi })\end{array}\right).$$

$$s(t;\mathbf{\theta })={\mathsf{a}}^{\mathsf{T}}\cdot \mathsf{h}(t;\mathbf{\xi }),$$

$$\Gamma (\mathsf{a},\mathbf{\xi })=\left(\begin{array}{cc}{\Gamma }_{\mathsf{aa}}(\mathbf{\xi })& {\Gamma }_{\mathsf{a}\mathbf{\xi }}(\mathsf{a},\mathbf{\xi })\\ {\Gamma }_{\mathsf{a}\mathbf{\xi }}(\mathsf{a},\mathbf{\xi }{)}^{\mathsf{T}}& {\Gamma }_{\mathbf{\xi }\mathbf{\xi }}(\mathsf{a},\mathbf{\xi })\end{array}\right),$$

$${\Gamma }_{\mathsf{aa}}(\mathbf{\xi })$$

$${\Gamma }_{\mathsf{a}\mathbf{\xi }}(\mathsf{a},\mathbf{\xi })$$

$${\Gamma }_{\mathbf{\xi }\mathbf{\xi }}(\mathsf{a},\mathbf{\xi })$$

$$M_{ij}(\mathbf{\xi }):=⟨h_i(t;\mathbf{\xi })h_j(t;\mathbf{\xi })⟩,i,j=1,\dots ,n.$$

$$F_{(\mathcal{A}){ij}}(\boldsymbol{\xi}):=\bigg{<}h_{i}(t;\boldsymbol{\xi})\,\frac{\partial h_{j}(t;\boldsymbol{\xi})}{\partial\xi_{\mathcal{A}}}\bigg{>},\quad\mathcal{A}=1,\ldots,m,\quad i,j=1,\ldots,n,$$

$$S_{(\mathcal{A}\mathcal{B}){ij}}(\boldsymbol{\xi}):=\bigg{<}\frac{\partial h_{i}(t;\boldsymbol{\xi})}{\partial\xi_{\mathcal{A}}}\,\frac{\partial h_{j}(t;\boldsymbol{\xi})}{\partial\xi_{\mathcal{B}}}\bigg{>},\quad\mathcal{A},\mathcal{B}=1,\ldots,m,\quad i,j=1,\ldots,n,$$

$$\widetilde{\Gamma}_{\mathcal{A}\mathcal{B}}(\boldsymbol{\xi}):=\frac{1}{n}\text{Tr}\Big{(}\mathsf{M}(\boldsymbol{\xi})^{-1}\cdot\mathsf{A}_{(\mathcal{A}\mathcal{B})}(\boldsymbol{\xi})\Big{)},\quad\mathcal{A},\mathcal{B}=1,\ldots,m,$$

$${\mathsf{A}}_{(\mathcal{A}\mathcal{B})}(\mathbf{\xi }):=S_{(\mathcal{A}\mathcal{B})}(\mathbf{\xi })-{\mathsf{F}}_{(\mathcal{A})}(\mathbf{\xi }{)}^{\mathsf{T}}\cdot \mathsf{M}(\mathbf{\xi }{)}^{-1}\cdot {\mathsf{F}}_{(\mathcal{B})}(\mathbf{\xi }),\mathcal{A},\mathcal{B}=1,\dots ,m.$$

$$M_g=M_o\cdot V^{\prime }\cdot \sqrt{E},$$

$$s(t)=h_o\cos (\omega t+\dot{\omega }{t}^2+{\varphi }_o)$$

$$x(t)=n(t)+s(t),$$

$$\mathcal{F}=\frac{|{\sum }_{t=1}^{N}x(t)\exp [i(\omega t+\dot{\omega }{t}^2)]{|}^2}{N{\sigma }^2},$$

$$\rho =\frac{h_o}{\sigma }\sqrt{\frac{N}{2}}$$

$$\Gamma=\rho^{2}\left(\begin{array}[]{cc}\frac{1}{12}N^{2}&\frac{1}{12}N^{3}\\[4.30554pt] \frac{1}{12}N^{3}&\frac{4}{45}N^{4}\end{array}\right).$$

$${\sigma }_{\omega }$$

$${\sigma }_{\dot{\omega }}$$

$$M_{o}=R\sqrt{3}\left(\begin{array}[]{cc}1&0\\ \frac{1}{2}&\frac{\sqrt{3}}{2}\end{array}\right),$$

Full text

††: CQG

1 Introduction

Historically the first gravitational-wave (GW) signal registered on Earth (Abbott et al., 2016a) was caused by a merger of a stellar-mass black hole binary system. Since then nine more such events were observed (Abbott et al., 2016b, 2017a, 2017b, 2017c, 2018b) by LIGO and Virgo laser interferometric detectors. Additionally, one binary neutron star (NS) merger (Abbott et al., 2017d) was detected.

Due to the unremitting efforts of the LIGO and Virgo (Acernese et al., 2014) Collaboration (LVC) in the detectors upgrades and in improvements of the data analysis methods, it will be possible to detect much subtler signals, in particular those emitted by rotating asymmetric neutron stars. Several searches for continuous gravitational waves (CGW) from isolated neutron stars have been carried out in LIGO and Virgo data (see Riles 2017 for a recent review). These searches have included coherent searches for gravitational radiation from known radio and X-ray pulsars, directed searches for known stars or locations having unknown signal frequencies, and spotlight or all-sky searches for stars with unknown signal frequency and sky location. Even though no statistically-significant signals were detected, interesting upper limits on the CGW emission were placed. CGW are expected to have strain amplitudes a few orders of magnitude weaker than signals produce by mergers. However in contrast to the merger signals, CGW are long-lasting, and therefore one can improve the sensitivity of the searches by increasing the observational time.

According to the current theoretical state-of-art, several mechanisms may be responsible for a long-lasting GW. In the case of young and isolated neutrons stars (NS), strong, evolving toroidal magnetic field and unstable oscillation modes (r-modes, Owen et al. 1998) may lead to detectable signals. NS in accreting binaries may accumulate accreted material on their surface. Deformation due to mass, temperature gradients, elastic or magnetic field strain result in a non-axisymmetric shape of NS and CGW emission (Bonazzola & Gourgoulhon, 1996; Bildsten, 1998; Ushomirsky et al., 2000). In general, any rotating non-axisymmetric NS (exhibiting time-varying mass quadrupole moment) will produce GW; for the review see Andersson et al. (2011); Lasky (2015); Riles (2017). Canonical emitters of CGW and main targets for the LVC searches are rapidly-rotating, non-axisymmetric NS in our Galaxy. Their model is motivated by the relation between GW strain amplitude $h_{0}$ and the spin frequency $f$ (Zimmerman & Szedenits, 1979):

[TABLE]

where $d$ is the distance to the source, $\epsilon=(I_{1}-I_{2})/I_{3}$ and $I_{1}$, $I_{2}$, $I_{3}$ denote moments of inertia along three axes, with the direction of $I_{3}$ aligned with the axis of the NS angular momentum.

One of the pipelines to perform all-sky search for gravitational waves from isolated rotating neutron stars is Time-Domain ${\mathcal{F}}$-statistic pipeline111Project’s repository: https://github.com/mbejger/polgraw-allsky. This pipeline implements a matched-filter statistic for detecting nearly periodic signals in the time-domain data, called the ${\mathcal{F}}$-statistic (Jaranowski et al., 1998; Jaranowski & Królak, 2009, 2010; Astone et al., 2010; Pisarski & Jaranowski, 2015). The Time-Domain ${\mathcal{F}}$-statistic pipeline is divided into several steps: generation of time domain data, construction of grid of templates, ${\mathcal{F}}$-statistic search for candidate signals, search for coincidences between candidates in different time segments, calculation of false alarm probability for coincidences obtained, estimation of sensitivity of the search if no significant signal is detected.

There are several other all-sky search pipelines that are used in the searches for gravitational waves from isolated neutron stars in LIGO and Virgo detectors data: Einstein@Home (Abbott et al., 2009) 222http://einstein.phys.uwm.edu, FrequencyHough (Astone et al., 2014), PowerFlux (Abbott et al., 2008), and SkyHough (Krishnan et al., 2004) (for a recent review see Bejger 2017).

In this paper we present a follow-up procedure which aims at verifying whether promising candidates obtained by the Time-Domain ${\mathcal{F}}$-statistic pipeline can be of astrophysical origin. The follow-up procedure is based on the assumption that gravitational wave signal is a coherent signal that is always present in the data. Consequently when coherence time for matched-filtering analysis increases the signal-to-noise ratio increases as square root of the observational time.

The article is composed as follows: in Sect. 2 we describe the model of the expected CGW signal and the ${\mathcal{F}}$-statistic method. In Sect. 3 we introduce the Followup procedure: description of steps (Sect. 3.1), construction of the optimal grid of templates (Sect. 3.2) and optimisation algorithms used to find maximum of the ${\mathcal{F}}$-statistic over the parameter space (Sect. 3.3). In Sect. 4 we present implementation of the Followup procedure for the of CGW signal buried in white Gaussian noise. We perform Monte Carlo simulations for two cases: the two-dimensional parameter space (Sect. 4.1), as well for general, four-dimensional parameter space (Sect. 4.2). Sect. 5, summarizes and discusses application of the Followup procedure to the Time-Domain ${\mathcal{F}}$-statistic pipeline.

2 CGW signal model and the Time-Domain ${\mathcal{F}}$-statistic method

The Followup algorithm is a part of the Time-Domain ${\mathcal{F}}$-statistic all-sky search pipeline. The pipeline performs a coherent search for CGW signals in the time-domain data segments using the ${\mathcal{F}}$-statistic (Jaranowski et al., 1998). Here, for the context, we provide a brief introduction to this method.

The time-domain response $s(t)$ of the interferometric detector to a CGW signal is given by a linear combination of four time-dependent components:

[TABLE]

where the functions $h_{i}$ ($i=1,\dots,4$) are of the form

[TABLE]

The functions $a(t)$ and $b(t)$ are the amplitude modulation functions that depend on the location and orientation of the detector on the Earth and on the position of the GW source in the sky, described in the equatorial coordinate system by the right ascension $\alpha$ and the declination $\delta$ angles. They are periodic functions of time with the period of one and two sidereal days. The analytic form of the functions $a(t)$ and $b(t)$ for the case of interferometric detectors is given by Eqs. 12, 13 of Jaranowski et al. (1998). The phase $\phi(t)$ is given by

[TABLE]

where ${\bf r}_{d}(t)$ is the vector that joins the solar-system barycenter (SSB) with the detector, and ${\bf n}$ is the unit vector pointing from the SSB to the source. In equatorial coordinates $(\delta$, $\alpha)$ we have ${\bf n}=(\cos\delta\cos\alpha,\cos\delta\sin\alpha,\sin\delta)$. In the following we assume that the frequency evolution is accurately described by one spindown parameter $\dot{\omega}$. The four amplitudes $a_{i}$ depend on amplitude $h_{0}$, phase $\phi_{0}$, polarization angle $\psi$, and inclination angle $\iota$ (see Eqs. 2.10 of Astone et al. 2010).

CGW signals are identified with the use of the ${\mathcal{F}}$-statistic introduced in Jaranowski et al. (1998). The ${\mathcal{F}}$-statistic is obtained by maximizing the likelihood function with respect to the four unknown parameters - $h_{0}$, $\phi_{0}$, $\iota$, and $\psi$ (which are henceforth called the extrinsic parameters). This leaves a function of only - in the most basic case - four remaining parameters: $\omega$, $\dot{\omega}$, $\delta$, and $\alpha$ (henceforth called the intrinsic parameters). The angular frequency $\omega$ and angular frequency derivative $\dot{\omega}$ are dimensionless and they are defined by $\omega=2\pi fdt$ and $\dot{\omega}=2\pi\dot{f}dt^{2}$, where $dt$ is sampling time. Thus the dimension of the parameter space that we need to search decreases from 8 to 4. Additionally, we introduce several simplifications. We set the observation time $T_{0}$ equal to the integer multiple of sidereal days. We also assume that the noise in the detector is white. This is a good approximation because we can assume that over a very narrow band of the signal spectral density of the noise of the detector is constant. We assume that data $x(t)$ is a discrete time series consisting of $N$ uniformly sampled data points, i.e., $t=1,...,N$. We also introduce a product $\left\langle\cdot\right\rangle$ defined as

[TABLE]

where $\sigma^{2}$ is the variance of the noise.

Under these assumptions the ${\mathcal{F}}$-statistic is given by

[TABLE]

where

[TABLE]

The signal-to-noise ratio $\rho$ is defined as

[TABLE]

For Gaussian noise $\rho$ determines probability of the detection of signal $s(t)$. Assuming Gaussian noise the Fisher matrix for signal $s(t)$ is given by

[TABLE]

where $\theta=(\theta_{1},\ldots,\theta_{M})$ are $M$ parameters of the signal $s(t)$. For sufficiently high signal-to-noise ratio the accuracy of the determination of parameters is approximately given by the covariance matrix which is equal to the inverse of the Fisher matrix. As we shall see in the following the Fisher matrix is also useful in the construction of the grid of templates. Often it is convenient to consider twice the value of the ${\mathcal{F}}$-statistic. This is because $2\times{\mathcal{F}}$ has a central $\chi^{2}$ distribution with 4 degrees of freedom when the signal is absent and non-central $\chi^{2}$ distribution with 4 degrees of freedom and non-centrality parameter $\rho^{2}$ when the signal is present (see Section IIIB of Jaranowski et al. 1998). Thus when the signal is present mean ($\mu$) and variance ($\Sigma^{2}$) of $2\times{\mathcal{F}}$ read

[TABLE]

with $n=4$.

A fully-coherent search over the whole set of data (typically of the length of a year) is computationally prohibitive. In order to alleviate this problem we apply a semi-coherent method, which consists of dividing the data into shorter time domain segments. The short time domain data are analysed coherently with the ${\mathcal{F}}$-statistic. Moreover, to reduce the computer memory required to do the search, the data are divided into narrow-band segments that are analysed separately. For a typical search we choose the length of the time-domain segment to be several sidereal days long and the bandwidth of the narrow band segment to be a fraction of $1$ Hz. Consequently the Time-Domain ${\mathcal{F}}$-statistic all-sky search pipeline consists of two parts. The first part is a coherent search of narrow-band time domain data segments, where we search a 4-parameter space defined by angular frequency $\omega$, angular frequency derivative $\dot{\omega}$, declination $\delta$, and right ascension $\alpha$. The search is performed on an effective 4-dimensional grid in the parameter space described in detail in Astone et al. (2010); Pisarski & Jaranowski (2015). We set a fixed threshold for the ${\mathcal{F}}$-statistic for each data segment. All the threshold crossings are recorded together with corresponding 4 parameters of the grid point and the signal-to-noise ratio $\rho_{c}$ at the threshold crossing. The signal-to-noise ratio $\rho_{c}$ that we record is defined in terms of the ${\mathcal{F}}$-statistic value ${\mathcal{F}}_{c}$ at the threshold crossing as

[TABLE]

The formula follows from Eq. 10 where for mean value $\mu$ we substitute $2\times{\mathcal{F}}_{c}$.

In this way for each time domain segment we obtain a set of candidates. The second part of the analysis consists of the search for coincidences among the candidates from different time segments. The coincidence procedure is described in detail in Section 8 of Aasi et al. (2014). We estimate the statistical significance of the coincidences by calculating the probability that a given coincidence is by chance only. The formula for this false alarm probability is given in the Appendix of Aasi et al. (2014). Whenever false alarm probability is sufficiently small we mark the coincidence as significant. Typically we choose the false alarm probability equal to $10^{-3}$. When we have a candidate coincident in $l$ time domain segments we estimate the parameters of the coincident candidate as the mean of the parameters of the $l$ candidates from individual time frames entering the coincidence. The significant candidates are then subject to the follow-up procedure which is presented in this paper.

So far, the pipeline has been used in several all-sky searches for CGW signals: Virgo VSR1 data (Aasi et al., 2014), LIGO O1 data (Abbott et al., 2017e, 2018a), LIGO O2 data (Abbott et al., 2019) and in the LIGO S6 Mock Data Challenge (Walsh, 2016).

In the case of an all-sky search, when large parameter space $(\omega,\dot{\omega},\delta,\alpha)$ has to be investigated scrupulously, optimal usage of the computational resources and reduction of the computational costs play an important role. On the one hand computational cost of the all-sky search depends on the observation time $T_{0}$ as ${\propto}$ $T_{0}^{5}\log T_{0}$ (Astone et al., 2010). On the other hand the signal-to-noise ratio scales as $\rho\,{\propto}\sqrt{T_{0}}$. These two scaling relations show how important balance between sensitivity and computational power is: in the too short data segments signal will be buried deeply in the noise and impossible to restore, while too long time series require unreachable computational power. Interplay between reduction of numerical cost and sensitivity loss indicate that search settings, like e.g. length of the data segments, sampling time, density of the search grid etc., have to be chosen carefully.

3 The Followup procedure

The aim of the Followup procedure is to verify whether significant coincident candidates are of astrophysical interest. We assume that the signal’s lifetime is longer than the duration of the observing run and consequently we assume that a signal of astrophysical origin will always be present in the data analysed. From the two-step search procedure of the Time-Domain ${\mathcal{F}}$-statistic all-sky pipeline we have candidate signals present in L time domain segments, each of D sidereal days long. As a coherent signal of true astrophysical origin should be present in all the data i.e. in L $\times$ D sidereal days long, the signal-to-noise ratio should increase by factor $\sqrt{L}$ with respect to the signal-to-noise in one segment. In addition the accuracy of the estimation of parameters of the signal should also increase considerably. The exact increase in accuracy is complicated by the fact that parameters of the signal are correlated. To facilitate the task of following up the signal we implement a hierarchical procedure. Namely we first search for a coherent signal (evaluate the ${\mathcal{F}}$-statistic) in two concatenated time-domain segments, then in four segments and so on.

3.1 Steps of the Followup procedure

The candidate signal summarized in Sect. 2 is described by four parameters: $\omega_{i}$, $\dot{\omega}_{i}$, $\delta_{i}$, $\alpha_{i}$.

The steps of the Followup procedure are the following:

Take two adjacent time domain segments of the band where the candidate signal is present. 2. 2.

For each segment construct an optimal grid (see Sect. 3.2) around the parameters of the candidate signal. 3. 3.

Perform a coarse search by evaluating the ${\mathcal{F}}$-statistic on the grid, find the maximum value of ${\mathcal{F}}$, and the corresponding coarse estimates of signal parameters: $\omega_{c}$, $\dot{\omega}_{c}$, $\delta_{c}$, $\alpha_{c}$. 4. 4.

Using an optimization algorithm (see Sect. 3.3) perform a fine search for the maximum of the ${\mathcal{F}}$-statistic and find the corresponding refined estimates of signal parameters. Take an average of the refined signal parameters from the two segments to obtain the fine estimates $\omega_{f}$, $\dot{\omega}_{f}$, $\delta_{f}$, $\alpha_{f}$ of four signal parameters. 5. 5.

Join the two segments together to obtain one segment of double length, and construct the optimal grid around the fine values of the parameters $\omega_{f}$, $\dot{\omega}_{f}$, $\delta_{f}$, $\alpha_{f}$ obtained in the previous step. Perform again a coarse and a fine search. This will result in the final Followup parameters $\omega_{fin}$, $\dot{\omega}_{fin}$, $\delta_{fin}$, $\alpha_{fin}$.

At the end of the procedure, the SNR of a signal of astrophysical origin should increase approximately by a factor of $\sqrt{2}$ with respect to its SNR in individual segments. Also the estimates $\omega_{fin}$, $\dot{\omega}_{fin}$, $\delta_{fin}$, $\alpha_{fin}$ should be more accurate estimates of the parameters of the signal than the initial estimates $\omega_{i}$, $\dot{\omega}_{i}$, $\delta_{i}$, $\alpha_{i}$. We can then iterate the above procedure to further join the segments. For a true signal the SNR should continue to increase with the above-mentioned factor.

3.2 Optimal grid

The grid used in the coherent ${\mathcal{F}}$-statistic search was introduced in Astone et al. (2010) and was further optimized in Pisarski & Jaranowski (2015). This grid is optimal for an approximate linear model of the CGW signal (see Sections IIIB and IV of Astone et al. 2010). The resulting grid is uniform, which simplifies the search procedure considerably. Moreover, the grid is constrained so that the $\omega$ grid points coincide with the Fourier frequencies. This enables the use of the FFT algorithm and results in computational speed-up.

In the Followup procedure we are less computationally-bound and we can use a fully optimal grid based on the reduced Fisher matrix introduced in Jaranowski & Królak (2005). Let us summarize the steps of the construction of the reduced Fisher matrix (for more details see Section 4.4.1 of Jaranowski & Królak 2005 and Chapter 6 of Jaranowski & Królak 2009). We present the construction for the general case of $n$ amplitude parameters and $m$ intrinsic parameters. First we introduce the following shorthand notation. Let us collect the $n$ amplitude parameters $a_{i}$ and the $n$ waveforms $h_{i}$ into column vectors.

[TABLE]

With this notation the signal $s(t)$ can compactly be written in the following form:

[TABLE]

where $\boldsymbol{\theta}$ is the set of signal parameters, $\mathsf{T}$ stands for the matrix transposition and $\cdot$ denotes matrix multiplication. The the $n+m$ signal parameters $\boldsymbol{\theta}$ consist of $n$ extrinsic amplitude parameters $\mathsf{a}$ and $m$ intrinsic parameters $\boldsymbol{\xi}$. For the case of the gravitational wave signal from a rotating neutron star presented in Section 2 there are $n=4$ amplitude parameters and 4 intrinsic parameters, $\boldsymbol{\xi}=(\omega,\dot{\omega},\delta,\alpha)$.

The Fisher matrix $\Gamma$ for the general signal $s(t;\boldsymbol{\theta})$ can be written in terms of block matrices as

[TABLE]

where $\Gamma_{\mathsf{a}\mathsf{a}}$ is an $n\times n$ matrix with components $\left\langle{\partial s}/{\partial a_{i}}\,{\partial s}/{\partial a_{j}}\right\rangle$ $(i,j=1,\ldots,n)$, $\Gamma_{\mathsf{a}\boldsymbol{\xi}}$ is an $n\times m$ matrix with components $\left\langle{\partial s}/{\partial a_{i}}\,{\partial s}/{\partial\xi_{\mathcal{A}}}\right\rangle$ $(i=1,\ldots,n,\,\mathcal{A}=1,\ldots,m)$, and $\Gamma_{\boldsymbol{\xi}\boldsymbol{\xi}}$ is an $m\times m$ matrix with components $\left\langle{\partial s}/{\partial\xi_{\mathcal{A}}}\,{\partial s}/{\partial\xi_{\mathcal{B}}}\right\rangle$ $(\mathcal{A},\mathcal{B}=1,\ldots,m)$. The explicit form of these matrices is

[TABLE]

The components of the matrix $M$ are given by

[TABLE]

The components of $m$ matrices $\mathsf{F}_{(\mathcal{A})}$ (note here the index $\mathcal{A}$ within parentheses has a meaning of the matrix label), and the components of $m^{2}$ matrices $\mathsf{S}_{(\mathcal{AB})}$ read

[TABLE]

[TABLE]

respectively. We then define a $m\times m$ square matrix $\widetilde{\Gamma}$ with components

[TABLE]

where $\mathsf{A}_{(\mathcal{A}\mathcal{B})}$ is an $n\times n$ matrix defined as

[TABLE]

We shall call $\widetilde{\Gamma}$ the reduced Fisher matrix. This matrix is a function of the intrinsic parameters alone and is used to construct a grid on the intrinsic parameter space.

The construction of a grid in the parameter space is equivalent to a mathematical problem known as the covering problem (Conway & Sloane, 1999). In the Euclidean space the covering problem is to cover the space by hyperspheres of radius $R$ in such a way that any point of the space belongs to at least one hypersphere. The radius $R$ is called the covering radius. The covering thickness $\Theta$ is defined as the average number of hyperspheres that contain a point in the space. The optimal covering has the minimal possible thickness. Thus the optimal covering is also referred to as the thinnest (see Section IV of Astone et al. 2010 for more details).

Let $M_{o}$ be the generator matrix of the thinnest covering in the Euclidean space. To obtain the generator matrix $M_{g}$ in the intrinsic parameter space we need to transform $M_{g}$ to space spanned by intrinsic parameters $\boldsymbol{\xi}$. This is obtained in the following way

[TABLE]

where $V$ is the matrix whose columns are eigenvectors of the reduced Fisher matrix $\widetilde{\Gamma}$, $V^{\prime}$ denotes transpose of matrix $V$ and $E$ is the diagonal matrix whose diagonal components are eigenvalues of $\widetilde{\Gamma}$ (see Section IV of Astone et al. 2010 for more details).

The constructed grid depends on the values of the intrinsic parameters, i.e. it changes throughout the parameter space. However, in a small neighbourhood around the parameters of the candidate signal the changes are small, and the grid constructed at the candidate’s location is used.

3.3 Optimisation procedure

We have implemented three algorithms to perform a fine search for the maximum of the ${\mathcal{F}}$-statistic in the Followup procedure described above: (i) Simplex method (Spendley et al., 1962; Nelder & Mead, 1965) - non-derivative optimisation method based on the simplex (triangle in arbitrary number of dimensions), which relocates and adapts to the local features of the parameter space, and contracts itself toward extremum direction. (ii) Mesh Adaptive Direct Search (MADS), introduced by Audet & Dennis (2006) - a grid-based algorithm, which explores parameter space, changing its grid-mesh size, depending on the distance from the extremum. Application of this algorithm to the CGW searches can be found in Shaltev & Prix (2013). (iii) Inverted Mesh Adaptive Direct Search (invMADS) algorithm - our modification of algorithm (ii) which makes optimisation procedure more suitable for the four-dimensional ${\mathcal{F}}$-statistic case.

A visualisation of the simplex algorithm can be found in on Figure 1 of Królak (1999). Schematic visualisation of the invMADS algorithm is shown on Fig. 1. The procedure starts from the point with the highest value of ${\mathcal{F}}$-statistic evaluated on the optimal grid. This initial seed point is denoted as a star marker on the plot. Around this point a four-dimensional (two-dimensional on the plot) hypercube is constructed, and the ${\mathcal{F}}$-statistic is evaluated on the edges, vertices, faces and cells of the hypercube. If the initial seed has a higher value of the ${\mathcal{F}}$-statistic than on the hypercube points, the mesh is expanded around the initial seed, and the evaluation procedure is repeated. If algorithm finds a higher value of the ${\mathcal{F}}$-statistic, the seed is changed to that point and the procedure is repeated, starting from the initial size of the hypercube.

Original MADS method as introduced by Audet & Dennis (2006) is based on dividing the parameter space into relatively large parts and on gradual shrinking toward the direction of the extremum. For the ${\mathcal{F}}$-statistic where by our follow-up procedure the initial seed point is relatively close to the extremum we have modified the MADS algorithm inverting its idea. The new procedure starts from a small mesh which progressively expands. The pseudo-code of the inverted MADS procedure is as follows:

Comparison of the three algorithms is shown on Fig. 2. The only criterion of choosing parameters fixed for each procedure, like for example Initial Mesh Size or Increase parameters for invMADS (see pseudo-code above), was the good convergence of the results to the theoretical predictions (Cramer-Rao bounds, see Sect. 4.2 and Fig. 9) for a single data segment. We find essential differences in the algorithms performance for one data segment (left column) and two data segments joined together (right column). For the standard and the inverted MADS the execution time increases several times when switching from the shorter to the longer time series. This is because for the longer time series the $\mathcal{F}$-statistic becomes more narrowly peaked, as was shown on Fig. 4 - narrower extrema are harder to find in the 4-dimensional space and it is more time-consuming to check whether the extremum is the global one. Execution time for Simplex algorithm is similar for single and concatenated data segments. This can be explained by the bottom panels on Fig. 2. The algorithm finds correct extrema (SNR loss is close to 0) for a single data segment however for the longer time series, where, as was mentioned above, it is harder to find global extremum, algorithm is as fast as previously but it gets stuck in some local peaks (SNR loss for Simplex in the Follow-up stage is very large ($\sim 0.8$) ). Standard MADS procedure does not fit our purposes: it generates relatively large SNR loss for a single data segment as well as in the Follow-up stage. The most promising algorithm in our test is inverted MADS. SNR loss for our inverted MADS both for shorter and concatenated data segments remains stable and is close to 0. Disadvantage of that method is that it is slow. However this is not a problem as the accuracy and not the computational speed is the main requirement for the follow-up procedure. Additionally, invMADS can be easily sped-up, e.g. by adding parallelisation to the procedure whereas Simplex is very hard to parallelise.

4 Tests in Gaussian noise

In this chapter we present implementation of the Followup procedure in computer codes and we carry out Monte Carlo simulations to present efficiency of our follow-up method. We present two simulations for the case of CWG signals buried in white Gaussian noise. First simulation is for a signal that depends only on 4 parameters: amplitude $h_{o}$, phase $\phi_{o}$, angular frequency $\omega$, and frequency derivative $\dot{\omega}$. The ${\mathcal{F}}$-statistic for such a signal depends only on two parameters - $\omega$ and $\dot{\omega}$. The second simulation is for the general signal presented in Section 2 that depends on eight parameters. The main principle of our follow-up procedure is to coherently join the data segments into a longer time series to increase detection sensitivity.

In Fig. 3 we plot sample data and the ${\mathcal{F}}$-statistic used in our follow-up procedure. The data consists of the signal given by Eq. 2 with the optimal signal-to-noise ratio $\rho=20$ (Eq.8) added to white Gaussian noise (upper panel). We see the signal is deeply buried in the noise. The lower panel shows the ${\mathcal{F}}$-statistic (Eq.6). We see that ${\mathcal{F}}$-statistic has several subsidiary maxima. The maxima are due to the amplitude modulation and are separated by multiples of sidereal day frequency. The amplitude modulation functions have 2 harmonics of sidereal frequency. The ${\mathcal{F}}$-statistic involves modulus squared of the amplitude demodulated data and thus it has in general 4 harmonics on both sides of the main peak. The relative amplitudes of the subsidiary maxima depend on the declination of the source and for the declination of the signal used in Fig. 3 only 5 subsidiary maxima are visible.

In Fig. 4 we show the ${\mathcal{F}}$-statistic as a function of frequency for one segment and two segments joined together. We see that for two segments joined together the amplitude of the ${\mathcal{F}}$-statistic increases and width of the function decreases; the ${\mathcal{F}}$-statistic becomes more narrowly peaked.

4.1 Two-dimensional case

In this case the signal $s(t)$ is given by

[TABLE]

We assume that the signal present in the data $x(t)$ is additive i.e.

[TABLE]

where $n(t)$ is Gaussian white noise of mean [math] and variance $\sigma^{2}$. The ${\mathcal{F}}$-statistic in this case is given by

[TABLE]

where $x_{k}(t),t=1,...,N$ are $N$ data points.

For signal (22), the signal-to-noise ratio $\rho$ is approximately given by

[TABLE]

and the Fisher matrix projected on the two-dimensional space spanned by angular frequency and the spin down parameter is given by

[TABLE]

Consequently the variances of the estimators of $\omega$ and $\dot{\omega}$ approximately read

[TABLE]

Next we present construction of the optimal grid outlined in Section 3.2 for the case of the signal given by Eq. 22. Here the intrinsic parameter space is two-dimensional and spanned by angular frequency $\omega$ and the spin down $\dot{\omega}$. The thinnest covering in the two-dimensional space is the hexagonal lattice denoted by $A_{2}^{*}$ (Conway & Sloane, 1999).

The generator matrix of this lattice is given by

[TABLE]

where $R$ is the covering radius of the lattice. The reduced Fisher matrix $\tilde{\Gamma}$ for signal (Eq. 22) reads

[TABLE]

The generator matrix $M_{2}$ transformed to the parameters $\omega_{0}$ and $\omega_{1}$ is obtained from the equation

[TABLE]

where $V$ is the matrix whose columns are eigenvectors of the reduced Fisher matrix $\tilde{\Gamma}$, $V^{\prime}$ denotes transpose of matrix $V$ and $E$ is the diagonal matrix whose diagonal components are eigenvalues of $\tilde{\Gamma}$. To characterize the density of the grid we introduce parameter $MM$ defined as

[TABLE]

The denser the grid the closer the value of the parameter $MM$ to $1$.

For our simulations we generate a set of data containing $2N$ data points. The data consists of Gaussian noise of mean 0 and variance 1 and signal (Eq. 22) added to it. We first divide the data into two segments of $N$ data points and perform the search and estimation of the parameters in the two segments for various signal-to-noise ratios according to the first step of the Followup procedure presented in Sect. 3.1. To obtain the signal with a specified SNR we scale the amplitude $h_{o}$ appropriately. We first verify whether the detectability of the signal agrees with the theoretical one. In the case of signal present $2\times{\mathcal{F}}$-statistic defined by Eq. 24 has a non-central $\chi^{2}$ distribution of 2 degrees of freedom with non-centrality parameter equal to $\rho^{2}$. Consequently the theoretical mean ($\mu$) and variance ($\Sigma^{2}$) for $2\times{\mathcal{F}}$ are given by Eqs.10 with $n=2$.

In Figure 5 we plot the mean and variance of ${\mathcal{F}}$-statistic obtained from our simulations as functions of the SNR and compare them with the theoretical ones. For each SNR we perform 250 simulations. The good agreement for the mean and variance after the Followup means that the Followup procedure achieves the theoretical increase in the signal-to-noise ratio and leads to best possible increase in signal detectability.

The results of the estimation of parameters are presented in Fig. 6. We calculate means and standard deviations for 250 estimators of parameters $\omega$ and $\dot{\omega}$ that we obtain and we compare them with the true values of these parameters and the standard deviations predicted by the Fisher matrix and given by Eqs. 27 and 28. The results of the comparison are presented on Fig. 6 and Fig. 7. The bias for a given parameter $\theta$ expressed in % is defined as

[TABLE]

where $\bar{\theta}$ is the mean value of the estimators of the parameter $\theta$ and $\theta$ is its true value.

We note a very good agreement of results our simulations with the Fisher matrix above $\rho=7$. We then carry out the second stage of the Followup procedure i.e. the search for the signal in the whole $2N$ data set and we again see a good agreement with the Fisher matrix predictions above $\rho=12$ (see Fig. 7). In both stages of the procedure we used the optimal grid constructed above with parameter $MM=0.98$.

4.2 General four-dimensional case

Our code allows us to inject artificial signal with an arbitrary amplitude or equivalently with an arbitrary SNR. To check performance of our follow-up code we injected signals with different values of SNR ranging from 8 to 20. Signal was buried in white Gaussian noise. To make our simulation more realistic and applicable to real search where we analyse many consecutive time segments we generated data consisting of six data segments. In each data segment we search for the signal using our two-step procedure: first a grid search using the optimal grid constructed in Sect. 3.2 and then a fine search using the Simplex algorithm. We than calculate the mean values of the parameters of the signal obtained from the fine search weighted by the estimated SNRs of the signal in each frame. These mean values constitute the initial values of the parameters for the follow-up procedure. We then apply the follow-up procedure to the first two data segments that we join together using the invMADS algorithm (Sect. 3.3) with the initial values of the parameters. We repeat whole procedures 250 times, for different realisations of the Gaussian noise. The results of the simulation are show in Fig. 8 and Fig. 9. In Fig. 8 we show the bias of the estimators of the four parameters and in Fig. 9 the standard deviation in comparison with the predictions of the Fisher matrix. For the follow-up case we have plotted the bias and the standard deviation as a function of signal-to-noise ratio for the two segments joined together which is a factor of $\sqrt{2}$ larger than for one segment. We note a satisfactory performance of the follow-up procedure. Namely the bias of the estimators decreases as we move from one segment to two segments joined together and standard deviations of the estimators improve in accordance with the Fisher matrix predictions.

Similarly, like in the two-dimensional case (Sect. 4.1), we also verify whether the detectability of the signal agrees with the theoretical one. In the 4-dimensional case $2\times{\mathcal{F}}$-statistic has a non-central $\chi^{2}$ distribution of 4 degrees of freedom with the non-centrality parameter equal to $\rho^{2}$. Consequently the theoretical mean ($\mu$) and variance ($\Sigma^{2}$) for $2\times{\mathcal{F}}$ are given by Eq. 10 with $n=4$. In Fig. 10 we present the comparison of our simulations with theoretical expectations for an array of signal-to-noise ratios.

We have also investigated how on average the SNR of the injected signal increases in the course of our follow-up procedure. The increase is measured by the quantity $\sqrt{\frac{\langle\rho_{2}^{2}\rangle}{\langle\rho_{1}^{2}\rangle}}$ where $\langle\rho_{1}^{2}\rangle$ and $\langle\rho_{2}^{2}\rangle$ are the averaged squared SNRs of the signal in two segments joined together and one segment respectively. The obtained results are shown in Fig. 11. We see a good agreement with the expected increase of SNR by $\sqrt{2}$.

5 Conclusions

In this paper we have presented a follow-up procedure for the analysis of the candidate signals obtained from the Time-Domain ${\mathcal{F}}$-statistic pipeline. The aim of the procedure is to verify whether a given candidate signal can be of astrophysical signal. The procedure is based on the assumption that a true signal is coherent and always present in the data and thus signal-to-noise ratio of the signal should increase as square root of the observation time.

The basic tool of the procedure is the ${\mathcal{F}}$-statistic. We have described data analysis methods and algorithms used in the procedure. They involve construction of optimal grids of templates over which the ${\mathcal{F}}$-statistic is evaluated and optimization algorithms to find accurately the global maximum of the ${\mathcal{F}}$-statistic. We have presented detailed steps of the procedure. The follow-up would proceed hierarchically - we first join two data segments and analyse them coherently, then concatenate data segments into four data segments and so on. In practice the candidate events from our pipeline are results of search for coincidences in $L$ data segments. A significant candidate need not be coincident in all $L$ data segments. Nevertheless by our basic assumption a candidate of astrophysical origin must be always present in the data. Thus in our follow-up we proceed with the steps outlined in Section 3.1 to search for signal in two consecutive segments as if the signal was present in both segments.

We have performed Monte Carlo simulations of the Followup by injecting artificial signal into white Gaussian noise. We have tested the procedure in two cases. Firstly, for a simple two-dimensional model when only two parameters - frequency $\omega$ and spindown $\dot{\omega}$ are unknown. Then we have performed tests for a general four-dimensional signal that is modulated by the motion of the detector and depends also on the position of the source in the sky. Our test involves comparison of the estimated parameters of the signal with predictions of the Fisher matrix. We find a satisfactory agreement. Thus we expect that the procedure will become an integral part of the Time-Domain ${\mathcal{F}}$-statistic pipeline and will be applied to the analysis of real data. As the methods presented here are quite general, they can be applied to any other continuous GW signal, e.g. in the searches for QCD axion clouds bound to black holes (Arvanitaki et al., 2015).

Acknowledgements

The work was supported in part by the Polish National Science Centre grants no. 2016/22/E/ST9/00037, 2017/26/M/ST9/00978 and 2018/28/T/ST9/00458.