Use of Excess Power Method and Convolutional Neural Network in All-Sky Search for Continuous Gravitational Waves

As stated in Sec. III, the time resampling technique can demodulate the phase, but complete demodulation is not available because of the limitation of computational resources. In our work, the time resampling technique is employed to eliminate only the effect caused by the Earth’s diurnal rotation, $\Phi_{\oplus}(t)$. Assuming a representative grid point $\bm{{n}}_{\mathrm{g}}$, we can rewrite the phase (6) as

	$\displaystyle\Phi(t)$	$\displaystyle=2\pi{f}_{\mathrm{gw}}t+\Phi_{\oplus}(t)+\Phi_{\odot}(t)$
		$\displaystyle=2\pi{f}_{\mathrm{gw}}\zeta+\delta\Phi_{\oplus}(t)+\delta\Phi_{\odot}(t)\,,$		(22)

where

	$\displaystyle\delta\Phi_{\oplus}(t)$	$\displaystyle:=2\pi f_{\mathrm{gw}}\frac{\bm{{r}}_{\oplus}(t)\cdot\Delta\bm{{n}}}{c}\,,$		(23)
	$\displaystyle\delta\Phi_{\odot}(t)$	$\displaystyle:=2\pi f_{\mathrm{gw}}\frac{\bm{{r}}_{\odot}(t)\cdot\Delta\bm{{n}}}{c}\,,$		(24)

and $\Delta\bm{{n}}:=\bm{{n}}_{\mathrm{s}}-\bm{{n}}_{\mathrm{g}}$. Since the residual phase varies with time, we will place grid points so that the amplitude of the residual phase in the worst case, i.e.,

$$\min_{\bm{{n}}_{\mathrm{g}}}\max_{t}|\delta\Phi_{\oplus}(t)|$$

becomes smaller than a threshold $\delta\Phi_{\epsilon}$ for any source direction $\bm{{n}}_{\mathrm{s}}$ within the area covered by the grid point ${\bm{{n}}}_{\mathrm{g}}$. To optimize the grid placement, we employ the method proposed in Ref. Nakano et al. (2003). The residual phase $\delta\Phi_{\oplus}$ is expanded up to the first order of $\Delta\alpha:=\alpha_{\mathrm{s}}-\alpha_{\mathrm{g}}$ and $\Delta\delta:=\delta_{\mathrm{s}}-\delta_{\mathrm{g}}$. Then, we get

	$\displaystyle\delta\Phi_{\oplus}\simeq$	$\displaystyle\frac{2\pi{f}_{\mathrm{gw}}}{c}R_{\text{E}}\cos\lambda\left\{-\Delta\delta\sin\delta_{\mathrm{g}}\cos(\alpha_{\mathrm{g}}-\varphi_{\oplus}-\Omega_{\oplus}t)\right.$
		$\displaystyle\left.-\Delta\alpha\cos\delta_{\mathrm{g}}\sin(\alpha_{\mathrm{g}}-\varphi_{\oplus}-\Omega_{\oplus}t)\right\}\,.$		(25)

Here, the constant term is neglected because it degenerates with the initial phase $\phi_{0}$. The maximum value of the residual phase is

	$\displaystyle\max_{t}|\delta\Phi_{\oplus}|=$	$\displaystyle\frac{2\pi{f}_{\mathrm{gw}}}{c}R_{\text{E}}|\cos\lambda|$
		$\displaystyle\times\sqrt{(\Delta\delta)^{2}\sin^{2}{\delta}_{\mathrm{g}}+(\Delta\alpha)^{2}\cos^{2}{\delta}_{\mathrm{g}}}\,.$		(26)

The grid points are to be determined to satisfy $\max_{t}|\delta\Phi_{\oplus}|<\delta\Phi_{\epsilon}$ for any source direction.

Because the residual phase (26) is symmetric under the transformation ${\delta}_{\mathrm{g}}\to-{\delta}_{\mathrm{g}}$, the placement of grids on the negative $\delta$ side can be generated by inverting the sign of the grids on the positive $\delta$ side. Therefore, we focus on the case with $0\leq\delta\leq\pi/2$.

Since the residual phase depends only on $\delta$ at $\delta_{\mathrm{g}}=\pi/2$, a single template can cover the neighbor of $\delta=\pi/2$. In fact, at $\delta=\pi/2$, Eq. (26) becomes

