Banks of templates for directed and all-sky narrow-band searches of continuous gravitational waves from spinning neutron stars with several spindowns

We construct efficient banks of templates suitable for searches of continuous gravitational waves from isolated spinning neutron stars. We assume that the search algorithm is based on the time-domain maximum-likelihood $\mathcal{F}$-statistic and we consider narrow-band searches with several spindown parameters included. Our template banks are suitable for both all-sky and directed searches and they enable the usage of the fast Fourier transform in the computation of the $\mathcal{F}$-statistic as well as, in the case of all-sky searches, efficient resampling of data to barycentric time.


I. INTRODUCTION
In the years 2015-2020 second-generation ground-based gravitational-wave (GW) detectors Advanced LIGO [1] and Advanced Virgo [2] completed three observational runs O1, O2, and O3, during which ninety GW transient signals from coalescence of pairs of black holes or neutron stars were registered [3][4][5][6].Many other possible astrophysical sources of gravitational waves remain undetected, among them long-lasting almost monochromatic gravitational waves coming from single rotating neutron stars located in our Galaxy (sources of continuous gravitational waves were recently reviewed in [7][8][9]).There exist several different data-analysis pipelines developed to search for this kind of continuous GW signals [10].Our paper discusses construction of efficient bank of templates needed to perform searches using the time-domain F-statistic pipeline.The F-statistic introduced in [11] is a detection statistic that follows from the maximum-likelihood principle of detection of deterministic signals buried in detector's noise (the maximum-likelihood detection in Gaussian noise is described e.g. in [12][13][14]).Implementation of the maximum-likelihood approach to all-sky narrow-band searches of almost monochromatic GW signals were done in details in the series of articles [11,[15][16][17][18].
In the case of all-sky searches (in which neither the location of the source in the sky nor its GW frequency is known), the time-domain F-statistic pipeline was first used in analysis of Virgo VSR1 data [19], and then also for LIGO/Virgo O1 [20,21], O2 [22], and O3 [23] data.In the search of VSR1 data [19], bank of templates build in section IV of [18] was used, whereas all-sky searches of O1-O3 data [20][21][22][23] employed banks of templates constructed in Ref. [24].Efficient banks of templates for F-statistic-based directed searches (in which the position of the source in the sky is known, but its GW frequency is not) were considered in [25].Our paper generalizes and unifies the results of both Refs.[25] and [24].In Ref. [25], banks of templates in two-dimensional parameter space (spanned by GW frequency and the first spindown parameter) were considered, they were suitable only for directed searches.Reference [24] constructed banks of templates suitable for all-sky searches in 4-dimensional parameter space (spanned by GW frequency, the first spindown, and two more parameters related to the unknown position of the source in the sky).The procedure of construction of template banks presented in section IV of this work can be used both for directed and all-sky searches, and for any number of spindown parameters included.We have implemented this procedure in a computer program that is publicly available (see section V below for more details), and we have tested it for parameter spaces with dimensions from two to six (inclusive).All our template banks enable the usage of the fast Fourier transform algorithm in the computation of the F-statistic and, in the case of all-sky searches, efficient resampling of data to barycentric time.This work completes the research on the construction of banks of templates developed in [24,25].We expect that the bank of templates constructed in the current work will be employed in the near future to analyze data from the LIGO/Virgo/KAGRA detectors [1,2,26] using the time-domain F-statistic pipeline.
As is explained in section IV below, our approach to template placement is based on usage of isoheights of the signal-free autocovariance function of the F-statistic.Another common approach in constructing bank of templates for GW searches is based on the notion of a metric in the space of signal's parameters.This approach was initiated in the works [27,28] for searches of gravitational waves from coalescence of compact binaries.For continuous-wave For the convenience of the reader, we present in this section a short derivation of the signal-free autocovariance function of the F-statistic, which we use to construct bank of templates (more details can be found in section 2 of [24]).We consider narrow-band searches for almost monochromatic GW signals.The searched signal is almost monochromatic in the sense that the modulus of the Fourier transform of the signal's waveform is well concentrated (for frequencies f ≥ 0) around some 'central' frequency f c within the bandwith of the search.The search is also narrow-band, which means that the frequency bandwith of the search is small enough to assume that the spectral density S n of detector's noise can be replaced by a constant S 0 := S n (f c ) within the bandwith.We also assume that the noise n in the detector is an additive, stationary, Gaussian, and zero-mean continuous stochastic process.Then the logarithm of the likelihood function is approximately given by (see, e.g., section 2 of [24] for more details) where x denotes the data from the detector, h is the deterministic signal we are looking for in the data, ⟨ • ⟩ is the time averaging operator defined by Here [t c − T o /2; t c + T o /2] denotes observational interval, so t c is the instant of time in the middle of this interval and T o is the length of observation time.
We employ a simplified linear1 model of the GW signal from a rotating neutron star, which can be written as

