{'fraction_non_alphanumeric': 0.03561555769848923, 'fraction_numerical': 0.0137790635379835, 'mean_word_length': 4.330820196481608, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Use of Excess Power Method and Convolutional Neural Network in All-Sky Search for Continuous Gravitational Waves

The signal of continuous gravitational waves has a longer duration than the observation period. Even if the waveform in the source frame is monochromatic, we will observe the waveform with modulated frequencies due to the motion of the detector. If the source location is unknown, a lot of templates having different sky positions are required to demodulate the frequency and the required large computational cost restricts the applicable parameter region of coherent search. In this work, we propose and examine a new method to select candidates, which reduces the cost of coherent search by following-up only the selected candidates. As a first step, we consider a slightly idealized situation in which only a single-detector having 100% duty cycle is available and its detector noise is approximated by the stationary Gaussian noise. We combine several methods: 1) the short-time Fourier transform with the re-sampled data such that the Earth motion for the source is canceled in some reference direction, 2) the excess power search in the Fourier transform of the time series obtained by picking up the amplitude in a particular frequency bin from the short-time Fourier transform data, and 3) the deep learning method to further constrain the source sky position. We compare the computational cost and the minimum amplitude of the detectable signal with the coherent matched filtering analysis. With a reasonable computational cost, we find that our method can detect the signal having only 32% larger amplitude than that of the coherent search with 95% detection efficiency.

The signal of continuous gravitational waves has a longer duration than the observation period. Even if the waveform in the source frame is monochromatic, we will observe the waveform with modulated frequencies due to the motion of the detector. If the source location is unknown, a lot of templates having different sky positions are required to demodulate the frequency and the required large computational cost restricts the applicable parameter region of coherent search. In this work, we propose and examine a new method to select candidates, which reduces the cost of coherent search by following-up only the selected candidates. As a first step, we consider a slightly idealized situation in which only a single-detector having 100% duty cycle is available and its detector noise is approximated by the stationary Gaussian noise. We combine several methods: 1) the short-time Fourier transform with the re-sampled data such that the Earth motion for the source is canceled in some reference direction, 2) the excess power search in the Fourier transform of the time series obtained by picking up the amplitude in a particular frequency bin from the short-time Fourier transform data, and 3) the deep learning method to further constrain the source sky position. We compare the computational cost and the minimum amplitude of the detectable signal with the coherent matched filtering analysis. With a reasonable computational cost, we find that our method can detect the signal having only 32% larger amplitude than that of the coherent search with 95% detection efficiency.

I. INTRODUCTION
Advanced LIGO and Advanced Virgo detected the first event of gravitational waves from a binary black hole merger in 2015 [1]. After the three observation runs, a lot of binary coalescence events are found [2,3]. In addition to Advanced LIGO and Advanced Virgo, KAGRA [4] and LIGO India [5] are planning to join the gravitational wave detector network [6]. The gravitational wave astronomy is expected to get fruitful results for improving our understanding about the astronomical properties of compact objects [7][8][9], the true nature of gravity [10][11][12], the origin of the Universe [13] and so on (see [14] for a review).
All gravitational wave signals which are detected so far have duration O(10 0−2 ) sec, which is much shorter than the observation period. By contrast, we also expect gravitational waves which last longer than the observation period. Such long-lived gravitational waves are called continuous gravitational waves (see [15,16] as textbooks). Continuous gravitational waves are defined by the following three properties: 1) small change rate of the amplitude, 2) almost constant fundamental frequency and 3) duration longer than the observation period. Rotating anisotropic neutron stars are typical candidate sources of continuous gravitational waves. In addition, there are several exotic objects proposed as possible candidates of the sources of continuous gravitational waves ( [17][18][19]).
Continuous gravitational waves are modeled by simpler waveforms than those of coalescing binaries. The parameters characterizing a typical waveform are its amplitude, the initial frequency and the frequency derivatives with time. Although the waveform generated by the source is analytically simple, the effect of the detector's motion makes the data analysis for continuous gravitational waves challenging. The detector's motion causes the modulation in the frequency and it results in the dispersion of the power in the frequency domain. If the source location is a priori known by electromagnetic observations, the modulation can be removed precisely enough. By contrast, for the unknown target search, we need to take correlation of the data with a tremendous amount of templates to cover the unknown source location on the sky. This severely restricts the applicability of the all-sky coherent search to long duration strain data. Therefore, semi-coherent methods, in which the strain data is divided into a set of segments and statistics calculated for respective segments are summed up appropriately, are often used. Various semi-coherent methods (e.g. Time-Domain F-statistic [20], SkyHough [21], Fre-quencyHough [22]) were proposed so far, and they are actually used to analyze LIGO and Virgo's data. Despite of tremendous efforts, up to now no continuous gravitational wave event is detected [23][24][25].
As another trend of the research, the deep learning method is introduced to the field of gravitational wave data analysis. After the pioneering work done by George & Huerta [26], there are many proposals to use deep learning for wide purposes, e.g., parameter estimation for binary coalescence [27][28][29][30], noise classification [31] and waveform modeling [32]. As for applications to the search for continuous gravitational waves, several groups already proposed deep learning methods. Dreissigacker et al. [33,34] applied neural networks to all-sky search of signals with the duration 10 5 sec and 10 6 sec. They used Fourier transformed strain as inputs. Their methods can be applied to the signal located in broad frequency bands and is extended to the case of multiple detectors and realistic noise. Also, it is shown that the synergies between the deep learning and standard methods or other machine learning techniques are also powerful [35,36].
In this paper, we propose a new method designed for detecting monochromatic waves, combining several transformations and the deep learning method. In Sec. II, the waveform model and some assumptions are introduced. The coherent matched filtering and the time resampling technique are briefly reviewed in Sec. III. In Sec. IV, we explain our strategy that combines several traditional methods such as the resampling, the shorttime Fourier transform and the excess power search with the deep learning method. After showing our setup in Sec. V, we show the results of the assessment of the performance of our new method in Sec. VI. Sec. VII is devoted to the conclusion.