$$\max_{t}|\delta\Phi_{\oplus}|=\frac{2\pi{f}_{\mathrm{gw}}}{c}R_{\text{E}}|\Delta\delta|\cos\lambda\,.$$

(27)

Therefore, the condition $\max_{t}|\delta\Phi_{\oplus}|\leq\delta\Phi_{\epsilon}$ gives the lower bound of $\delta_{1}$ such that the region $\delta_{1}\leq\delta\leq\pi/2$ can be covered by a signle patch represented by $\{({\alpha}_{\mathrm{g}},{\delta}_{\mathrm{g}})=(0,\pi/2)\}$, to find

$$\delta_{1}:=\frac{\pi}{2}-\delta\Phi_{\epsilon}\times\frac{c}{2\pi{f}_{\mathrm{gw}}}\frac{1}{R_{\text{E}}\cos\lambda}\,.$$

(28)

Plural patches are necessary to cover the strip of a constant $\delta$ in the other range. We introduce a 2-dimensional metric corresponding to the residual phase (26),

$$d\sigma^{2}=\cos^{2}\delta d\alpha^{2}+\sin^{2}\delta d\delta^{2}\,.$$

(29)

In general, a metric in a 2-dimensional manifold can be transformed into a conformally flat metric by an appropriate coordinate transformation. When the space is conformally flat, the curve of a small constant distance measured from an arbitrary chosen point can be approximated by a circle. Therefore, a template spacing in the 2-dimensional parameter space becomes relatively easy. By defining new variables $X:=\alpha$ and $Y:=-\log|\cos\delta|$, the metric can be transformed into

$$d\sigma^{2}=e^{-2Y}(dX^{2}+dY^{2})\,.$$

(30)

Along with Nakano et al. (2003), we can construct the sky patches covering the half-sky region with $0\leq\delta\leq\delta_{1}$. Figure 1 shows a part of grid points constructed under the condition

$$\delta\Phi_{\epsilon}=0.058\,,$$

(31)

which we adopt throughout this paper. The total number of grid points to cover the whole sky is

$$N_{\rm grid}=352,436\,,$$

(32)

for $f_{\mathrm{gw}}$ = 100Hz.

As we choose $\delta\Phi_{\epsilon}$ to be sufficiently small, we neglect $\delta\Phi_{\oplus}$ in the following discussion. Then, after subtracting the phase modulation due to the Earth’s rotation, the phase of the gravitational wave (22) becomes

$$\Phi(t)=2\pi{f}_{\mathrm{gw}}\zeta+\delta\Phi_{\odot}(t).$$

(33)

We apply the short-time Fourier transform (STFT) to the time-resampled strain,

$$s(\zeta)={h}_{\mathrm{obs}}(\zeta)+n(\zeta).$$

(34)

In the rest of the paper, we treat only the time-resampled data. Therefore, without confusion, the time-resampled data in Eq. (34) can be denoted by the same character as the original one. The strain is divided into ${N}_{\mathrm{seg}}$ segments having the duration ${T}_{\mathrm{seg}}$ and their start times are denoted by $\zeta_{j}:=jT_{\mathrm{slide}},\ (j=0,1,\cdots,{N}_{\mathrm{seg}}-1)$. $T_{\mathrm{slide}}$ is not necessary to be equal to $T_{\mathrm{seg}}$. The output of STFT with the window function $w(\zeta)$ is defined by

$$s^{\mathrm{STFT}}_{j,k}=h^{\mathrm{STFT}}_{j,k}+n^{\mathrm{STFT}}_{j,k}\,,$$

(35)

where

$$h^{\mathrm{STFT}}_{j,k}=\frac{1}{T_{\mathrm{seg}}}\int^{\zeta_{j}+T_{\mathrm{seg}}}_{\zeta_{j}}d\zeta^{\prime}\ w(\zeta^{\prime}-\zeta_{j}){h}_{\mathrm{obs}}(\zeta^{\prime})e^{-2\pi if_{k}\zeta^{\prime}}\,,$$

(36)

$$n^{\mathrm{STFT}}_{j,k}=\frac{1}{T_{\mathrm{seg}}}\int^{\zeta_{j}+T_{\mathrm{seg}}}_{\zeta_{j}}d\zeta^{\prime}\ w(\zeta^{\prime}-\zeta_{j})n(\zeta^{\prime})e^{-2\pi if_{k}\zeta^{\prime}}\,,$$