.3b)
Here h 0 and Φ 0 are constant amplitude and initial phase of the signal, respectively, and the vector ξ collects s + 3 phase parameters, ξ = (ω 0 , . . ., ω s , α 1 , α 2 ). (2.4) The dimensionless frequency ω 0 and spindown ω k parameters are defined as where f 0 and f (k) 0 (k = 1, . . ., s) is respectively an instantaneous frequency and the kth time derivative of the instantaneous frequency of gravitational wave at the solar system barycenter (SSB), evaluated at t = 0.The parameters α 1 and α 2 depend on the position of the GW source in the sky: α 1 := 2πf 0 (sin α cos δ cos ε+sin δ sin ε), α 2 := 2πf 0 cos α cos δ, where α is the right ascension and δ is the declination of the source, ε is the obliquity of the ecliptic.The functions µ 1 and µ 2 are known functions of time, which depend on the motion of the detector with respect to the SSB (see section 2 of [24] and section II of Ref. [11] for more details).
Let us emphasize that we use the simplified model of the signal defined above only to construct a grid of parameters for which we want to calculate the F-statistic.However, in the F-statistic-based pipeline [19][20][21][22][23], to obtain a numerical value of the F-statistic, we use an exact, not an approximate, signal model.In Ref. [15], it was verified that the linear model reproduces well the covariance matrix, defined as the inverse of the Fisher matrix, for the maximum-likelihood estimators of the exact GW signal, both for directed and all-sky searches.Moreover, this reproduction is better the longer the observation time (provided one uses large enough number of spindown parameters).So, the linear model is a good approximation for computing the Fisher matrix, and in the rest of this work we use the Fisher matrices to construct banks of templates.
After introducing new parameters h 1 := h 0 cos Φ 0 and h 2 := h 0 sin Φ 0 (so that ), the signal (2.3) can be written as (2.6) For observations which last at least several hours and for GW frequencies of the order of tens of Hertz or higher, to a good approximation ⟨sin 2Φ⟩ ∼ = 0 and ⟨cos 2Φ⟩ ∼ = 0. Then the time average ⟨h2 ⟩ approximately equals Substituting Eqs.(2.6) and (2.7) into (2.1)we get the following formula for the log likelihood ratio of the signal (2.6): We also compute, making use of (2.7), the optimal signal-to-noise ratio ρ for the signal (2.6): (2.9) In the next step we maximize ln Λ with respect to the parameters h 1 and h 2 .To do this we solve equations ∂ ln Λ[x; h 1 , h 2 , ξ]/∂h i = 0 (i = 1, 2).Their unique solution h 1 ∼ = 2⟨x sin Φ⟩ and h 2 ∼ = 2⟨x cos Φ⟩ gives the maximumlikelihood estimators of the parameters h 1 and h 2 .To get the F-statistic we replace in (2.8) the parameters h 1 and h 2 by their estimators h 1 and h 2 : To construct efficient banks of templates we consider the expectation value3 of the F-statistic (2.10) in the case when the data x contains some GW signal h, i.e. when x(t) = n(t) + h(t; θ ′ ), where θ ′ = (h ′ 1 , h ′ 2 , ξ ′ ) collects the parameters of the signal present in the data.This expectation value can be written as where ρ(h ′ 0 ) is the signal-to-noise ratio (2.9) computed for the signal h(t; θ ′ ) (so 2 ) and C 0 (ξ, ξ ′ ) is the autocovariance function of the signal-free F-statistic. 4In section IV of Ref. [16] it was shown that the autocovariance function C 0 for the narrow-band signal of the form (2.6) can be approximated by Because the phase Φ of the signal (2.6) is a linear function of ξ, C 0 depends only on τ := ξ − ξ ′ , so that one can rewrite Eq. (2.12) as (2.13) C 0 attains its maximal value 1 for τ = 0, i.e. for ξ = ξ ′ .The final approximation we employ relies on assuming that |τ | ≪ 1.After expanding the right-hand side of (2.13) in Taylor series around τ = 0 up to terms quadratic in τ , we obtain the approximate autocovariance function C a :5 where Γ is the (s + 3)-dimensional reduced Fisher matrix with elements equal to The consequence of linearity of the signal's phase [see Eq. (2.3b)] is that the elements Γkl do not depend on the values of the parameters τ .

III. REDUCED FISHER MATRICES
In this section we give explicit formulas for the reduced Fisher matrices defined in Eqs.(2.15), which are needed for constructions of template banks suitable for both directed and all-sky searches.For convenience we replace the time t by the dimensionless variable χ and introduce the dimensionless quantity χ c , Then the observational interval t ∈ [t c −T o /2; t c +T o /2] is transformed into the interval χ ∈ [−1/2; 1/2] of unit length.We can also introduce averaging with respect to the variable χ defined as ⟨⟨g(χ)⟩⟩ := For any function f , ⟨f (t)⟩ = ⟨⟨f (t(χ))⟩⟩, where, according to Eq. (3.1), t(χ) = (χ + χ c )T o .