II. WAVEFORM MODEL
We consider a monochromatic gravitational wave. We denote by f gw its frequency constant in time. With the assumption that the source is at rest with respect to the solar system barycenter (SSB), a complex-valued gravitational waveform in the source frame h source (τ ) will be simply written as where τ is called SSB time, h 0 and φ 0 are the amplitude and the initial phase, respectively. In this work, for simplicity, φ 0 is fixed to zero. The phase of a gravitational wave is modulated due to the detector motion and the modulation depends on the source location. The normal vector pointing from the Earth's center to the sky position specified by a right ascension α and a declination angle δ is defined by with the tilt angle between the Earth's rotation axis and the orbital angular momentum . Here, we work in the SSB frame, in which the z-axis is along the Earth's orbital angular momentum and the x-axis points towards the vernal equinox. Defining the detector time t so as to satisfy we obtain the waveform in the detector frame with Φ(t) = 2πf gw t + 2πf gw r(t) · n(α s , δ s ) c . (2.5) A subscript "s" indicates the quantity related to the gravitational wave source. Namely, (α s , δ s ) means the sky position of the source. In the following, we use the notation n s := n(α s , δ s ).
For the modeling of the detector motion, we adopt a little simplification, which we believe will not affect our main result. We assume that the position vector of the detector can be written by a sum of the Earth's rotation part r ⊕ (t), and the Earth's orbital motion part r (t). The Earth is assumed to take a circular orbit on xyplane. Then, we can write r (t) as where R ES , Ω and ϕ are the distance between the Earth and the Sun, the angular velocity of the orbital motion and the initial phase, respectively. The detector motion due to the Earth's rotation can be described as where R E , λ, Ω ⊕ and ϕ ⊕ are the radius of the Earth, the latitude of the detector, the angular velocity of the Earth's rotation and the initial phase, respectively. The modulated phase Φ(t) can be decomposed into (2.10) Finally, we take into account the amplitude modulation due to the detector's motion, which can be described by the antenna pattern function. In this work, the inclination angle and the polarization angle are assumed to be zero for simplicity. Then the gravitational wave to be observed by a detector can be written as (2.12) The definitions of F + (t) and F × (t) are the same as those used in Jaranowski et al., [20]. In this work, the antenna pattern function of LIGO Hanford is employed.
The strain data is written as where n(t) is the detector noise. We assume that the strain data from the detector has no gaps in time and the detector noise is stationary and Gaussian.