(37)

and $f_{k}:=k\Delta f=k/T_{\mathrm{seg}}$ is the frequency of the $k$-th element of STFT. Let us focus on the positive frequency modes, i.e., $f_{k}>0$. Then, the second term of Eq. (13) can be neglected and Eq. (36) can be approximated by

	$\displaystyle{h}^{\mathrm{STFT}}_{j,k}\simeq$	$\displaystyle\frac{1}{{T}_{\mathrm{seg}}}\int^{\zeta_{j}+{T}_{\mathrm{seg}}}_{\zeta_{j}}d\zeta^{\prime}\ \bigg{\{}w(\zeta^{\prime}-\zeta_{j})$
		$\displaystyle\times G(t(\zeta^{\prime}))e^{2\pi i\delta f_{k}\zeta^{\prime}}e^{i\delta\Phi_{\odot}(\zeta^{\prime})}\bigg{\}}.$		(38)

with $\delta f_{k}:={f}_{\mathrm{gw}}-f_{k}$. In the expression of $G(t(\zeta))$, the SSB time $\zeta$ appears only through the combination $\Omega_{\oplus}t(\zeta)$. The difference between $\Omega_{\oplus}t(\zeta)$ and $\Omega_{\oplus}\zeta$ is negligiblly small. Therefore, in Eq. (38), $G(t(\zeta))$ can be replaced by $G(\zeta)$. The duration ${T}_{\mathrm{seg}}$ is chosen so that $G(t(\zeta))$ can be approximated by a constant in each segment. With this choice of ${T}_{\mathrm{seg}}$, the factor $e^{i\delta\Phi_{\odot}(t)}$ also can be seen as a constant in each segment because it varies slower than the antenna pattern function. Therefore, Eq. (38) can be approximated by

$$h^{\mathrm{STFT}}_{j,k}\simeq h_{0}e^{i\delta\Phi_{\odot}(\zeta_{j})}G(\zeta_{j})W_{k}(\zeta_{j})\,,$$

(39)

where

$$W_{k}(\zeta_{j}):=\frac{1}{T_{\mathrm{seg}}}\int^{\zeta_{j}+T_{\mathrm{seg}}}_{\zeta}d\zeta^{\prime}\ w(\zeta^{\prime}-\zeta_{j})e^{2\pi i\delta f_{k}\zeta^{\prime}}\,.$$

(40)

In this work, we use the tukey window,

	$\displaystyle w(\zeta)=$
	$\displaystyle\begin{cases}\frac{1}{2}-\frac{1}{2}\cos\left(\frac{2\pi\zeta}{\alpha T_{\mathrm{seg}}}\right)\,,&\left(0\leq\frac{\zeta}{T_{\mathrm{seg}}}<\frac{\alpha}{2}\right)\,,\\ 1\,,&\left(\frac{\alpha}{2}\leq\frac{\zeta}{T_{\mathrm{seg}}}\leq 1-\frac{\alpha}{2}\right)\,,\\ \frac{1}{2}-\frac{1}{2}\cos\left(\frac{2\pi(T_{\mathrm{seg}}-\zeta)}{\alpha T_{\mathrm{seg}}}\right)\,,&\left(1-\frac{\alpha}{2}<\frac{\zeta}{T_{\mathrm{seg}}}\leq 1\right)\,.\end{cases}$		(41)

We set the parameter $\alpha$ to 0.125. With $\beta_{k}:=T_{\mathrm{seg}}\delta f_{k}$, Eq. (40) can be calculated as

	$\displaystyle W_{k}(\zeta_{j})=$	$\displaystyle\ e^{2\pi i\delta f_{k}\zeta_{j}}$
		$\displaystyle\times\frac{(1+e^{i\pi\alpha\beta_{k}})(1-e^{2\pi i\beta_{k}(1-\alpha/2)})}{4\pi i\beta_{k}(\alpha^{2}\beta_{k}^{2}-1)}\,.$		(42)

Using the Jacobi-Anger expansion, we can expand the factor $e^{i\delta\Phi_{\odot}(\zeta_{j})}$ that appears in Eq. (39) as

$$e^{i\delta\Phi_{\odot}(\zeta_{j})}=\sum_{\ell=-\infty}^{\infty}i^{\ell}J_{\ell}(X)e^{i\ell\Omega_{\odot}\zeta_{j}}e^{i\ell(\varphi_{\odot}-\phi_{X})}$$