A. Directed searches with three spindown parameters
To save space we first present the reduced Fisher matrix Γdir s=3 for the case of directed searches with three spindown parameters included.In this case the phase Φ of the signal is of the form [see Eq. (2.3b)] where τ = (ω 0 , ω 1 , ω 2 , ω 3 ).The reduced Fisher matrix Γ for the phase Φ defined above reads (the matrix is symmetric, therefore we display below only its diagonal and upper-diagonal elements) The elements of the matrix Γdir s=3 depend only on the dimensionless parameter χ c .

B. Directed searches with five spindown parameters
In the case of directed searches with five spindown parameters included, the phase Φ of the signal is of the form [see Eq. (2.3b)] where τ = (ω 0 , ω 1 , ω 2 , ω 3 , ω 4 , ω 5 ).The reduced Fisher matrix Γdir s=5 for the phase Φ defined above is equal to (we display again only diagonal and upper-diagonal elements of the matrix) where Γdir s=3 is given in Eq. (3.4) and where ) ) ) ) ) In the case of all-sky searches with three spindown parameters included, the phase Φ of the signal is of the form [see Eq. (2.3b)] where τ = (ω 0 , ω 1 , ω 2 , ω 3 , α 1 , α 2 ).The 6-dimensional reduced Fisher matrix Γ with elements Γij (i, j = 1, . . ., 6) given in Eq. (2.15) and for the phase Φ defined in Eq. (3.8) equals where Γdir s=3 is given in Eq. (3.4) and where Γall-sky Γall-sky Γall-sky Γall-sky Γall-sky Γall-sky Γall-sky Γall-sky Γall-sky Γall-sky The elements of the matrix Γall-sky s=3 depend on the parameter χ c and also on the time-dependent functions µ 1 and µ 2 introduced in Eq. (2.3b), which are determined by the motion of the detector with respect to the SSB.

D. Network of detectors
In this subsection, we derive the reduced Fisher matrix that takes into account the information gathered by a network of N GW detectors.To do this we employ the same simplified model of the GW signal we used above for the computation of single-detector Fisher matrices.Thus the signal related with the I-th where h 0I is a constant amplitude and Φ 0I is a constant initial phase of the signal.The phase parameters ξ are common for all detectors, but the phase functions Φ I are in general different for different detectors (they are different for all-sky searches and identical for all detectors in the case of directed searches).The log likelihood ratio for the network, assuming that noises in different detectors are independent stochastic processes, reads: Here the vector x collects all data streams, x = (x 1 , . . ., x N ), and the vector A consists of amplitude parameters, A := (A 1 , . . ., A N ), where A I := (h 0I , Φ 0I ) (I = 1, . . ., N ).After introducing h 1I := h 0I cos Φ 0I and h 2I := h 0I sin Φ 0I , one easily solves equations ∂ ln Λ net /∂h 1I = 0, ∂ ln Λ net /∂h 2I = 0 for maximum-likelihood estimators of h 1I and h 2I ; the solution reads ĥ1I = 2⟨x I sin Φ I ⟩ and ĥ2I = 2⟨x I cos Φ I ⟩.The F-statistic for the network one gets by replacing in (3.12) the amplitude parameters A by their estimators Â: Let us now consider the situation when the data streams x I in all detectors contain the GW signals h I coming from the same astrophysical source.The signals have amplitude parameters A ′ I (different for different detectors) and phase parameters ξ ′ (the same for all detectors): One can then show that the expectation value of the network F-statistic (3.13) reads where ρ I (h ′ 0I ) is the signal-to-noise ratio in the I-th detector, and C 0I (ξ, ξ ′ ) is the signal-free [i.e., defined for x I (t) = n I (t)] autocovariance function of the F-statistic computed for the I-th detector [see Eq. (2.12)]: (3.17) Making use of (3.16), Eq. (3.15) can be rewritten as One can define the signal-free auotcovariance function of the network F-statistic (3.13) as where m net 0 is the signal-free expectation value of F net : We have found that C net 0 (ξ, ξ ′ ) is just the sum of the auotcovariance functions of the F-statistics computed for individual detectors, In the I-th detector the phase Φ I reads (note that the functions µ 1I and µ 2I are different for different detectors): Φ I depends linearly on the parameters ξ, consequently the one-detector autocovariance functions C 0I (ξ, ξ ′ ) depend on τ := ξ − ξ ′ and can be written as Next we replace the functions C 0I (τ ) by the approximation derived in Eq. (2.14), so where ΓI is the (s + 3)-dimensional reduced Fisher matrix of the Ith detector with elements equal to It is clear that the right-hand side of Eq. (3.18) can not be directly expressed by the network autocovariance function C net 0 , because of the way it depends on the individual amplitudes h ′ 0I of the GW signals and on noise spectral densities S nI (f c ) of particular detectors.In the next step we replace in Eq. (3.18) all amplitudes h ′ 0I by the same for all detectors averaged amplitude h ′ 0 av .This procedure can be justified as follows.The signal-to-noise ratio (SNR) ρ I (h ′ 0I ) is computed here for the simplified model (3.11) of the GW signal and it represents some effective or averaged SNR.The SNR for non-simplified GW signal actually observed in the detector was studied in detail in section III C of [11].The "exact" SNR ρ I for each individual detector obviously depends on the observation time T o and on the parameters of both the GW source and the detector.It depends on right ascension α and declination δ of the GW source, polarization angle ψ of the wave, and inclination angle ι, i.e. the angle between the total angular momentum of the star and the direction from the star to the Earth; it also depends on the position of the detector on Earth, and on orientation of its arms with respect to local geographical directions; finally it depends on the spectral density S nI (f c ) of the detector's noise.One can average the squared of the SNR with respect to the angles α, δ, ψ, and ι.The result is given in Eqs. ( 92) and (93) in [11] and it has the form where h av 0 does not depend on the position of the detector on Earth and orientation of its arms, so it is the same for all detectors.
After replacing in Eq. (3.18) all amplitudes h ′ 0I by the averaged amplitude h ′ 0 av one gets Making use of the approximation (3.24), the above equation becomes (τ = ξ − ξ ′ and we employ matrix notation): It is convenient to introduce the harmonic mean of the detectors' noise spectral densities (evaluated at the same central frequency f c ), (3.29) Making use of the above definition Eq. (3.28) can be rewritten as This equation suggests that the reduced Fisher matrix for the network of detectors can be defined as Let us observe that if the reduced Fisher matrices ΓI for individual detectors are identical for all detectors (this is the case for directed searches, see sections III A and III B above), i.e.ΓI = Γ for I = 1, . . ., N , then the reduced Fisher matrix for the network equals the one-detector Fisher matrix, Γnet = Γ.