III. COHERENT SEARCH METHOD
Before explaining our method, we briefly review the coherent search method and the time resampling technique [20].
If the expected waveforms can be modeled precisely and the noise is Gaussian, the matched filtering is the optimal method for the detection and the parameter estimation, besides the computational cost. The noise weighted inner product is defined by where S n (f ) is the power spectral density of the detector noise. A signal-to-noise ratio (SNR) can be calculated with the inner product between a strain s(t) and a template h temp (t) as Theoretically predicted waveforms h temp (t) have various parameters characterizing the source properties and the geometrical information. A set of waveforms having different parameters is called a template bank. For each template in a template bank, we can assign the value of SNR calculated by Eq. (3.2). If the maximum value of SNR in the template set exceeds a threshold value, it is a sign that an actual signal may exist and the parameter inference is also obtained from the distribution of SNR. Due to a long duration and a narrow frequency band of continuous gravitational waves, the inner product (3.1) can be recast into the time-domain expression as where T obs is the observation time. For a monochromatic source, the waveform can be modeled by Eq. (2.1).
The modulation due to the detector motion is only in the phase of the waveform. The time resampling technique nullifies the phase modulation by redefining the time coordinate. If the position of the source is a priori known, the exact relation (2.3) can be obtained. Therefore, the phase modulation can be completely removed and the monochromatic waveform is applicable to the time-domain matched filtering (3.3). Calculation of the inner product (3.3) between the resampled signal and a monochromatic waveform is equivalent to the Fourier transform. Thus, the fast algorithm (i.e., Fast Fourier transform) can be employed to rapidly search the gravitational wave frequency, f gw . When the source location is unknown, we need to search all-sky by placing a set of grid points {n to keep the maximum loss of SNR within the acceptable range. For each grid point, we carry out the Fourier transform after the transformation Here, we omit the superscript (i), for brevity. The necessary number of grid points N grid can be estimated by the angular resolution of gravitational wave sources. The angular resolution of gravitational wave sources can be roughly given by the ratio between the wavelength of the gravitational wave and the diameter of the Earth's orbit, i.e., Here, we adopt 100Hz as the fiducial value for cλ −1 gw . Thus, the required number of grid points is, at least, The time resampling and the Fourier transform are applied to each grid point. The number of floating point operations required for carrying out FFT is ∼ 1.7 × 10 12 per grid point, with the signal of a duration 10 7 sec and a sampling frequency 1024 Hz. Even if we have a 1PFlops machine, the computational time becomes 2.2 × 10 8 sec, which is longer than the signal duration. For this reason, the coherent all-sky search is not realistic even for monochromatic sources yet.

IV. OUR METHOD
A. Subtracting the effect due to the Earth's rotation 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, Φ ⊕ (t). Assuming a representative grid point n g , we can rewrite the phase (2.5) as where and ∆n := n s − n 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 ng max t |δΦ ⊕ (t)| becomes smaller than a threshold δΦ for any source direction n s within the area covered by the grid point n g . To optimize the grid placement, we employ the method proposed in Ref. [37]. The residual phase δΦ ⊕ is expanded up to the first order of ∆α := α s − α g and ∆δ := δ s − δ g . Then, we get Here, the constant term is neglected because it degenerates with the initial phase φ 0 . The maximum value of the residual phase is The grid points are to be determined to satisfy max t |δΦ ⊕ | < δΦ for any source direction. Because the residual phase (4.5) is symmetric under the transformation δ g → −δ g , the placement of grids on the negative δ side can be generated by inverting the sign of the grids on the positive δ side. Therefore, we focus on the case with 0 ≤ δ ≤ π/2.
Since the residual phase depends only on δ at δ g = π/2, a single template can cover the neighbor of δ = π/2. In fact, at δ = π/2, Eq. (4.5) becomes Therefore, the condition max t |δΦ ⊕ | ≤ δΦ gives the lower bound of δ 1 such that the region δ 1 ≤ δ ≤ π/2 can be covered by a signle patch represented by {(α g , δ g ) = (0, π/2)}, to find Plural patches are necessary to cover the strip of a constant δ in the other range. We introduce a 2-dimensional metric corresponding to the residual phase (4.5), (4.8) 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 := α and Y := − log | cos δ|, the metric can be transformed into Along with [37], we can construct the sky patches covering the half-sky region with 0 ≤ δ ≤ δ 1 . Figure 1 shows a part of grid points constructed under the condition δΦ = 0.058 , (4.10) which we adopt throughout this paper. The total number of grid points to cover the whole sky is for f gw = 100Hz. {(α, δ)|δ < δ1} is covered by a single template (αg, δg) = (0, π/2) and the shape of the patch is square on this plane.

B. Modeling the effect due to the Earth's orbital motion
As we choose δΦ to be sufficiently small, we neglect δΦ ⊕ in the following discussion. Then, after subtracting the phase modulation due to the Earth's rotation, the phase of the gravitational wave (4.1) becomes (4.12) We apply the short-time Fourier transform (STFT) to the time-resampled strain, In the rest of the paper, we treat only the time-resampled data. Therefore, without confusion, the time-resampled data in Eq. (4.13) can be denoted by the same character as the original one. The strain is divided into N seg segments having the duration T seg and their start times are denoted by ζ j := jT slide , (j = 0, 1, · · · , N seg − 1). T slide is not necessary to be equal to T seg . The output of STFT with the window function w(ζ) is defined by (4.16) and f k := k∆f = k/T 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. (2.11) can be neglected and Eq. (4.15) can be approximated by In the expression of G(t(ζ)), the SSB time ζ appears only through the combination Ω ⊕ t(ζ). The difference between Ω ⊕ t(ζ) and Ω ⊕ ζ is negligiblly small. Therefore, in Eq. (4.17), G(t(ζ)) can be replaced by G(ζ). The duration T seg is chosen so that G(t(ζ)) can be approximated by a constant in each segment. With this choice of T seg , the factor e iδΦ (t) also can be seen as a constant in each segment because it varies slower than the antenna pattern function. Therefore, Eq. (4.17) can be approximated by In this work, we use the tukey window, We set the parameter α to 0.125. With β k := T seg δf k , Eq. (4.19) can be calculated as . (4.21) Using the Jacobi-Anger expansion, we can expand the factor e iδΦ (ζj ) that appears in Eq. (4.18) as where J (z) is the Bessel function of the first kind and Therefore, Eq. (4.18) can be expressed as The Fourier transform of h STFT j,k with a fixed integer k is defined by We refer to H ,k as the -domain signal. To understand the pattern hidden in H j,k , we set aside the factor G(ζ j ) for a while. Then, Eq. (4.26) can be estimated as with ∼ + 2πΩ −1 δf k . Because of the fact that J (z) 0 for | | |z| and X O(10 3 ) for f gw = 100Hz, the signal is localized within the region where only few thousand bins in the -domain. Putting back the antenna pattern G(ζ j ), we expect that -domain signals lose their amplitude and their localizations become worse than those for the idealized cases. Figure 2 shows an example of the -domain signal.

C. Excess power method for finding candidates
By the method shown in the previous subsection, for every grid point n g and every frequency bin f k , we obtain an -domain strain defined by There are T obs /T seg ∼ O(10 6 ) data points in a singledomain strain and we know that the signal in -domain will be localized within a small region ∼ O(10 3 ). Thus, the excess power method [38] is useful for selecting the candidates with a minimal computational cost. We here divide an -domain signal into short chunks so that each chunk has the length δ and neighbored segments have an overlap by δ /2, which is one of the simplest choices but not the optimal one. Then, we obtain N chunk/signal = 2(N seg − δ )/δ chunks from one -domain signal. The excess power statistic for the grid point n g , the frequency bin f k , and the c-th chunk (c = 0, 1, . . . , N chunk/signal −1) is defined by The variance of noise in -domain,σ k , is estimated as The derivation of Eq. (4.32) is summarized in Appendix . We define the SNR of the excess power by and are, respectively, the expectation value and the standard deviation of E when only noise exists. We select the candidate set of parameter values {n g , f k , c}, when is satisfied with a threshold valueρ EP . Strictly speaking, since the excess power statistic E is the sum of 2δ squared Gaussian random variables with the variance 1/2 √ δ , E follows a chi square distribution with the degree of freedom 2δ . However, since here we choose δ to be large, the distribution of E can be approximated by a Gaussian distribution with the average 2δ and the standard deviation 2 √ δ . Therefore, in the absence of gravitational wave signal, the probability distribution of ρ EP is a Gaussian distribution with zero mean and unit variance.
Also in the presence of some signal, the excess power statistics ρ EP is given by a sum of many statistical variables. Thus, the statistical distribution of ρ EP can be approximated by the Gaussian distribution whose mean and variance are calculated as where we define, The false alarm rate and the detection efficiency will be assessed with this Gaussian approximation.

D. Neural network for localizing
Deep learning is one of the approaches for finding features being hidden in the data. Artificial neural networks (ANNs) are the architectures playing the central roll in deep learning. An ANN consists of consecutive layers and each layer is formed by a lot of units (neurons). Each layer takes inputs from the previous layer and processed data is passed to the next layer. As a simple example, the process occurring in each layer can be written as the combination of affine transformation and a non-linear transformation, i.e., 39) where x ( ) is a set of input data on the -th layer and g is a nonlinear function, which is called an activation function. We use a ReLU function, defined by (4.40) The parameters w and b are respectively called weights and biases. They are tunable parameters and optimized to capture the features of data. The process to optimize weights and biases is called training. Frequently, the affine transformation and the non-linear transformation are divided into two layers, called a linear layer and a non-linear transformation layer, respectively.
In addition to the layers as given by Eq. (4.39), many variants are proposed so far. In this work, we use also one-dimensional convolutional layers and max-pooling layers. The input of a convolutional layer, denoted by x c i , is a set of vectors. For example, in the case of color images, each pixel has three channels corresponding to three primary colors of light. Therefore, the input data is a set of three two-dimensional arrays. The discrete convolution, which is represented as is calculated in a convolutional layer. Here, x is the input and o is the output data of the layer. C and K are, respectively, the number of channels and the width of the kernel. Each pixel of the data is specified by an index i. The parameters f and b are optimized during the training. A max pooling layer, whose operation can be written as with the kernel size K and the stride s, reduces the length of the data and hence the computational cost.
In supervised learning, a given dataset consists of many pairs of input data and target values. An ANN learns the relation between input data and target values from the dataset and predicts values corresponding to newly given input data. In order to train an ANN, the deviation between the predicted values and the target value is quantified by a loss function. For a regression problem, the mean square loss, is often employed. Here, y and t are a set of predicted values and that of target values, respectively, and they are expressed as d-dimensional vectors. The prediction depends on the weights of the neural network, which are denoted by a single symbol w. An ANN is optimized so as to minimize the loss function for a given dataset, which is the sum of the loss functions for all data contained in the training dataset. Because the complete minimization using all dataset cannot be done, the iterative method is used. The weight w is updated by the replacement algorithm given by where N train is the number of data contained in the dataset and η is called learning rate and characterizes the strength of each update. The algorithm shown in Eq. (4.44) is called gradient descent, which is the simplest procedure to update the weights, and many variants (e.g., momentum, RMS prop, Adam) are proposed so far. Regardless of the choice of the update algorithm, the gradients of a loss function is required and they can be quickly calculated by the backpropagation scheme. In Eq. (4.44), all data in the dataset are used for each iteration. In practice, the loss function for a subset of the dataset is calculated. The subset is called a batch and chosen randomly in every iteration. This procedure is called a mini-batch training.
In the training process, we optimize a neural network so that the loss function is minimized for a dataset. However, this strategy cannot be straightforwardly applied to practical situations. First, the trained neural network may fall in overfitting. Then, the neural network does not have an expected ability to correctly predict the label for a newly given input data which is not used for training. Second, we have to optimize the neural network model and the update procedure, too. For this purpose, we have to appropriately select the hyperparameters, such as the number of neurons of the -th layer (N ( ) ) and the learning rate (η). They are not automatically tuned during the training process.
To solve these problems, we prepare a validation dataset which is independent from the training dataset. The weights of the neural network are optimized so that the loss function for the training dataset is minimized. The validation data is used for monitoring the training process and assessing which model is better for the problem that the user wants to solve. To prevent the overfitting, the training should be stopped when the loss for the validation dataset tend to deviate from that for the training dataset (early stopping). To optimize the hyperparameters, many neural network models having various structures are trained with different training schemes. Among them, we choose the one performing with the smallest loss for the validation dataset.
The whole architecture of the neural network we used is shown in Table I. The input data of the neural network is the complex valued numbers taken from a short chunk of the -domain signal, and the output is the predicted sky position. The -domain waveform H ,k is determined mainly by the residual phase δΦ , which depends on the sky position (α s , δ s ) through the vector ∆n. Because z = 0, only x and y components of ∆n affect δΦ . Therefore, we label each waveform with the values of ∆n x and ∆n y , which are the targets of the prediction of the neural network.
We apply the neural network to each candidate, selected by the excess power method, in order to narrow down the possible area in which the source is likely to be located. As we mentioned, in order to train the neural network, we need the train and the validation datasets. In Sec. V, we show how these datasets are generated.