(43)

where $J_{\ell}(z)$ is the Bessel function of the first kind and

	$\displaystyle X:=\frac{2\pi{f}_{\mathrm{gw}}{R}_{\mathrm{ES}}}{c}\sqrt{(\Delta n_{x})^{2}+(\Delta n_{y})^{2}},$		(44)
	$\displaystyle e^{i\phi_{X}}:=\Delta n_{x}+i\Delta n_{y}.$		(45)

Therefore, Eq. (39) can be expressed as

	$\displaystyle{h}^{\mathrm{STFT}}_{j,k}$	$\displaystyle\simeq h_{0}G(\zeta_{j})W_{k}(0)e^{2\pi i\delta f_{k}\zeta_{j}}$
		$\displaystyle\times\sum_{\ell=-\infty}^{\infty}i^{\ell}J_{\ell}(X)e^{i\ell(\varphi_{\odot}-\phi_{X})}e^{i\ell\Omega_{\odot}\zeta_{j}}\,.$		(46)

The Fourier transform of ${h}^{\mathrm{STFT}}_{j,k}$ with a fixed integer $k$ is defined by

$$\mathsf{H}_{\ell,k}:=\frac{1}{{N}_{\mathrm{seg}}}\sum_{j=0}^{{N}_{\mathrm{seg}}-1}{h}^{\mathrm{STFT}}_{j,k}e^{-2\pi ij\ell/{N}_{\mathrm{seg}}}.$$

(47)

We refer to $\mathsf{H}_{\ell,k}$ as the $\ell$-domain signal. To understand the pattern hidden in $\mathsf{H}_{j,k}$, we set aside the factor $G(\zeta_{j})$ for a while. Then, Eq. (47) can be estimated as

$$\mathsf{H}_{\ell,k}\sim h_{0}W_{k}(0)i^{\ell^{\prime}}J_{\ell^{\prime}}(X)e^{i\ell^{\prime}(\varphi_{\odot}-\phi_{X})},$$

(48)

with $\ell^{\prime}\sim\ell+2\pi\Omega_{\odot}^{-1}\delta f_{k}$. Because of the fact that $J_{\ell}(z)\simeq 0$ for $|\ell|\gtrsim|z|$ and $X\lesssim O(10^{3})$ for ${f}_{\mathrm{gw}}=100$Hz, the signal is localized within the region where only few thousand bins in the $\ell$-domain. Putting back the antenna pattern $G(\zeta_{j})$, we expect that $\ell$-domain signals lose their amplitude and their localizations become worse than those for the idealized cases. Figure 2 shows an example of the $\ell$-domain signal.

By the method shown in the previous subsection, for every grid point $\bm{{n}}_{\mathrm{g}}$ and every frequency bin $f_{k}$, we obtain an $\ell$-domain strain defined by

$$\mathsf{S}_{\ell,k}:=\mathsf{H}_{\ell,k}+\mathsf{N}_{\ell,k},$$

(49)

with

$$\mathsf{N}_{\ell,k}:=\frac{1}{{N}_{\mathrm{seg}}}\sum_{j=0}^{{N}_{\mathrm{seg}}-1}{n}^{\mathrm{STFT}}_{j,k}e^{-2\pi ij\ell/{N}_{\mathrm{seg}}}.$$

(50)

There are ${T}_{\mathrm{obs}}/{T}_{\mathrm{seg}}\sim O(10^{6})$ data points in a single $\ell$-domain strain and we know that the signal in $\ell$-domain will be localized within a small region $\sim O(10^{3})$. Thus, the excess power method Anderson et al. (2001) is useful for selecting the candidates with a minimal computational cost. We here divide an $\ell$-domain signal into short chunks so that each chunk has the length $\delta\ell$ and neighbored segments have an overlap by $\delta\ell/2$, which is one of the simplest choices but not the optimal one. Then, we obtain ${N}_{\mathrm{chunk/signal}}=2({N}_{\mathrm{seg}}-\delta\ell)/\delta\ell$ chunks from one $\ell$-domain signal. The excess power statistic for the grid point ${\bm{{n}}}_{\mathrm{g}}$, the frequency bin $f_{k}$, and the $c$-th chunk ($c=0,1,\dots,{N}_{\mathrm{chunk/signal}}-1$) is defined by