A. Covering problem
Construction of a bank of templates begins with the fixing of the minimum value C min of the autocovariance function C 0 (ξ, ξ ′ ) of the signal-free F-statistic.One then looks for such a grid (i.e. a discrete set) of points in the space spanned by the parameters ξ, that for any GW signal with parameters ξ ′ that may be present in the data, there exists a grid node with parameters ξ such that C 0 (ξ, ξ ′ ) ≥ C min .We employ the approximation C 0 (ξ, ξ ′ ) ∼ = C a (ξ, ξ ′ ), so by virtue of Eq. (2.14) this inequality can be written as For the fixed ξ and 0 < C min < 1, Eq. (4.1) defines an hyperellipsoid in (s + 3)-dimensional Euclidean space with the center at ξ and with the size determined by C min .The Euclidean volume of this hyperellipsoid reads 6 (here d = s + 3) For the fixed (and bounded) region of the parameter space, we want to find such a grid fulfilling the requirement (4.1), which consists of possibly smallest number of points.Therefore the problem of finding the optimal grid is a covering problem, i.e. the problem to cover the d-dimensional Euclidean space R d by the smallest number of identical hyperellipsoids.It is convenient to replace the problem of finding the optimal covering of space by identical hyperellipsoids by the problem of finding the optimal covering of space by identical hyperspheres of unit radius, we do this below in section IV A 2. A thorough presentation of the problem of covering d-dimensional Euclidean space R d with identical hyperspheres (with not necessarily a unit radius) can be found in Chap. 2 of Ref. [38].Let us now recall the basic notions related to this problem.
We consider only d-dimensional lattice coverings, i.e. sets of points which are linear combinations with integer coefficients of some basis vectors (P 1 , . . ., P d ), P i ∈ R d (for i = 1, . . ., d).Let the coordinates of the basis vectors be: P 1 = (P 11 , . . ., P 1d ), . . ., P d = (P d1 , . . ., P dd ).The matrix whose rows are the coordinates of vectors P i , is called a generator matrix for the lattice.For each lattice one can define its fundamental region, which when repeated fills the space with one lattice point in each copy.An example of this is a fundamental parallelotope collecting the points of the form The quality of a covering is given by the covering thickness Θ defined as the average number of hyperspheres that contain a point in the space.For lattice coverings Θ is the ratio (volume of one hypersphere)/(volume of fundamental region).Thickness of lattice covering of R d with identical hyperspheres of radius R reads where M is the generator matrix for the lattice.
For any discrete collection of points Q = {Q 1 , Q 2 , . ..} ⊂ R d one defines its covering radius R, as the least upper bound for the distance from any point x ∈ R d to the closest point of Q, R := sup x∈R d inf Qi∈Q |x − Q i |.The identical hyperspheres of radius R centered at the points of Q will cover R d , and no hyperspheres of radius smaller than R will cover it.