E. Heterodyning and down sampling
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. (3.5) 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 down sampling can significantly reduce the data length and hence the computational cost [39]. Let us assume that we have a candidate labeled with {n g , f k , c}. If the candidate is the true event, the gravitational wave frequency f gw should take the value in the narrow frequency band indicated by By multiplying the factor e −2πif k ζ to the resampled strain, we can convert the gravitational wave signal frequency to near DC components (heterodyning). After that, the gravitational wave signal has lower frequency than 1/2T seg Hz. Therefore, downsampling by appropriately averaging the resampled strain data with a sampling frequency ∼ 1/T seg reduces the number of data points without loss of the significance of the gravitational wave signal.

V. SETUP OF PERFORMANCE TEST
We analytically estimate the performances of the excess power method and the follow-up search by the coherent matched filtering method. The deep learning method is tested by using mock data, which is generated by injecting an -domain waveform given in Eq. (4.26) into a simulated Gaussian noise. We assume that we use only a single detector, and use the geometry information (e.g. latitude of the detector) of LIGO Hanford in calculating the antenna pattern function as an example. In this work, we focus on one sky patch covered by a single grid point and a frequency bin fixed at f k = 100 Hz, since the scaling to the search over the whole sky and the wider frequency band is straightforward. The sources are randomly distributed within the sky patch. The parameters β k are randomly sampled from a uniform distribution on [−0.5, 0.5]. The original strain has the duration 2 24 sec and the sampling frequency 1024 Hz. In assessing the ability of our detection strategy, we use the normalized gravitational wave amplitude divided by the amplitude of the noise spectral density at the reference frequency, say f ref . Namely, the normalized amplitude is given bŷ Here, we set f ref = f k . In order to train the neural network, we generate 200000 and 10000 waveforms for training and validation datasets, respectively. At each training step, the amplitude whose logarism is randomly chosen from the uniform distribution on −2.3 ≤ log 10ĥ0 ≤ −1.2 is multiplied to the waveforms, and they are injected into the random Gaussian noise. In the same manner, the test dataset containing 10000 waveforms is also prepared and used in Sec. VI C. For training, we employ the mini-batch training. We set the batch size to 256. The Adam [40] is used for the update algorithm. We implement with the Python library PyTorch [41] and use a GPU GeForce 1080Ti. The parameter values we used are listed in Table II.