$$\mathcal{E}({\bm{{n}}}_{\mathrm{g}},f_{k},c):=4\sum_{\ell=c\delta\ell/2}^{(c+2)\delta\ell/2-1}\frac{|\mathsf{S}_{\ell,k}|^{2}}{\tilde{\sigma}_{k}^{2}}\,,$$

(51)

where

$$\langle\mathsf{N}_{\ell,k}\mathsf{N}^{\ast}_{\ell^{\prime},k}\rangle=:\frac{1}{2}\tilde{\sigma}_{k}^{2}\delta_{\ell\ell^{\prime}}\,.$$

(52)

The variance of noise in $\ell$-domain, $\tilde{\sigma}_{k}$, is estimated as

$$\tilde{\sigma}_{k}^{2}=\frac{{S}_{\mathrm{n}}(f_{k})}{{N}_{\mathrm{seg}}{T}_{\mathrm{seg}}}\times\mathcal{W}\,,$$

(53)

where $\mathcal{W}$ is the factor coming from the window function and defined by

$$\mathcal{W}:=\int^{0.5}_{-0.5}dx[w(x)]^{2}\,.$$

(54)

The derivation of Eq. (53) is summarized in Appendix A.

We define the SNR of the excess power by

$${\rho}_{\mathrm{EP}}({\bm{{n}}}_{\mathrm{g}},f_{k},c):=\frac{\mathcal{E}({\bm{{n}}}_{\mathrm{g}},f_{k},c)-\langle\mathcal{E}\rangle_{\mathrm{n}}}{\sigma_{\mathrm{n}}(\mathcal{E})}\,,$$

(55)

where

$$\langle\mathcal{E}\rangle_{\mathrm{n}}=2\delta\ell\,,$$

(56)

and

$$\sigma_{\mathrm{n}}(\mathcal{E}):=\sqrt{\langle(\mathcal{E}-\langle\mathcal{E}\rangle_{\mathrm{n}})^{2}\rangle_{\mathrm{n}}}=2\sqrt{\delta\ell}\,,$$

(57)

are, respectively, the expectation value and the standard deviation of $\mathcal{E}$ when only noise exists. We select the candidate set of parameter values $\{{\bm{{n}}}_{\mathrm{g}},f_{k},c\}$, when

$${\rho}_{\mathrm{EP}}({\bm{{n}}}_{\mathrm{g}},f_{k},c)>{\hat{\rho}}_{\mathrm{EP}}$$

is satisfied with a threshold value ${\hat{\rho}}_{\mathrm{EP}}$. Strictly speaking, since the excess power statistic $\mathcal{E}$ is the sum of $2\delta\ell$ squared Gaussian random variables with the variance $1/2\sqrt{\delta\ell}$, $\mathcal{E}$ follows a chi square distribution with the degree of freedom $2\delta\ell$. However, since here we choose $\delta\ell$ to be large, the distribution of $\mathcal{E}$ can be approximated by a Gaussian distribution with the average $2\delta\ell$ and the standard deviation $2\sqrt{\delta\ell}$. Therefore, in the absence of gravitational wave signal, the probability distribution of ${\rho}_{\mathrm{EP}}$ is a Gaussian distribution with zero mean and unit variance.

Also in the presence of some signal, the excess power statistics ${\rho}_{\mathrm{EP}}$ is given by a sum of many statistical variables. Thus, the statistical distribution of ${\rho}_{\mathrm{EP}}$ can be approximated by the Gaussian distribution whose mean and standard deviation are calculated as

$${\mu}_{\mathrm{EP}}(\bm{{\xi}})=\frac{2P_{k}(\bm{{\xi}})}{\tilde{\sigma}^{2}_{k}\sqrt{\delta\ell}}\,,$$

(58)

and

$${\sigma}_{\mathrm{EP}}(\bm{{\xi}})=\sqrt{1+\frac{4P_{k}(\bm{{\xi}})}{\tilde{\sigma}_{k}^{2}\delta\ell}}\,,$$

(59)

Here, we define

$$P_{k}(\bm{{\xi}}):=\sum_{\ell}|\mathsf{H}_{\ell,k}(\bm{{\xi}})|^{2}\,,$$

(60)

and we define $\bm{{\xi}}$ as a set of parameters $\bm{{\xi}}:=(h_{0},\vec{\xi})=(h_{0},{f}_{\mathrm{gw}},{\alpha}_{\mathrm{s}},{\delta}_{\mathrm{s}})$. The false alarm rate and the detection efficiency will be assessed with this Gaussian approximation.