Constraints
We consider searches with very large number of grid points in the parameter space, so the time needed to compute the F-statistic for all grid nodes is long and we want to speed up this computation by employing the fast Fourier transform (FFT) algorithm.As the FFT algorithm computes the values of the discrete Fourier transform (DFT) of a time series for a certain set of discrete frequencies (the Fourier frequencies), it will be possible to use the FFT algorithm if the frequency coordinates of grid points will all coincide with the Fourier frequencies.We will construct grids which fulfill this requirement.
If the data form a sequence of N samples x u , u = 1, . . ., N , with the sampling-in-time period ∆t (so N ∆t = T o ), then the frequency resolution of the DFT is ∆f = 1/(N ∆t) and the resolution of the parameter ω 0 [defined in Eq. (2.5)] is One can improve the frequency resolution of the DFT by adding N FFT − N zeros to the N samples (so one takes the Fourier transform of N FFT data points), then the frequency resolution (4.5) changes to7 We thus need such grid that all nodes can be arranged along straight lines parallel to the ω 0 -axis and the distance between neighboring nodes along these lines must be equal to ∆ω 0 .To fulfil this condition it is enough to require that the first basis vector has components8 In the case of all-sky searches there is another important constraint related to resampling of the data to the so called barycentric time (see e.g.section III D in [11], section V A in [18], [39], and section 6.2 in [19]).Because numerically accurate resampling is computationally demanding, it is crucial to do resampling only once per sky position for all spindown values.Let us arrange the parameters ξ in all-sky searches in such a way that the parameters α 1 and α 2 (related to the location of the source in the sky) are the last two components of the vector ξ, see Eq. (2.4).Then, to fulfill the above-mentioned requirement, it is enough to require that the last two components of the second basis vector vanish: P 2 = (P 21 , P 22 , . . ., P 2d−2 , 0, 0) T . (4.8)

Replacing hyperellipsoid-coverings by hypersphere-coverings
We replace the problem of finding the optimal covering of space by identical hyperellipsoids by the problem of finding the optimal covering of space by identical hyperspheres of unit radius.Let us denote the space of original parameters τ by Ω (Ω ⊂ R d , τ ∈ Ω) and the space of the transformed parameters τ ′ by Ω ′ (Ω ′ ⊂ R d , τ ′ ∈ Ω ′ ).We consider a linear transformation from Ω to Ω ′ defined be some matrix F, We require that (4.9) transforms the hyperellipsoid defined by the constant value C min of the autocovariance function into the hypersphere of unit radius.The hyperellipsoid with its center at the origin of the coordinate system is determined in the Ω space by the equation [see Eq. (4.1)] whereas the hypersphere of unit radius in the Ω ′ space is described by the relation τ ′ ⊤ • τ ′ = 1, which, after making use of (4.9), can be written as Equations (4.10) and (4.11) describe the same hyperellipsoid if and only if where R av := √ 1 − C min is the average radius of the hyperellipsoid (4.10).The elements of the matrix Γ depend on the parameter χ c , therefore Eq. (4.12) implies that the elements of the matrix F depend on the parameters χ c and C min .Additionally, in the case of all-sky searches, the elements of F depend, through the evaluation of time averages needed to compute the elements of the Fisher matrix [see section III C], on the position vector of the detector with respect to the SSB during observational interval.
All considered in section III reduced Fisher matrices Γ are symmetric and strictly positive definite, i.e. τ ⊤ • Γ•τ > 0 for any τ ̸ = 0, therefore Eq. (4.12) can be interpreted as the Cholesky decomposition of the matrix Γ/R 2 av .Consequently there exists the unique upper triangular matrix F, with strictly positive diagonal elements, fulfilling Eq. (4.12).Making use of the explicit formulas for the matrices Γ from section III, one easily checks that for all matrices Γ considered in section III the (1, 1) element of the corresponding matrix F equals

.13)
Now we transform the two constraints (4.7) and (4.8) into the Ω ′ space.We start with transforming the first basis vector P 1 [with Ω-space components given in Eq. (4.7)], Because the matrix F is upper triangular, this leads to where, by means of Eq. (4.13), the length ∆ω ′ 0 of the vector P ′ 1 equals

.16)
The second basis vector P 2 has Ω-space components given in (4.8).Its image in the in Ω ′ space also has its two last components equal to zero: When the basis vectors of the grid in Ω ′ space are found, one transform them into Ω space by means of the inverse matrix F −1 .If S ′ is the generator matrix of the grid in Ω ′ space, then the generator matrix S of the corresponding grid in Ω space can be computed as All generator matrices S ′ we construct below are lower triangular.It implies that also generator matrices S are lower triangular (note that the matrix F −1 is upper triangular so its transpose is lower triangular).