A. Computational cost for conversion from resampled strains to -domain strains
Before applying the excess power method, we carry out STFT for each resampled strain. Using the values shown in Table II, we estimate the computational cost required in the unit of the number of floating point operations. To obtain the -domain strains, another Fourier transform is carried out to the STFT output in each frequency bin. Therefore, an additional cost, is required for obtaining -domain strains. As a result, the computational cost for the conversion from resampled strains to -domain strains is the sum of Eqs. (6.1) and (6.2).

B. Candidate selection by excess power method
For our test dataset, the average value of P k , which is defined in Eq.(4.38), is By replacing P k in Eqs. (4.36) and (4.37) with P k , the average and the standard deviation of ρ EP can be estimated as First, let us consider the case where the signal is absent. The false alarm probability (FAP) for the excess power method is defined by the probability that the excess power of a chunk obtained from a noise data exceeds the threshold amplitude: where Prob[A] is the probability that the event A occurs. By using the statistics of ρ EP shown in Sec. IV C, the FAP of the excess power method with the thresholdρ EP is calculated by The threshold value with the specified value of the FAP can be derived by inverting the above equation. Figure 3 shows the relation between the FAP and the threshold valueρ EP . We leave FAP EP as a variable to be determined later in Sec. VI E. Next, we consider the case where the signal with the amplitudeĥ 0 exists. The false dismissal probability (FDP) with the thresholdρ EP is defined by From the discussion in Sec. IV C, the FDP can be calculated by Here, we need a careful treatment of the detection efficiency. In our analysis, the excess power method is followed by the deep learning analysis, which further restricts the source location of candidates. Let us assume the excess power method can detect the signal having an amplitude h * with 100x-% efficiency and the sky region estimated by the neural network corresponds to 100y-% credible region. Combining them, we can estimate the efficiency of our method is 100xy-% for the signal of the amplitude h * . In this work, as a measure of the depth of the search, we focus on the amplitude of the signal that can be detected with 95% efficiency by our entire analysis. It is conservative to set To quantify the ability of the method, the amplitude parameterĥ 95% EP is defined by the minimum value of the normalized amplitudeĥ 0 that can be detected with 95% efficiency by our entire analysis, not only by the excess power method. From the above discussion, we choosê h 95% EP so as to satisfy Prob ρ EP (ĥ 95% EP ) >ρ EP = 97.5% . (6.11) In terms of the FDP, the condition given in Eq. (6.11) can be written as FDP EP (ρ EP ;ĥ 95% EP ) = 2.5%. The minimum detectable amplitude with 95% efficiency with changingρEP. A smaller value of false alarm probability corresponds to a bigger value of the thresholdρEP. Therefore, for example, the dot with the largest value ofρEP is given by setting the false alarm rate to 10 −8 . We can detect smaller signal asρEP decreases.
The amplitude parameterĥ 95% EP depends on how large false alarm rate we accept or, equivalently, how small the thresholdρ EP we are allowed to set. Figure 4 showŝ h 95% EP as a function of the SNR thresholdρ EP . In general, if we require a smaller false alarm rate, the thresholdρ EP must be larger. As a remark, in literature, the sensitivity depth defined by D 95% := (ĥ 95% ) −1 is frequently used. It is easy to convert a normalized amplitude to a sensitivity depth.
We compare the excess power method with the coherent matched filtering, while all-sky search by the coherent method is not realistic because of its huge computational cost. For simplicity, the template waveform, h temp (t) in Eq. (3.2), is set to be equal to the observed waveform h obs (t). When the signal is absent, a SNR of the coherent matched filtering follows a Gaussian distribution with mean 0 and variance 1. Therefore, the threshold ρ MF for specified FAP is given by inverting the following equation, Unlike the excess power method case, the FAP of the coherent matched filtering is specified so that the number of false positives is suppressed to O(1). Because the number of grid points for the coherent matched filtering is ∼ 1.3 × 10 11 as mentioned in Sec. III and the number of frequency bins is f gw T obs 1 × 10 9 , the number of trials is ∼ 1.3 × 10 20 . Therefore, it is reasonable to set where The antenna pattern functions F + (t) and F × (t) have the form and respectively (see Ref. [20]). The coefficients {a n } 4 n=0 and {b n } 4 n=1 depend on the latitude of the detector's site λ, the orientation of the detector's arms with respect to local geographical directions γ and the declination angle δ. Substituting the information of LIGO Hanford to λ and γ and taking the average over δ ∈ [−π/2, π/2], we get With the threshold given by Eq. (6.13), the FDP is given by The parameterĥ 95% MF is defined by For the case of the coherent matched filtering, we fixed the value of FAP MF as shown in Eq. (6.14). Therefore, the parameterĥ 95% MF is no longer a variable. Using the relations shown in Eqs. (6.13), (6.14), (6.21) and (6.22), we obtain log 10ĥ 95% MF = −2.27286 .
(6.23) Table III shows the comparison between the excess power method and the matched filtering. Here, we vary FAP EP from 10 −2 to 10 −8 and the ratio of the amplitude parameters,ĥ 95% EP /ĥ 95% MF , changes from 1.1 to 1.5. This means that the excess power method can detect the signals having the amplitude larger than that of the coherent search only by 15(52)% when we choose the threshold for the excess powerρ EP so that FAP EP = 10 −2 (10 −8 ).
Depending on the FAP of the excess power method, the expected number of candidates is estimated as 24) where N chunk is the number of all chunks obtained from all -signals. An -domain signal contains 2 20 data points, and the windows of chunks are generated by successively sliding by 2 10 points. Therefore, N chunk/signal = 1024 chunks are obtained per grid point in the sky and per frequency bin. With the values N grid and N bin listed on Table II, the number of all chunks can be estimated as  One of the key points of our work is to use a convolutional neural network for the purpose of reducing the computational cost required for the follow-up analysis by narrowing down the possible sky area quickly. Using the test dataset, we assess how accurately the trained neural network can determine the sky position. The error in the estimated parameter Q is defined by where Q pred and Q target are the predicted and the true values, respectively. For simplicity, the (α, δ)-plane is regarded as a two-dimensional Euclidean space and the shape of the predicted region is assumed to be a disk on the (α, δ)-plane. For each data, the origin of the disk is set to the point (α pred , δ pred ). The area of the disk is determined so that the probability that the disk contains the true location of the source is equal to 97.5%. To assess the plausible value of the disk area, we evaluate the empirical distribution of the Euclidean distance between the predicted and the true values of (α, δ), which is given by The distribution of D is expected to depend on the signal amplitude. Employing the test dataset which is prepared in Sec. V, we generate seven datasets by multiplying the different values of the amplitudeĥ 0 and adding the Gaussian noise to the waveforms. Figure 5 shows the cumulative histograms of D for various amplitudes of the test data. We define D * (ĥ 0 ) as the distance D where the cumulative fraction of events exceeds 97.5%. The disk area, which is the area of the 97.5% confidence region, is   Figure 6 shows how large the error circle of the sky position estimated by the neural network is. It shows that the neural network achieves the O(10 −7 ) rad 2 localization. It is smaller by two order of magnitude than the area of the patch 6.36 × 10 −5 rad 2 . Therefore, the number of grids required for the follow-up coherent search is reduced by 10 −2 . It takes 14 sec for processing ten thousand signals with a single GPU (GeForce 1080Ti). The computational time for processing all candidates is estimated as Thus, it is reasonable to require the false alarm probability for the excess power method to be less than 10 −3 . If multiple GPUs are available, the computational time of the deep learning can be reduced by employing them for a parallel computing.