After selecting candidates and narrowing down the possible area at which the source is likely to be located, we apply the coherent matched filtering for the follow-up analysis. The grid points with the resolution shown in Eq. (20) are placed to cover the selected area. Assuming a grid point, we can carry out the demodulation of the phase by using the time resampling technique. If the deviation between the directions of the grid point and the source is smaller than the resolution, the residual phase remaining after the time resampling is sufficiently small to avoid the loss of SNR.

In this operation, heterodyning and downsampling can significantly reduce the data length and hence the computational cost Dupuis and Woan (2005). Let us assume that we have a candidate labeled with $\{{\bm{{n}}}_{\mathrm{g}},f_{k},c\}$. If the candidate is the true event, the gravitational wave frequency ${f}_{\mathrm{gw}}$ should take the value in the narrow frequency band indicated by

$$f_{k}-\frac{1}{2{T}_{\mathrm{seg}}}\leq{f}_{\mathrm{gw}}\leq f_{k}+\frac{1}{2{T}_{\mathrm{seg}}}\,.$$

(69)

By multiplying the factor $e^{-2\pi if_{k}\zeta}$ to the resampled strain, we can convert the gravitational wave signal frequency to near DC components (heterodyning). After that, the gravitational wave signal has a lower frequency than $1/2{T}_{\mathrm{seg}}$ Hz. Therefore, downsampling by appropriately averaging the resampled strain data with a sampling frequency $\sim 1/{T}_{\mathrm{seg}}$ reduces the number of data points without loss of the significance of the gravitational wave signal.

The coherent matched filtering follows the heterodyning and the downsampling processes. As stated in Eqs. (2) and (12), we fix $\psi=0$, $\cos\iota=1$, and $\phi_{0}=0$ and treat them as known parameters. Also, we assume that the signal waveform and the template completely match. The definition of a match is already given in Eq. (17). The gravitational waveform is written as

$$h(\bm{{\xi}})=h_{0}\cdot{h}_{\mathrm{temp}}(\vec{\xi})\,.$$

(70)

Among these parameters, the amplitude $h_{0}$ can be analytically marginalized to maximize the likelihood. Then, we obtain the signal-to-noise ratio in Eq. (17) and use it as the detection statistic. When only the detector noise dominates the strain data, ${\rho}_{\mathrm{MF}}$ follows the standard normal distribution. On the other hand, if the signal exists, the SNR follows a Gaussian distribution with a mean

$${\mu}_{\mathrm{MF}}(\bm{{\xi}})=h_{0}\cdot\sqrt{({h}_{\mathrm{temp}}(\vec{\xi})|{h}_{\mathrm{temp}}(\vec{\xi}))}\,,$$

(71)

and a unit variance.

Our procedure is characterized by three parameters $\bm{{\rho}}:=(\mathrm{FAP}_{\mathrm{EP}},{r}_{\mathrm{NN}},{\hat{\rho}}_{\mathrm{MF}})$, i.e.,

• •

  $\mathrm{FAP}_{\mathrm{EP}}$, false alarm probability for each chunk,

• •

  ${r}_{\mathrm{NN}}$, the radius of the predicted region to which the follow-up analysis is applied,

• •

  ${\hat{\rho}}_{\mathrm{MF}}$, the threshold of the SNR of the coherent matched filtering.

${\mathcal{N}}_{\mathrm{EP}}$ denotes the computational cost of the excess power method. $2{N}_{\mathrm{seg}}$ multiplications and $2{N}_{\mathrm{seg}}$ additions of real numbers are required to calculate the excess powers for all chunks in one $\ell$-domain signal. The computational cost for calculating the excess powers for all chunks can be estimated as

$${\mathcal{N}}_{\mathrm{EP}}=4{N}_{\mathrm{seg}}\times{N}_{\mathrm{grid}}\times{N}_{\mathrm{bin}}\sim 4.7\times 10^{15}\,,$$

(72)

in the unit of the number of floating point operations. As we see in the following, this cost can be neglected. Next, we check the computational time of the neural network analysis. We estimate the computational time of the neural network by measuring the elapsed time for analyzing ten thousand data. Because the elapsed time is 1.4sec, the total computational time of the neural network is estimated as