B. Construction of grids
In this subsection we describe constructions of grids, which fulfil the constraints given by Eqs.(4.15)-(4.16)and (4.17).All grids we build in Ω ′ space, they only depend on the value of the parameter ∆ω ′ 0 , that is, through Eq. (4.16), they depend on the values of the search parameters ∆ω 0 and C min .The grids we construct are based on simple deformations of the d-dimensional lattice coverings A ⋆ d of space R d by hyperspheres of unit radius.We construct two families of grids, called S 1 and S 2 .The construction of the S 1 grids boils down to squeezing the A ⋆ d lattice in one appropriately selected direction, while the construction of the S 2 grids is more complex, although it also consists in deforming the fundamental region of the A ⋆ d lattice.The grids S 2 have mostly smaller than the grids S 1 thicknesses for smaller values of the parameter ∆ω ′ 0 .

S1 grids
In this subsection we describe the construction of the family S 1 of grids.We start by finding such a vector q of the lattice A ⋆ d joining the origin of the space R d with a lattice node, the Euclidean length |q| of which is as close as possible to, but not less than ∆ω ′ 0 [∆ω ′ 0 is the length of the vector P ′ 1 , see Eq. ( 4. 15)].The smaller the difference, the better-the lattice thickness after its deformation, i.e. compression, will be closer to the thickness of the lattice A ⋆ d .In the limiting case ∆ω ′ 0 = |q|, the grid S 1 coincides with the lattice A ⋆ d .We will illustrate the procedure of construction of generator matrix for the S 1 grid in more detail in the (d = 5)-dimensional case.We denote coordinates of any point in R 5 by . In 5 dimensions, the generator matrix of the lattice A ⋆ 5 with covering radius equal to 1 can be calculated from Eq. (A4), using Eqs.(A5e) and (A6).The explicit form of this matrix reads To find the vector q, it is convenient to modify the generator matrix of the lattice A ⋆ d in such a way as to search only the half-space defined by nonnegative values of the coordinate ω ′ 0 .In the 5-dimensional case the generator matrix of A ⋆ 5 that provides this reads The matrix (4.20) can be obtained from the matrix M 5,1 by making reflection of the lattice generated by M 5,1 with respect to the hyperplane ω ′ 0 = 0 and then replacing the first basis vector with the opposite vector (let us recall that components of the basis vectors form the rows of the generator matrix).
We look for the vector q in the form of a linear combination of the basis vectors v a of the lattice A ⋆ d : where i 1 , . . ., i d are integers.Further, we require that the chosen vector q together with some d − 1 vectors taken out of d basis vectors v a , form a new basis for the lattice A ⋆ d .This way we construct a new generator matrix with the components of the vector q forming, say, the first raw of the matrix.Next we will rotate the lattice several times to make the vector q parallel to the ω ′ 0 -axis.We will display below the necessary rotation matrices in the 5-dimensional case.
Let us denote the components of the chosen vector q (of the length |q| closest to ∆ω ′ 0 ) in 5 dimensions by q = (q 1 , q 2 , q 3 , q 4 , q 5 ) ⊤ , |q| The new generator matrix of A ⋆ 5 looks like this where (p 1 , . . ., p 5 ) ⊤ , (r 1 , . . ., r 5 ) ⊤ , (s 1 , . . ., s 5 ) ⊤ , and (t 1 , . . ., t 5 ) ⊤ are the components of four out of five basis vectors {v 1 , . . ., v 5 } selected in such a way as to create together with the vector q a new basis for the A ⋆ 5 lattice.We then first rotate the lattice in the (α ′ 1 , α ′ 2 ) plane by the angle β 1 , using the rotation matrix R 1 (β 1 ): where the angle β 1 satisfies the equations We obtain the new generator matrix M 3 5,1 , Next, we perform rotation in the (ω ′ 2 , α ′ 1 ) plane using the rotation matrix where the angle β 2 satisfies the equations We get In the next step we make rotation in the (ω ′ 1 , ω ′ 2 ) plane by the angle β 3 such, that This rotation is described by the matrix and the resulting generator matrix reads Finally, the lattice generated by the matrix (4.32) we rotate in the (ω ′ 0 , ω ′ 1 ) plane using the rotation matrix where the angle β 4 is defined through equations This way we obtain the generator matrix of the lattice A ⋆ 5 in the form Next we apply the above described procedure to the basis vector from the second raw of the matrix (4.35).We rotate the lattice generated by the matrix M 6  5,1 successively in the , and (ω ′ 1 , ω ′ 2 ) planes, to replace the components p ′′ 3 , p ′′ 4 , and p ′′ 5 by zero.We can continue this procedure for the remaining basis vectors and we are always able to obtain a lower-triangular generator matrix, regardless of the dimension d of the grid.After applying the above described procedure to all basis vectors of the lattice A ⋆ 5 in 5 dimensions, we obtain the generator matrix of the form To get the generator matrix for the grid S 1 , we have to deform the lattice A ⋆ 5 generated by the matrix (4.36) by squeezing it in the direction of the ω ′ 0 -axis.We squeeze the lattice by the factor s ≤ 1 equal to The generator matrix for the grid S 1 finally equals where the diagonal matrix S reads s 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 The first two rows of the matrix S 1 are components of the basis vectors which fulfill the two constraints (4.15) and (4.17).By means of Eqs.(4.4) and (A7) one easily computes the thickness of the grid S 1 generated by the matrix S 1 .It equals or, expressing ∆ω ′ 0 by search parameters ∆ω 0 and C min , Constructions of the grids S 1 and S 2 , described in section IV B, we have implemented in the computer program GridsGenerator, freely available on GitHub at https://github.com/apisarski/gridgen.Most of the results discussed in the current section were obtained by means of this program.With this program it is possible to compute generator matrices of grids suitable for both direct and all-sky searches, and with different numbers of spindown parameters included.The program uses Fisher matrices derived in section III, in each case finds the matrix F [see Eq. (4.9)] which translates the problem of covering by hyperellipsoids into the problem of covering by unit hyperspheres, then constructs the generator matrix of covering by unit hyperspheres fulfilling the constraints (4.15) and (4.17), and finally, using the matrix F −1 , translate the generator matrix into the original parameter space [see Eq. (4.18)].
We have constructed the grids S 1 and S 2 for different values of the parameter ∆ω ′ 0 and in different dimensions: the grids S 1 in all dimensions from 2 to 6 and the grids S 2 only in dimensions 2 and 4. In figures 1-5 we depict covering thicknesses Θ of the grids S 1 and S 2 as functions of the quantity ∆ω ′ 0 .Figure 1 shows 2-dimensional grids S 1 and S 2 which can be applied for directed searches with one spindown parameter included.3-dimensional grids S 1 are suitable for directed searches with two spindown parameters included and for all-sky searches with no spindown parameters included, they are depicted in figure 2. Figure 3 shows 4-dimensional grids S 1 and S 2 applicable for directed searches with three spindown parameters included and all-sky searches with one spindown parameter included 9 .Figures 4 and 5 10 depict the grids S 1 in 5 and 6 dimensions, respectively.5-dimensional grids can be applied for directed searches with four spindown parameters included and for all-sky searches with two spindown parameters included, whereas 6-dimensional grids are suitable for directed searches with five spindown parameters included and for all-sky searches with three spindown parameters included.
In figures 1-5, we have marked by horizontal lines the value of the thickness of the A ⋆ d lattice, and in figures 1-3, we have also marked the thickness of the corresponding hypercubic Z d lattice11 (the thicknesses of the Z 5 and Z 6 hypercubic lattices are not shown in figures 4 and 5, respectively, as their values are beyond the displayed range of the vertical axis).The A ⋆ d lattices are the thinnest known lattices in dimensions 2 ≤ d ≤ 5, and they are close to the thinnest known lattices in dimensions 6 ≤ d ≤ 15 (compare table 2.1 in [38] with table 2 in [40]).
Allen in Ref. [36] claims, that the thinnest template bank gives the most constraining upper "strict" limit (strict means that it applies at every point in parameter space), but it does not maximize the expected number of detections under assumptions that both the number of templates and the total volume of the parameter space covered by them, are fixed.The expected number of detections is maximized by grids which are "optimal quantizers" [36].For small mismatches the optimal quantizer minimizes the scale-invariant second moment G of the lattice [defined in Eq. (5.8) of [36]], which measures the average squared distance from the nearest template (or the average expected mismatch).Let us note that the lattices A ⋆ 2 and A ⋆ 3 are also the best lattice quantizers, and from table I of [36] it follows, that the lattices A ⋆ 4 -A ⋆ 6 used by us have G only 1.2%-3.0%larger than the best known quantizers in dimensions 4-6 (the same table shows that A ⋆ d lattices in dimensions 7-15 have G larger by 0.6%-7.8%compared to the best known quantizers).At this point, it should also be emphasized that the grids S 1 and S 2 constructed in this work meet the constraints related to the possibility of using the FFT algorithm and fast resampling to the barycentric time.There is a huge difference, in terms of the computational cost needed to calculate the F-statistic for all grid's nodes, between the situation when the applied data analysis pipeline uses the fulfillment of these constraints and when it does not.This means that for a given total computational cost, the constraint-compliant pipeline will enable searching a much larger portion of the parameter space.Considering this, and noting that our grids, which are generally small deformations of the A ⋆ d lattices, have, firstly, near-optimal thicknesses12 , and, secondly, have values of G also not much different from the best known values 13 , using these grids for data analysis is very well justified.
The grids constructed in this work are, and we expect that they will be employed in the first part of continuous-wave searches using the time-domain F-statistic pipeline [19] (the second part is the search for coincidences between the candidates obtained in the first part).The first part is the coherent search of time-domain but narrowband in frequency data segments.In all-sky F-statistic-based searches of O1-O3 data [20][21][22][23], the data were split into frequency bands 0.25 Hz long.For each band, the data were inverse Fourier transformed to get a time series.In O1-O3 searches, data segments of 2 to 24 sidereal days were analyzed (with longer segments for lower frequency ranges), and after performing the search 10 10 -10 12 candidates were obtained, which were then screened for coincidences.In all these searches, to compute the F-statistic, an exact model of the GW signal with two spindowns was used, and to construct banks of templates, the approximate linear model with also two spindowns was employed.The study performed in Ref. [15] showed that these models were adequate for lengths of observational intervals and frequency/first frequency derivative ranges used in the searches of O1-O3 data.However, it can be expected that longer data segments will be coherently analyzed in future searches, making it necessary to use signal models and template banks with two or more spindowns included.Such banks of templates are constructed in the current work.