D. Follow-up by coherent matched filtering
The coherent matched filtering is applied to all the candidates selected by the excess power method with the information about the sky position estimated by the convolutional neural network. As stated in Sec. III, the coherent matched filtering is equivalent to the Fourier transform after resampling. Since the process of Fourier transform is computationally most expensive, the computational cost for each candidate is estimated by that for the fast Fourier transform (FFT). the number of floating point operations required for FFT per candidate can be estimated by 31) where N length is the number of data points contained in a heterodyned and downsampled strain. From the discussion in Sec. IV E, it is reasonable to set a sampling frequency of a downsampled strain to T −1 seg = 2 −5 Hz. Because a duration of a strain is set to T obs = 2 24 sec (see Table. II), we get If we accept relatively large number of false alarms in the excess power method, then the signals having smaller amplitudes would be detected. In this sense, we want to set FAP EP to a larger value. However, the computational costs for the deep learning and the follow-up search by the coherent matched filtering drastically increase because the number of candidates selected by the excess power method increases. The value of FAP EP (or a threshold valueρ EP ) is determined by the trade-off between the sensitivity for a smaller signal and the computational cost. Figure 7 shows the computational cost as a function of FAP EP . When FAP EP becomes larger than 10 −5 , the computational cost of the follow-up coherent matched filtering surpasses that of the signal conversion to -domain. It increases as the value of FAP EP increases. On the other hand, when FAP EP is lower than 10 −6 , the computational cost of the signal conversion to -domain dominates the total computational cost and it cannot be reduced by decreasing FAP EP . Therefore, at least, we need computational resources that are enough to finish the signal conversion to -domain within the realistic computational time. Figure 8 shows the relation between the computational cost and the minimum amplitude of the detectable signals. From this figure, We can read how small signal can be detected with the acceptable computational cost. For example, assuming that computational resources of O(1)TFLOPS are available, the computational cost N comp ∼ 10 19 is acceptable with an allowed computational time O(10 7 ) sec. In this case, it is reasonable to set FAP EP ∼ 10 −4 and our method can achieve log 10ĥ 95% 0 = −2.15 with 95% detection efficiency. This value of the amplitude is larger than that of the coherent matched filtering only by 32%.