$${\mathcal{T}}_{\mathrm{NN}}\simeq\frac{1.4\mathrm{sec}}{10^{4}\mathrm{data}}\times{N}_{\mathrm{candidate}}\,,$$

(73)

where ${N}_{\mathrm{candidate}}$ is the number of candidates which are selected by the excess power method and is estimated as

	$\displaystyle{N}_{\mathrm{candidate}}$	$\displaystyle={N}_{\mathrm{chunk}}\times{\mathrm{FAP}}_{\mathrm{EP}}$
		$\displaystyle={N}_{\mathrm{grid}}\cdot{N}_{\mathrm{bin}}\cdot\frac{{N}_{\mathrm{seg}}}{{(\delta\ell)}_{\mathrm{slide}}}\times{\mathrm{FAP}}_{\mathrm{EP}}$
		$\displaystyle\simeq 4.6\times 10^{10}\left(\frac{\mathrm{FAP}_{\mathrm{EP}}}{10^{-2}}\right)\,.$		(74)

Substituting it, we obtain

$${\mathcal{T}}_{\mathrm{NN}}\simeq 6.4\times 10^{6}\mathrm{sec}\left(\frac{\mathrm{FAP}_{\mathrm{EP}}}{10^{-2}}\right)\,.$$

(75)

Therefore, we focus on the case $\mathrm{FAP}_{\mathrm{EP}}\leq 10^{-2}$. The computational cost of our analysis is dominated mainly by the preprocessing of the observed strain data and the follow-up analysis. The computational cost of the entire analysis is denoted by ${\mathcal{N}}_{\mathrm{comp}}$ and is approximately calculated by

$${\mathcal{N}}_{\mathrm{comp}}={\mathcal{N}}_{\mathrm{preprocess}}+{\mathcal{N}}_{\mathrm{follow-up}}\,,$$

(76)

where ${\mathcal{N}}_{\mathrm{preprocess}}$ and ${\mathcal{N}}_{\mathrm{follow-up}}$ are the computational cost of the preprocessing and the follow-up analysis, respectively. In this work, we fix the STFT segment duration and the length of the chunk. Thus, the computational cost of the preprocessing is a constant.

$${\mathcal{N}}_{\mathrm{preprocess}}={N}_{\mathrm{grid}}\times({\mathcal{N}}_{\mathrm{STFT}}+{\mathcal{N}}_{\mathrm{FFT}}\times{N}_{\mathrm{bins}})\,.$$

(77)

The computational cost of the STFT is

$${\mathcal{N}}_{\mathrm{STFT}}={N}_{\mathrm{seg}}\cdot 5{T}_{\mathrm{seg}}{f}_{\mathrm{s}}\log_{2}{T}_{\mathrm{seg}}{f}_{\mathrm{s}}\,,$$

(78)

and the computational cost of FFT is

$${\mathcal{N}}_{\mathrm{FFT}}=5{N}_{\mathrm{seg}}\log_{2}{N}_{\mathrm{seg}}\,.$$

(79)

With ${T}_{\mathrm{obs}}=2^{24}$sec, ${T}_{\mathrm{seg}}=2^{5}$sec, ${f}_{\mathrm{s}}=2^{10}$Hz, and ${N}_{\mathrm{grid}}\sim 3.5\times 10^{5}$, the computational cost of the preprocess is estimated as

$${\mathcal{N}}_{\mathrm{preprocess}}\simeq 2.3\times 10^{18}\,.$$

(80)

On the other hand, the computational cost of the follow-up analysis is determined by the combination of $\mathrm{FAP}_{\mathrm{EP}}$ and ${r}_{\mathrm{NN}}$. The computational cost of the follow-up analysis is

$${\mathcal{N}}_{\mathrm{follow-up}}={N}_{\mathrm{candidate}}\times\frac{\pi({r}_{\mathrm{NN}})^{2}}{{(\delta\theta)}_{\mathrm{coh}}^{2}}\times{\mathcal{N}}_{\mathrm{FFT,coh}}\,.$$

(81)

Here, ${(\delta\theta)}_{\mathrm{coh}}^{2}$ is the typical area of region where each grid point of the coherent analysis covers (see Eq. (20)). The computational cost of taking match is dominated by the Fourier transform and calculated as