FIG. 3.
Covering thicknesses Θ of the 4-dimensional grids S1 (black diamonds) and S2 (grey circles) as functions of the quantity ∆ω ′ 0 .The upper horizontal line corresponds to the thickness of the 4-dimensional hypercubic lattice Z 4 (it equals ∼ = 4.9348) and the lower horizontal line denotes the thickness of the optimal lattice A ⋆ 4 (it equals ∼ = 1.7655).The grids can be applied for directed searches with three spindown parameters included and for all-sky searches with one spindown parameter included.The graph made with black diamonds and the part of the graph made with grey circles for ∆ω ′ 0 ≳ 0.5 are identical to the corresponding graphs in Fig. 4 from work [24] (where they were presented as possible to be used only in all-sky searches with one spindown parameter included).

FIG. 4.
Covering thicknesses Θ of the 5-dimensional grids S1 (black diamonds) as a function of the quantity ∆ω ′ 0 .The lower horizontal line denotes the thickness of the optimal lattice A ⋆ 5 (it equals ∼ = 2.1243).The thickness of the 5-dimensional hypercubic lattice Z 5 equals ∼ = 9.1955, so it is beyond the displayed range of the vertical axis.The grids can be applied for directed searches with four spindown parameters included and for all-sky searches with two spindown parameters included.grid G 2,0 done in Ref. [25].For completeness, this construction is repeated below in section B 1 using the notation from this work.In section B 2, we describe the construction of the grid S 2 in the case ∆ω ′ 0 ∈ [2, 4), not considered in Ref. [25].
1. Grid S2 for ∆ω ′ 0 < 2 We start from constructing three circles of unit radii which all cross each other at some point, see figure 6.We build then the three vectors v 1 , v 2 , and v 3 , all of unit length (|v 1 | = |v 2 | = |v 3 | = 1).The vector v 1 begins at the centre of one of the circles and ends at the point which is common to all three circles, the vectors v 2 and v 3 are constructed in a similar way, see figure 6.The components of the vectors v 1 , v 2 , and v 3 are v 1 = (p, q), v 2 = (p, −q), v 3 = (0, 1), (B1) where p and q are positive numbers fulfilling the condition From figure 6 it can be seen that Because the length of the vector P ′ 0 has to be equal to ∆ω ′ 0 , we get and

FIG. 1 .
FIG. 1.Covering thicknesses Θ of the 2-dimensional grids S1 (black diamonds) and S2 (grey circles) as functions of the quantity ∆ω ′ 0 .The upper horizontal line corresponds to the thickness of the 2-dimensional cubic lattice Z 2 (it equals ∼ = 1.5708) and the lower horizontal line denotes the thickness of the optimal hexagonal lattice A ⋆ 2 (it equals ∼ = 1.2092).The grids can be applied for directed searches with one spindown parameter included.

FIG. 5 .FIG. 6 .
FIG. 5. Covering thickness Θ of the 6-dimensional grids S1 (black diamonds) as a function of the quantity ∆ω ′ 0 .The lower horizontal line denotes the thickness of the lattice A ⋆ 6 (it equals ∼ = 2.5511).The thickness of the 6-dimensional hypercubic lattice Z 6 equals ∼ = 17.4410, so it is beyond the displayed range of the vertical axis.The grids can be applied for directed searches with five spindown parameters included and for all-sky searches with three spindown parameters included. ) C. All-sky searches with three spindown parameters