VII. CONCLUSION
In this paper, we proposed a new method of all-sky search for continuous gravitational waves, combining the excess power method and the deep learning method. The time resampling and the STFT are used for localizing the signal into a relatively small number of elements in the whole data. Then, the excess power method selects the candidates of the grid point in the sky and the frequency bins where the signal likely exists. The deep neural network narrows down the region where to be explored by the follow-up search by two orders of magnitude than the original area of the patch. Before the follow-up coherent search, the heterodyning and the downsampling can be employed for reducing the computational cost. We calculated the computaional cost of our method. Most of the computational costs are spent by the STFT and the follow-up coherent matched filtering search. The computational costs of the excess power method and the deep learning are negligibly small. We compared the detection abilities of our method and the straightforward coherent matched filtering. With a reasonable computational power, the minimum amplitude which can be detected by our method is larger than that of the coherent matched filtering only by 32% with 95% detection efficiency. Our training data, which is used for training the neural network, span the restricted parameter region. Namely, the gravitational wave frequencies of the training data are distributed within the small frequency band centered at 100 Hz of width ±1/(2T seg ) and the source locations are sampled from very narrow regions around the fixed grid point. Nevertheless, we can expect our method can be applied to the all-sky search and the frequency band below 100 Hz. If the gravitational wave frequency becomes lower, the strength of the phase modulation becomes weaker (see Eq. (2.5)). Therefore, even if f gw < 100 Hz, the signal power in -domain would still be concentrated in a narrow region and it can be expected that the efficiency of the excess power method is maintained We can employ the similar discussion also for the dependency on the source location. The power concentration in -domain still be valid even if we take into account the dependency of the source location, while it causes the variation of the signal amplitude. From the above discussion, only slight modifications of the construction of the training data and the structure of our neural network are enough to apply our strategies to all-sky search of monochromatic sources having the frequency lower than ∼ 100 Hz.
In addition to the above points, there are several rooms for improving our method. We fixed various parameters such as the width of the STFT T seg and the length of each chunk δ in a little hand-waving manner. Surveying and optimizing these parameter may improve the detection efficiency of our method. Especially, the sampling frequency when downsampling might reduce the computational cost significantly. As can be seen from Eq. (4.27), the deviation δf k causes the translation of the signal in the -domain. It is expected that we can further constrain the gravitational wave frequency than ∼ T −1 seg . Taking into account this effect, we can set the sampling frequencies of downsampled strains to lower value than that of our current choice. This would result in the further reduction of the computational time of the follow-up coherent search.
In the present paper, we assumed that the Gaussianity and the stationarity of the detector noise and 100% duty cycle. We also simplified the waveform model, e.g., the frequency change df /dt is not incorporated. In spite of these simplifications, the obtained results can be regarded as a proof-of-principle and are enough to convince that our method has the potential for improving the allsky search for continuous gravitational waves with the duration of O(10 7 ) sec. Relaxing these assumptions is beyond the scope of this paper and left as future work.

Acknowledgments
We thank Yousuke Itoh for fruitful discussions. This work was supported by JSPS KAKENHI Grant Num-ber JP17H06358 (and also JP17H06357), A01: Testing gravity theories using gravitational waves, as a part of the innovative research area, "Gravitational wave physics and astronomy: Genesis", and JP20K03928. Some part of calculation has been performed by using a GeForce 1080Ti GPU at Nagaoka University of Technology.

Appendix: Noise statistics in -domain
In general, the power spectral density of a stochastic process n(t) is defined by Therefore, we get