$${\mathcal{N}}_{\mathrm{FFT,coh}}=5({T}_{\mathrm{obs}}{f}_{\mathrm{s,coh}})\log_{2}({T}_{\mathrm{obs}}{f}_{\mathrm{s,coh}})\simeq 5.0\times 10^{7}\,,$$

(82)

where ${f}_{\mathrm{s,coh}}=1/{T}_{\mathrm{seg}}=2^{-5}$Hz. Therefore, we estimate the computational cost of the follow-up analysis as

$${\mathcal{N}}_{\mathrm{follow-up}}=9.0\times 10^{22}\left(\frac{{r}_{\mathrm{NN}}}{10^{-3}\mathrm{rad}}\right)^{2}\left(\frac{\mathrm{FAP}_{\mathrm{EP}}}{10^{-2}}\right)\,.$$

(83)

Substituting Eqs. (80) and (83) into Eq. (76), we can assess the computational cost of the entire analysis as a function of ${r}_{\mathrm{NN}}$ and $\mathrm{FAP}_{\mathrm{EP}}$. Figure. 3 shows the computational cost for various combinations of $\mathrm{FAP}_{\mathrm{EP}}$ and ${r}_{\mathrm{NN}}$. One can read a feasible combination of $\mathrm{FAP}_{\mathrm{EP}}$ and ${r}_{\mathrm{NN}}$ depending on one’s available computational resources.

The false alarm probability of the entire process (see Wette (2012); Dreissigacker et al. (2018)) is

	$\displaystyle{p}_{\mathrm{fa}}(\bm{{\rho}})=$	$\displaystyle\left\{1-\left(\mathrm{Prob}\left[{\rho}_{\mathrm{EP}}<{\hat{\rho}}_{\mathrm{EP}}\mid{\rho}_{\mathrm{EP}}\sim\mathcal{N}(0,1)\right]\right)^{{N}_{\mathrm{chunk}}}\right\}$
		$\displaystyle\times\left\{1-\left(\mathrm{Prob}\left[{\rho}_{\mathrm{MF}}<{\hat{\rho}}_{\mathrm{MF}}\mid{\rho}_{\mathrm{MF}}\sim\mathcal{N}(0,1)\right]\right)^{{N}_{\mathrm{t}}}\right\}\,,$		(84)

where ${N}_{\mathrm{t}}$ is the number of required templates for the coherent search. It can be estimated by

$${N}_{\mathrm{t}}={N}_{\mathrm{candidate}}\times\frac{\pi({r}_{\mathrm{NN}})^{2}}{(\delta\theta^{2})_{\mathrm{coh}}}\times{N}_{\mathrm{bin,coh}},$$

(85)

where ${N}_{\mathrm{bin,coh}}$ is the number of the frequency bins of the coherent search. Using the value listed in Table. 2, we obtain

$${N}_{\mathrm{t}}\simeq 5.6\times 10^{21}\left(\frac{{r}_{\mathrm{NN}}}{10^{-3}\mathrm{rad}}\right)^{2}\left(\frac{\mathrm{FAP}_{\mathrm{EP}}}{10^{-2}}\right)\,.$$

(86)

Because the false alarm probability of the follow-up stage determines that of the entire process, we can approximate it as

$${p}_{\mathrm{fa}}(\bm{{\rho}})\simeq\left\{1-\left(\mathrm{Prob}\left[{\rho}_{\mathrm{MF}}<{\hat{\rho}}_{\mathrm{MF}}\mid{\rho}_{\mathrm{MF}}\sim\mathcal{N}(0,1)\right]\right)^{{N}_{\mathrm{t}}}\right\}\,.$$

(87)

Furthermore, because ${N}_{\mathrm{t}}\gg 1$, we can approximate

$${p}_{\mathrm{fa}}(\bm{{\rho}})\simeq{N}_{\mathrm{t}}\cdot\mathrm{Prob}\left[{\rho}_{\mathrm{MF}}>{\hat{\rho}}_{\mathrm{MF}}\mid{\rho}_{\mathrm{MF}}\sim\mathcal{N}(0,1)\right]\,.$$

(88)

In this work, the threshold ${\hat{\rho}}_{\mathrm{MF}}$ is chosen so that the false alarm probability of the entire process is 0.01. As shown in Eq. (86), the number of templates depends on $({r}_{\mathrm{NN}},\mathrm{FAP}_{\mathrm{EP}})$, and the same is true for ${\hat{\rho}}_{\mathrm{MF}}$.