TY - JOUR AU - Covas, P, B AB - ABSTRACT All-sky and directed continuous gravitational-wave searches look for signals from unknown asymmetric rotating neutron stars or from neutron stars of unknown rotational frequency. These searches do not take into account the proper motion of the neutron star, assuming that the loss of signal-to-noise ratio caused by this is negligible and that no biases in parameter estimation are introduced. In this paper, we study the effect that proper motion has on continuous wave searches, and we show for what regions of parameter space (frequency, proper motion, sky position) and observation times this assumption may not be valid. We also calculate the relative uncertainty on the proper motion parameter estimation that these searches can achieve. gravitational waves, proper motions, stars: neutron 1 INTRODUCTION Continuous waves (CWs) are long-lasting and almost monochromatic gravitational waves that can be emitted by rotating neutron stars if they are asymmetric around their rotation axis. These asymmetries can be supported either by elastic or magnetic deformations, as recently summarized in Sieniawska & Bejger (2019). Many searches for CWs have been done in the past, looking both for CWs from known pulsars and from unseen neutron stars in our galaxy, both from known locations such as the Galactic Centre or from all the sky (Abbott et al. 2019b, a; Covas & Sintes 2020). These searches have not reported a CW detection, placing bounds on the maximum gravitational-wave amplitude. The optimal frequentist technique to uncover a signal buried in Gaussian noise is the matched filtering, where the data obtained by ground-based detectors such as Advanced LIGO (Aasi et al. 2015), Advanced Virgo (Acernese et al. 2014), and KAGRA (Akutsu et al. 2019) is correlated with a theoretical waveform. These waveforms are generated after a signal model has been assumed, and when this model does not accurately describe the true waveform the signal may not be found. The typical CW signal model takes into account the Doppler modulation produced by Earth’s rotation and orbit around the Solar system barycentre (SSB), and the spin-down of the neutron star produced by the emission of electromagnetic and gravitational radiation. In order to describe this model, four amplitude parameters (the amplitude h0, initial phase ϕ0, inclination angle ι, and polarization angle ψ) and 3 + s phase parameters (initial frequency f0, Right Ascension α and Declination δ, and s spin-down parameters such as the first- and second-order frequency derivatives f1 and f2) are used. When the neutron star is in a binary system, more parameters are needed in order to take into account the Doppler modulation produced by the motion around the binary barycentre, where for the general case, five additional parameters are needed (three for the circular orbit case). Searches for CWs from known pulsars only need to perform the matched filtering once, since all the phase parameters that describe the waveform (without taking into account the so-called amplitude parameters) are previously known. On the other side, searches for CWs from unknown neutron stars have to calculate the matched filter over many different waveforms, which correspond to different combinations of the unknown phase parameters describing the source. When the signal model does not completely describe the signal (such as when the spin-down of the source is neglected), two different effects will take place: The mismatch (loss of signal-to-noise ratio, S/N) produced by using an incorrect signal model will lower the probabilities of detection. Even if the signal is not missed, the estimated parameters will be somewhat biased, which may difficult further confirmation of the source such as from a complementary electromagnetic detection. There are several physical processes that when unaccounted for may render the usual CW signal model incomplete, such as spin-wandering (Mukherjee, Messenger & Riles 2018), the presence of glitches (Ashton, Prix & Jones 2017), timing noise (Ashton, Jones & Prix 2015), or proper motion. In this paper we aim to quantify the effects produced on CW searches when proper motion is neglected. Neutron stars are known to be high-velocity objects (Hobbs et al. 2005), and their proper motion has been measured only for less than 400 pulsars. Searches from known pulsars take into account the proper motion information if available, while searches for unknown neutron stars do not search over the two parameters that characterize the proper motion, assuming that the mismatch produced by this is negligible. In this study, we derive analytical expressions that are able to estimate the mismatch produced by this assumption, and we find that for high frequencies and integration times longer than a year this may cause a large loss of S/N. Prior to this paper, only one attempt to quantify the ability to measure proper motion by CW searches was reported in (Jaranowski & Królak 1999), where a single example was treated: an integration time of four months, with only a single value of the proper motion and frequency. In that study it was obtained that the median relative uncertainty of the proper motion estimators was around 40 per cent. In this paper, we quantify the ability to measure proper motion with more examples and by using a Bayesian MCMC procedure, instead of the analytical Fisher matrix. This paper is structured as follows: In Section 2, we develop the phase model and its dependence on the proper motion of the neutron star; in Section 3, we give a summary of the measured proper motions of known pulsars; in Section 4, we calculate the loss of S/N and biases when the proper motion parameters are neglected, and we present analytical equations that can predict the mismatch; in Section 5, we show the relative uncertainty that the proper motion estimators can achieve; and, in Section 6, we present our conclusions. 2 PHASE MODEL A standard CW signal is described by the following equation (Jaranowski, Królak & Schutz 1998): $$\begin{eqnarray*} h(t) = h_0[F_+(t,\psi , \hat{n})\frac{1+\cos ^2{\iota }}{2} \cos {\phi (t)} + F_{\times }(t,\psi ,\hat{n}) \cos {\iota } \sin {\phi (t)} ], \nonumber \\ \end{eqnarray*}$$(1) where F+ and F× are the antenna patterns of the detectors (which can be found in Jaranowski et al. 1998) for the two different gravitational-wave polarizations, t is the time at the detector frame, the inclination angle ι is the angle between the total angular momentum vector of the star and the direction from the star to the observer, ψ is the wave polarization angle, ϕ(t) is the phase of the signal, h0 is the amplitude of the signal, and |$\hat{n}$| is the position of the source in the sky. The signal given by equation (1) is described by four amplitude parameters (h0, ι, ψ, ϕ0) and 3 (f0, α, δ) + s phase parameters, where s is the number of spin-down/up parameters. The rotational phase of a neutron star is usually described with a Taylor approximation around a reference time, where the different orders of the approximation represent frequency derivatives that are present due to the emission of electromagnetic and gravitational waves. For most of the known pulsars, only one frequency derivative is needed to describe this phase. We assume that the gravitational-wave phase equals two times the rotational phase, thus being described by $$\begin{eqnarray*} \phi (\tau) = \phi _0 + 2\pi \sum _{k=0}^{s} f_k \frac{(\tau -t_{\rm r})^{k+1}}{(k+1)!}, \end{eqnarray*}$$(2) where we define fk as the (k + 1)th-order time derivative of the phase (divided by 2π) given at reference time tr, while ϕ0 is an initial phase. To relate the phase in the source frame ϕ(τ) to the phase in the detector frame ϕ(t), a timing relation that takes into account relativistic effects is developed in Jaranowski et al. (1998). There it is shown that after performing another Taylor approximation, the most important terms affecting the phase are $$\begin{eqnarray*} \phi (t) \cong \phi _{0} + 2 \pi \sum _{k=0}^{s} f^{\prime }_k \frac{(t-t_{\rm r})^{k+1}}{(k+1) !} + \frac{2 \pi }{c} \hat{n} \cdot \vec{r} \sum _{k=0}^{s-1} f^{\prime }_k \frac{(t-t_{\rm r})^{k}}{k !}, \end{eqnarray*}$$(3) where |$\vec{r}$| is the position of the detector with respect to the SSB and |$\hat{n}$| is the position of the source in the sky, given by |$\hat{n}(t)=[\cos {\alpha (t)}\cos {\delta (t)},\sin {\alpha (t)}\cos {\delta (t)},\sin {\delta (t)}]$| where the two sky coordinates are described by (in equatorial coordinates) $$\begin{eqnarray*} \alpha (t) &= \alpha _0 + \mu _{\alpha }(t-t_{\rm r}), \end{eqnarray*}$$(4) $$\begin{eqnarray*} \delta (t) &= \delta _0 + \mu _{\delta }(t-t_{\rm r}), \end{eqnarray*}$$(5) where α0 and δ0 are the sky positions at reference time tr, and μα and μδ are the proper motions in the Right Ascension and Declination. As explained in Jaranowski et al. (1998), the frequencies |$f_k^{\prime }$| appearing in equation (3) are not equal to the frequencies fk in the source frame, differing by a constant offset. The source vector |$\hat{n}(t)$| can be approximated by a Taylor expansion around tr up to first order in time: $$\begin{eqnarray*} &&\hat{n}(t) \approx \hat{n}(t_{\rm r}) + \hat{\dot{n}}(t_{\rm r})(t-t_{\rm r}) \nonumber \\ &&= [\cos {\alpha _0}\cos {\delta _0},\sin {\alpha _0}\cos {\delta _0},\sin {\delta _0}] \nonumber \\ && + (t-t_{\rm r})[-\mu _{\alpha }\sin {\alpha _0}\cos {\delta _0} - \mu _{\delta }\cos {\alpha _0}\sin {\delta _0}, \nonumber \\ && \mu _{\alpha }\cos {\alpha _0}\cos {\delta _0} - \mu _{\delta }\sin {\alpha _0}\sin {\delta _0}, \mu _{\delta }\cos {\delta _0}] . \end{eqnarray*}$$(6) It can be seen that for values of proper motion smaller than ∼10−14 rad s−1, higher order corrections do not have an important contribution for integration times of the order of a few years, since they include terms such as |$\mu _{\alpha }^2 (t-t_{\rm r})^2$|⁠, which are many orders of magnitude smaller than the first-order terms. The detector position is given by the sum of an Earth barycentre component (assumed to be circular) and the barycentre-to-detector component: |$\vec{r} (t) = \vec{r}_{\rm O}(t) + \vec{r}_{\rm d}(t)$|⁠, described by (in equatorial coordinates): $$\begin{eqnarray*} \vec{r}_{\rm O} (t) &=& R_{\rm ES}[\cos {\left(\phi _{\rm O} + \Omega _{\rm O} (t-t_{\rm r})\right)}, \nonumber \\ && \cos {\epsilon }\sin {\left(\phi _{\rm O} + \Omega _{\rm O} (t-t_{\rm r})\right)}, \sin {\epsilon }\sin {\left(\phi _{\rm O} + \Omega _{\rm O} (t-t_{\rm r})\right)}] \nonumber \\ \vec{r}_{\rm d} (t) &=& R_{\rm E}[\cos {\lambda }\cos {(\phi _{\rm r} + \Omega _{\rm r} (t-t_{\rm r}))}, \nonumber \\ &&\cos {\lambda }\sin {(\phi _{\rm r} + \Omega _{\rm r} (t-t_{\rm r}))},\sin {\lambda }], \end{eqnarray*}$$(7) where ΩO is the orbital angular velocity, Ωr is the rotational angular velocity, RES is the mean distance between the SSB and Earth’s barycentre, RE is the distance between Earth’s barycentre and the detector, ϵ is the ecliptic angle, λ is the latitude of the detector, and ϕO and ϕr are initial phases. For the following analytical calculations, we will assume that the vector |$\vec{r}(t)$| only consists on the first term |$\vec{r}_{\rm O}(t)$|⁠, since RES ≫ RE and the effects produced by the rotational term can be neglected, as shown in the Appendix (although all the codes used to obtain the results for this paper use the complete |$\vec{r}$|⁠). 3 PROPER MOTION OF NEUTRON STARS Pulsars are known to have high spatial velocities, reaching up to 1500 km s−1 (Hui & Becker 2006). These proper motions are measured through electromagnetic detections of neutron stars, mainly by using three different mechanisms: pulsar timing (Edwards, Hobbs & Manchester 2006; Matthews et al. 2016); comparison between sky positions at different epochs (Kaplan et al. 2008); scintillation (Cordes 1987; Reardon et al. 2019). Proper motion has been measured only for 344 pulsars, shown in Fig. 1. It can be seen that |μα| is above 10−14 rad/s for 30 pulsars and above 10−15 rad s−1 for 212 pulsars, while |μδ| is above 10−14 rad s−1 for 18 pulsars and above 10−15 rad s−1 for 203 pulsars. This figure also shows that the majority (but not all) of pulsars with measured proper motion values with gravitational-wave frequencies higher than 50 Hz are located in a binary system. Figure 1. Open in new tabDownload slide Proper motion in Right Ascension and Declination for 344 pulsars. The orange circles show 190 pulsars with gravitational-wave frequency less than 50 Hz, while blue circles show the other 154 pulsars. Pulsars in a binary system are shown with a red point. Data taken from Manchester et al. (2005) and downloaded with Pitkin (2018). Figure 1. Open in new tabDownload slide Proper motion in Right Ascension and Declination for 344 pulsars. The orange circles show 190 pulsars with gravitational-wave frequency less than 50 Hz, while blue circles show the other 154 pulsars. Pulsars in a binary system are shown with a red point. Data taken from Manchester et al. (2005) and downloaded with Pitkin (2018). The low number of pulsars with proper motions greater than 10−14 rad s−1 should be taken with caution, since it is known that these measurements suffer from selection effects that bias them towards neutron stars with lower velocities, as discussed in Chatterjee et al. (2005). Neutron stars (which are usually born in the Galactic plane) with high velocities might be able to move away from the Galactic plane, which is the area where the electromagnetic surveys that try to detect unknown neutron stars are more sensitive to. For this reason, even if the majority of the pulsars in the distribution seem to have values lower than 10−14 rad s−1, neutron stars with higher values cannot be dismissed. The highest velocities shown in Fig. 1 are much higher than the usual space velocities of the progenitor stars of type O and B. Orbital velocities in binary systems only reach up to around 200 km s−1, implying that an additional mechanism needs to be present. Several mechanisms have been proposed in order to explain these high velocities (although the dominant mechanism is not known): A post-natal electromagnetic rocket mechanism (due to asymmetric electromagnetic radiation when the magnetic dipole is displaced from the centre of the star), which requires a high initial rotational frequency (Harrison & Tademaru 1975). Asymmetric radiation during the supernova, either of neutrinos (Lai & Qian 1998), of hydrodynamical nature due to asymmetries in the mass ejection (Lai & Goldreich 2000), or due to an asymmetric explosion of gamma-ray bursts (Cui et al. 2007). Dynamical interactions such as the ones produced with the supermassive black hole in the Galactic Centre or within globular clusters. When detected, a proper motion measurement from a continuous wave signal could help to improve our understanding of the mechanism that gives rise to the observed proper motions, since gravitational waves may access a part of the galactic neutron star population that is hidden from the electromagnetic surveys. For example, a very high transverse velocity measurement can be used to constrain the physics of supernova core collapse. Furthermore, many proper motion measurements from CW detections could help to solve the hypothesis of whether there is a correlation between the proper motion vector and the angular momentum vector (spin-kick alignment) (Noutsos et al. 2012), since the latter can also be measured with CWs (represented by the ψ and ι parameters). This correlation would be helpful in determining a specific mechanism for the production of such high proper motion values, since not all of the previous mentioned mechanisms predict the spin-kick alignment (Lai, Chernoff & Cordes 2001). Moreover, measuring the proper motion of a neutron star allows the estimation of the birth site (after an age estimate is done), which can be used to associate the neutron star with a supernova remnant or a nebula. 4 MISMATCH As previously mentioned, all searches for CWs from unknown neutron stars (both all-sky and from known sky positions such as the Galactic Centre) assume that the proper motion of the source is zero. In this section, we calculate what is the expected loss of S/N due to this assumption. In order to detect a signal, CW searches for unknown neutron stars calculate a detection statistic for all the different templates (combinations of phase parameters) that are searched. This detection statistic sorts the templates by the probability that a true astrophysical signal described by those parameters is present in the data. One of the most used detection statistics is the |$\mathcal {F}$|-statistic, which is the frequentist likelihood ratio maximized over the four amplitude parameters that define the CW signal, firstly developed in Jaranowski et al. (1998). A signal is said to be detected (or saved for follow-up) if the |$\mathcal {F}$|-statistic value for some template is above a certain threshold, calculated from the false alarm probability that is created from the background noise. The expected value of the |$\mathcal {F}$|-statistic is related to the S/N of the signal (Pletsch 2011): $$\begin{eqnarray*} \left\langle 2\mathcal {F} \right\rangle = 4 + \rho ^2_0, \end{eqnarray*}$$(8) where |$\rho ^2_0$| is the squared S/N when there is no mismatch, i.e. when the searched parameters are exactly equal to the astrophysical parameters. Due to prohibitively high computational costs, all-sky and (some) directed searches calculate a semicoherent detection statistic, where the data is separated in N shorter segments and phase coherence is only demanded within each of these segments, but not between them. The expected value of the semicoherent |$\mathcal {F}$|-statistic is (Pletsch 2011) $$\begin{eqnarray*} \left\langle 2\mathcal {\tilde{F}} \right\rangle = 4N + \tilde{\rho }^2_0, \end{eqnarray*}$$(9) where |$\tilde{\rho }^2_0$| is the sum of the squared S/N of each segment. Since the values of the searched parameters will never be exactly equal to the parameters of the astrophysical signal, a fraction of the S/N is not recovered. The mismatch m describes the amount of squared S/N that is lost due to not searching at exactly the signal parameters, and it is given by $$\begin{eqnarray*} m = \frac{\rho ^2_0 - \rho ^2(\Theta _{\rm s},d\Theta)}{\rho ^2_0}, \end{eqnarray*}$$(10) ranging from 0 (fully recovered S/N) to 1 (no recovered S/N), where Θs represents the different signal parameters, such as frequency or sky positions, and dΘ represents the difference between the true and searched parameters. The mismatch lowers the obtained |$\mathcal {F}$|-statistic value, implying that a signal that would be detectable without mismatch may not be recovered. The mismatch can be estimated by doing a Taylor expansion of the likelihood ratio around the signal parameters, where it attains a maximum. Usually, only the second-order term is kept: $$\begin{eqnarray*} m \approx g_{ij} (\Theta _{\rm s}) d\Theta ^i d\Theta ^j + \mathcal {O}(d\Theta ^3), \end{eqnarray*}$$(11) where gij is the parameter space metric (i and j run over the dimensions, given by the number of parameters). This approximated mismatch is unbounded and can be higher than 1, and from previous studies, it is known that this approximation highly overestimates the actual mismatch for mismatches higher than ∼0.3 (Prix 2007), a fact that was also studied in Wette (2016), which further analysed where the metric approximation breaks down. 4.1 Parameter bias and expected mismatch Since all-sky and directed searches do not search over the proper motion parameters, some mismatch will always be present (even when an infinitely fine grid over the other searched parameters is used), but usually it is assumed that this mismatch is much lower than the mismatch produced by the other parameters. A similar situation was discussed in Ashton et al. (2017), where the mismatch produced by the presence of glitches in the signal was studied. As noted there, the template that attains the minimum mismatch will not be located at the true signal parameters, since there will be a shifted template combination that minimizes the effect of the missing proper motion parameters. The minimum mismatch and displaced parameters can be estimated by minimizing the mismatch function given by equation (11) (where we have separated the parameters between searched parameters λ and non-searched proper motion parameters Λ): $$\begin{eqnarray*} m = g_{ij} \Delta \lambda ^i \Delta \lambda ^j + g_{IJ} \Delta \Lambda ^I \Delta \Lambda ^J + 2 g_{iI} \Delta \lambda ^i \Delta \Lambda ^I, \end{eqnarray*}$$(12) where now the lowercase indices only go through the searched parameters and the uppercase indices I and J go from 0 to 1 (the two proper motion parameters). The minimum mismatch will be obtained at non-zero displacements of the searched λ parameters. These displacements can be estimated by minimizing the previous mismatch equation: $$\begin{eqnarray*} \frac{ \partial m}{\partial \Delta \lambda ^i} = 0 \longrightarrow \Delta _{\rm min} \lambda ^j = - g_{ij}^{-1} g_{iI} \Delta \Lambda ^I, \end{eqnarray*}$$(13) $$\begin{eqnarray*} m_{\rm min} = \Delta \Lambda ^I \Delta \Lambda ^J (g_{IJ} - g_{iI} g^{-1}_{ij} g_{jJ}). \end{eqnarray*}$$(14) As mentioned in Ashton et al. (2017), this expression is only valid for displacements that would produce a mismatch lower than ∼0.3, since for higher mismatches, the second-order Taylor approximation is not valid. For example, if the unknown searched parameter was f0, and μα was the unknown non-searched parameter, these expressions would be $$\begin{eqnarray*} m &=& g_{f_0f_0}(\Delta f_0)^2 + g_{\mu _{\alpha } \mu _{\alpha }}(\Delta \mu _{\alpha })^2 + 2 g_{f_0 \mu _{\alpha }} \Delta f_0 \Delta \mu _{\alpha } \nonumber \\ \Delta _{\rm min} f_0 &=& - \frac{g_{f_0 \mu _{\alpha }}}{g_{f_0 f_0}} \Delta \mu _{\alpha } \nonumber \\ m_{\rm min} &=& \left(g_{\mu _{\alpha } \mu _{\alpha }} - \frac{g^2_{f_0 \mu _{\alpha }}}{g_{f_0 f_0}} \right) (\Delta \mu _{\alpha })^2. \end{eqnarray*}$$(15) It can be seen that the mismatch is reduced due to the second negative term of the last equation, as compared to the simple case where Δf0 = 0: $$\begin{eqnarray*} m_N = g_{\mu _{\alpha } \mu _{\alpha }}(\Delta \mu _{\alpha })^2. \end{eqnarray*}$$(16) We remark that these are the minimum mismatches that would be obtained if we searched over an infinitely finely spaced template bank over the frequency, spin-down, and sky positions. In a more realistic scenario, the mismatch will always be bigger than this minimum mismatch. In order to calculate these parameter displacements and minimum mismatches, we use a modified (which includes the two proper motion components) version of the UniversalDopplerMetric code from the LALSuite repository (LIGO Scientific Collaboration 2018), which is able to calculate the metric components by computing equation (87) from Prix (2007). After calculating all the metric components, we can obtain the parameter displacements and the minimum mismatch. The parameter biases are shown in Fig. 2, together with the fraction of minimum mismatch compared to the mN mismatch. We have simulated signals from neutron stars with isotropic sky positions and orientations, with frequencies from 100 to 1500 Hz. It can be seen that the fraction between the minimum and mN mismatches highly depends on the total observing time Tobs. When there is more reduction in the mismatch, the recovered parameters deviate more from the true parameters: The sky positions can differ from the true sky positions by more than five bins. For the 1-yr case, the sky positions are the parameters that are more biased, while for the 2-yr search, the first frequency derivative is more biased. For the directed search, the second frequency derivative has the highest bias. These results have been obtained by setting the reference time tr to the middle of the observation time. The minimum mismatch mmin is an invariant quantity with respect to tr, but the sizes of the parameter bias are greatly incremented when using other reference times, such as the initial or ending times of the observation. We remark the fact that the biases shown in this figure for f0, f1, and f2 are between the recovered value and the modified primed frequencies that appear in equation (3), not between the recovered values and the source-frame frequencies. Figure 2. Open in new tabDownload slide Absolute shift of the different parameters compared to the ratio of minimum mismatch given by equation (14) and simple mismatch mN given by equation (12) when Δλi = 0 for all i. The length of each bin has been defined as |$\sqrt{m/g_{\lambda _i \lambda _i}}$| using m = 0.1 (defining the bin with a lower m value would imply higher values for the displayed shifts in bins). This bin definition comes from the mismatch equation (11), where the covariant terms are neglected (this way of defining the bins has been used in many past studies and searches, such as Abbott et al. (2017). These results are for a fully coherent search, where the reference time has been defined as the middle time tmid = t0 + Tobs/2, where t0 is the starting time. Figure 2. Open in new tabDownload slide Absolute shift of the different parameters compared to the ratio of minimum mismatch given by equation (14) and simple mismatch mN given by equation (12) when Δλi = 0 for all i. The length of each bin has been defined as |$\sqrt{m/g_{\lambda _i \lambda _i}}$| using m = 0.1 (defining the bin with a lower m value would imply higher values for the displayed shifts in bins). This bin definition comes from the mismatch equation (11), where the covariant terms are neglected (this way of defining the bins has been used in many past studies and searches, such as Abbott et al. (2017). These results are for a fully coherent search, where the reference time has been defined as the middle time tmid = t0 + Tobs/2, where t0 is the starting time. The figure also shows than the reduction of minimum mismatch is smaller for directed searches, since the sky position is fixed (although a second spin-down parameter is also searched). This figure clearly shows that biases created by assuming the proper motion to be zero can be much larger than the typical resolution of the search. In order to confirm these calculations, we have compared the obtained results with the mismatch obtained when calculating the |$\mathcal {F}$|-statistic values at both the signal and the displaced parameters, using the lalapps_ComputeFstatistic_v2 code (also part of the LALSuite repository). This procedure has returned the same mismatch results as obtained with the UniversalDopplerMetric code. We have repeated these calculations by varying the frequency and total proper motion components, in order to study at which regions of parameter space will the minimum mismatch mmin exceed a certain threshold value. The results are shown in Fig. 3, where two different plots are shown for two different integration times. Each cell of the plot is made by averaging the results from 100 signals distributed with an isotropic sky position and random amplitude parameters (producing S/Ns between 10 and 1000). It can be seen that the mismatch increases with the frequency and with the total proper motion value, and also with the coherent integration time. The maximum proper motion value in these plots is 2.9 × 10−14 rad s−1, although as discussed previously unknown neutron stars could attain even higher proper motion values. For observing times smaller than 1 yr, no minimum mismatches above 0.01 have been obtained. From these plots, it can be seen that when doing a search with a coherent time longer than a year and not searching the proper motion parameters, there is a non-negligible probability of having a high mismatch and missing a signal for gravitational-wave frequencies greater than ∼600 Hz. At lower frequencies, if the total proper motion is higher than 3 × 10−14 rad s−1, it can be inferred that for coherent integration times longer than 2 yr the mismatch could also be non-negligible. Figure 3. Open in new tabDownload slide These plots show the average (between 100 signals for each cell) minimum mismatch given by equation (14) as a function of frequency and total proper motion for two different coherent times: 1.5 (left-hand panel) and 2 yr (right-hand panel). The reference time for these searches has been selected as the middle of the observing time. Figure 3. Open in new tabDownload slide These plots show the average (between 100 signals for each cell) minimum mismatch given by equation (14) as a function of frequency and total proper motion for two different coherent times: 1.5 (left-hand panel) and 2 yr (right-hand panel). The reference time for these searches has been selected as the middle of the observing time. The mismatches for the semicoherent case are shown in Fig. 4. These results belong to one single cell of Fig. 3, but very similar results with the same scaling are obtained for all other cells. The semicoherent metric components are obtained by averaging the different coherent integrations, where, for each of them, the starting time ti will be different. The figure shows that for less than five segments the mismatch is comparable to the fully coherent case, but when there are more segments the mismatch quickly reduces to negligible amounts. Figure 4. Open in new tabDownload slide Scaling of the minimum mismatch with respect to the number of segments in a semicoherent search. Circles show the maximum mismatch, while crosses show the average result from 100 signals belonging to one cell of Fig. 3 (with 1σ error bars). Blue points show the results for an observation time of 2 yr, while orange points show the results for 1 yr. Figure 4. Open in new tabDownload slide Scaling of the minimum mismatch with respect to the number of segments in a semicoherent search. Circles show the maximum mismatch, while crosses show the average result from 100 signals belonging to one cell of Fig. 3 (with 1σ error bars). Blue points show the results for an observation time of 2 yr, while orange points show the results for 1 yr. The previous calculations show that we expect high mismatches only for observation times longer than a year. Although the past O1 and O2 observing runs have been of approximately four and nine months, respectively, the newest O3 run has lasted for about a year, and future observing runs are planned to be longer than a year, as discussed in Abbott et al. (2018). Furthermore, sometimes data from different observing runs has been combined in order to follow-up candidates from a search, as shown in Abbott et al. (2017). These considerations show that for future analysis an explicit search over the proper motion parameters might be needed in order to safely lower the probabilities of missing a signal. Since high mismatches are obtained only for small number of segments, we expect that an explicit search over the proper motion parameters might only be needed at the last stages of a typical CW follow-up procedure, where the initial stages have a large number of segments and subsequent stages reduce the number of segments. The number of candidates is reduced at each stage of the follow-up, and at the last stages, a very small amount of candidates remains above the threshold. For this reason, the increase in computational cost produced by searching over two extra parameters would not highly increase the final cost of a follow-up procedure, since it would only be needed at the last stages. A caveat of these results is that they have been obtained for duty cycles of 1 (i.e. simulated Gaussian data without gaps), while realistic data from gravitational-wave detectors always has duty cycles smaller than 1. We leave for future work an estimation of the effect that this would have on our results. 4.2 Derivation of the proper motion coherent metric components In order to quickly estimate the mismatch that will be present when the proper motion is neglected, an analytical equation is needed. From the previous subsection, it can be seen that we need to calculate all the components of the parameter space metric that are related to the proper motion, such as |$g_{\mu _{\alpha } \mu _{\alpha }}$|⁠. These metric components can be used to estimate the mismatch or to construct a bank of templates with a desired resolution. The phase metric approximation is used to obtain these components, where the amplitude parameters are taken as constant and only the phase parameters are taken into account (Prix 2007). Within this approximation, the metric components are given by $$\begin{eqnarray*} g_{ij}(\Theta) = \langle \partial _i \phi (\Theta) \partial _j \phi (\Theta) \rangle - \langle \partial _i \phi (\Theta) \rangle \langle \partial _j \phi (\Theta) \rangle , \end{eqnarray*}$$(17) where $$\begin{eqnarray*} \partial _i \phi (\Theta) = \frac{ \partial \phi }{\partial \Theta _i} \quad \text{ and } \quad \langle \phi \rangle = \frac{1}{T} \int _{t_0}^{t_0+T} \phi (t) {\rm d}t, \end{eqnarray*}$$(18) where t0 is the starting time of the integral and Θi represents the different parameters describing the signal. The |$g_{\mu _{\alpha } \mu _{\alpha }}$| metric component is explicitly derived in the Appendix. Here we show all the metric components related to the proper motion, where we only keep the terms that depend on the highest order of the coherent integration time T (when the reference time is defined as tmid): $$\begin{eqnarray*} g_{\mu _{\alpha }\mu _{\alpha }} &\approx &\frac{4\pi ^2 R_{\rm ES}^2 f_0^2 T^2 }{24 c^2} [ \cos ^2{\delta _0}\sin ^2{\alpha _0} + \cos ^2{\epsilon }\cos ^2{\alpha _0}\cos ^2{\delta _0}], \nonumber \\ g_{\mu _{\delta }\mu _{\delta }} &\approx &\frac{4\pi ^2 R_{\rm ES}^2 f_0^2 T^2}{24c^2} [ \cos ^2{\alpha _0}\sin ^2{\delta _0} + \cos ^2{\epsilon }\sin ^2{\delta _0}\sin ^2{\alpha _0} \nonumber \\ & + & \sin ^2{\epsilon }\cos ^2{\delta _0} - 2\cos {\epsilon }\sin {\epsilon }\cos {\delta _0}\sin {\delta _0}\sin {\alpha _0}], \nonumber \\ g_{\mu _{\alpha }\mu _{\delta }} &\approx & \frac{4\pi ^2 R_{\rm ES}^2 f_0^2 T^2}{24c^2} [(1-\cos ^2{\epsilon })\cos {\alpha _0}\sin {\delta _0}\sin {\alpha _0}\cos {\delta _0} \nonumber \\ &+& \cos {\epsilon }\sin {\epsilon }\cos {\alpha _0}\cos ^2{\delta _0}], \nonumber \\ g_{\mu _{\alpha }\alpha _0} &\approx & \frac{4\pi ^2R_{\rm ES}^2f_0^2}{2 \Omega _{\rm O} c^2} \sin {\phi _{\rm O}} \cos {\phi _{\rm O}} \cos {T\Omega _{\rm O}} [\cos ^2{\delta _0}\sin ^2{\alpha _0} \nonumber \\ &-&\cos ^2{\epsilon }\cos ^2{\alpha _0}\cos ^2{\delta _0} - \cos {\alpha _0}\sin {\alpha _0}\cos ^2{\delta _0}\cos {\epsilon } (1 \nonumber \\ &-& 2\sin ^2{\frac{T\Omega _{\rm O}}{2}} \frac{\sin {\phi _{\rm O}}}{\cos {\phi _{\rm O}}} - 2\cos ^2{\frac{T\Omega _{\rm O}}{2}} \frac{\cos {\phi _{\rm O}}}{\sin {\phi _{\rm O}}})], \nonumber \end{eqnarray*}$$ $$\begin{eqnarray*} g_{\mu _{\delta }\alpha _0} &\approx & \frac{4\pi ^2R_{\rm ES}^2f_0^2}{2 \Omega _{\rm O} c^2} \sin {\phi _{\rm O}} \cos {\phi _{\rm O}} \cos {T\Omega _{\rm O}} [\cos {\delta _0}\sin {\delta }\cos {\alpha }\sin {\alpha _0} \nonumber \\ &+ &\cos ^2{\epsilon }\sin {\alpha _0}\cos {\alpha _0}\cos {\delta _0}\sin {\delta _0} - \cos {\epsilon }\sin {\epsilon }\cos {\alpha _0}\cos ^2{\delta _0} \nonumber \\ &+ &(\sin ^2{\alpha _0}\cos {\delta _0}\sin {\delta _0}\cos {\epsilon } - \sin {\alpha _0}\cos ^2{\delta _0}\sin {\epsilon } \nonumber \\ &- &\cos ^2{\alpha }\cos {\delta _0}\sin {\delta _0}\cos {\epsilon }) \nonumber \\ &&(\frac{1}{2} - \sin ^2{\frac{T\Omega _{\rm O}}{2}} \frac{\sin {\phi _{\rm O}}}{\cos {\phi _{\rm O}}} - \cos ^2{\frac{T\Omega _{\rm O}}{2}} \frac{\cos {\phi _{\rm O}}}{\sin {\phi _{\rm O}}})], \nonumber \\ g_{\mu _{\alpha }\delta _0} &=& g_{\mu _{\delta }\alpha _0}, \nonumber \\ g_{\mu _{\delta }\delta _0} &\approx & \frac{4\pi ^2R_{\rm ES}^2f_0^2}{2 \Omega _{\rm O} c^2} \sin {\phi _{\rm O}} \cos {\phi _{\rm O}} \cos {T\Omega _{\rm O}} [\cos ^2{\alpha _0}\sin ^2{\delta _0} \nonumber \\ &- &\cos ^2{\epsilon }\sin ^2{\alpha _0}\sin ^2{\delta _0} - \sin ^2{\epsilon }\cos ^2{\delta _0} \nonumber \\ &+ &(\cos {\alpha _0}\sin {\alpha _0}\sin ^2{\delta _0}\cos {\epsilon } - \cos {\alpha _0}\cos {\delta _0}\sin {\delta _0}\sin {\epsilon }) \nonumber \\ &&(1 - 2\sin ^2{\frac{T\Omega _{\rm O}}{2}} \frac{\sin {\phi _{\rm O}}}{\cos {\phi _{\rm O}}} - 2\cos ^2{\frac{T\Omega _{\rm O}}{2}} \frac{\cos {\phi _{\rm O}}}{\sin {\phi _{\rm O}}})], \nonumber \end{eqnarray*}$$ $$\begin{eqnarray*} g_{\mu _{\alpha }f_0} &\approx & \frac{4\pi ^2 R_{\rm ES}f_0T}{2 c\Omega _{\rm O}} \sin {\frac{\Omega _{\rm O} T}{2}} [-\cos {\delta _0}\sin {\alpha _0} \nonumber \\ & +& \cos {\epsilon }\cos {\delta _0}\cos {\alpha _0}], \nonumber \\ g_{\mu _{\alpha }f_1} &\approx & \frac{4\pi ^2 R_{\rm ES}f_0T^2}{4 c\Omega _{\rm O}} \cos {\frac{\Omega _{\rm O} T}{2}} [-\cos {\delta _0}\sin {\alpha _0}(1+\frac{1}{3}\sin {\phi _{\rm O}}) \nonumber \\ &-& \cos {\epsilon }\cos {\delta _0}\cos {\alpha _0} (1+\frac{1}{3}\cos {\phi _{\rm O}})], \nonumber \\ g_{\mu _{\delta }f_0} &\approx & \frac{4\pi ^2 R_{\rm ES}f_0T}{2 c\Omega _{\rm O}} \sin {\frac{\Omega _{\rm O} T}{2}} [-\cos {\alpha _0}\sin {\delta _0} \nonumber \\ &-& \cos {\epsilon }\sin {\delta _0}\sin {\alpha _0} + \sin {\epsilon } \cos {\delta _0} ], \nonumber \\ g_{\mu _{\delta }f_1} &\approx & \frac{4\pi ^2 R_{\rm ES}f_0T^2}{4 c\Omega _{\rm O}} \cos {\frac{\Omega _{\rm O} T}{2}} [-\cos {\alpha _0}\sin {\delta _0}(1+\frac{1}{3}\sin {\phi _{\rm O}}) \nonumber \\ &+& (1+\frac{1}{3}\cos {\phi _{\rm O}})(\cos {\epsilon }\sin {\delta _0}\sin {\alpha _0} - \sin {\epsilon }\cos {\delta _0}) ]. \end{eqnarray*}$$(19) It can be noticed that these metric components depend on the sky position of the source and on its frequency, in a very similar way to the sky position metric components. For a search that has to cover all the sky, this would produce difficulties in the template bank construction (as explained in Wette 2014), but, as argued before, the proper motion components only produce a noticeable mismatch for observation times longer than a year and with less than five segments. All-sky and directed searches only allow such long coherent times at the last stages of the follow-up procedure. In these stages the sky position of the source has already been determined with enough accuracy that it can be used as a constant input to the proper motion metric components. In order to validate the previous metric components, we calculate the relative error ϵr between the true and predicted mismatches, in order to study the behaviour of the mismatch: $$\begin{eqnarray*} \epsilon _r = 2\frac{m_0 - m}{m_0 + m}, \end{eqnarray*}$$(20) where m0 is the true mismatch obtained by calculating the |$\mathcal {F}$|-statistic values and m is the mismatch predicted by the phase metric components. The relative error is negative when the predicted mismatch is higher than the true mismatch, meaning that we are overestimating the mismatch. The relative error is shown in Fig. 5 for a different number of coherent integration times (106 signals have been used for each coherent time). Overall, the errors that we obtain have similar values to the errors obtained in other papers where the metric components are estimated, such as Wette (2014). It can be seen that the relative error decreases as the coherent time is increased, due to the approximated metric components being more accurate as the neglected terms are less important. Figure 5. Open in new tabDownload slide Relative error given by equation (20) as a function of the coherent integration time. For each time, the maximum, mean, and minimum points are shown, along with an error bar comprising one standard deviation. Figure 5. Open in new tabDownload slide Relative error given by equation (20) as a function of the coherent integration time. For each time, the maximum, mean, and minimum points are shown, along with an error bar comprising one standard deviation. 4.3 Mismatch for frequency searches Many semicoherent CW searches (such as the two Hough algorithms used in Abbott et al. 2019a) just track the frequency-time pattern of the signal, instead of searching the full signal given by equation (1). For these searches, we can estimate if proper motion will produce non-negligible mismatch by calculating the difference between the true and searched tracks, a method that has been previously used to estimate the mismatch produced by higher order spin-down terms or by neglected eccentricity in binary systems (Krishnan et al. 2004; Covas & Sintes 2019). The frequency-time pattern (assuming s = 1) is given by $$\begin{eqnarray*} f_{\rm P}(t) &=& \frac{1}{2\pi } \frac{{\rm d}\phi (t)}{{\rm d}t} \nonumber \\ &=& f^{\prime }_0 + f^{\prime }_0\frac{\vec{v}(t)\cdot \hat{n}(t) + \vec{r}(t)\cdot \hat{v}_{\rm s}(t)}{c} + f^{\prime }_1 t \nonumber \\ &\approx & f^{\prime }_0 + f^{\prime }_0\frac{\vec{v}(t)\cdot \hat{n}_0}{c} \nonumber \\ &+& f^{\prime }_0\left(\frac{\vec{v}(t)\cdot \hat{\dot{n}}(t_{\rm r}) (t-t_{\rm r})}{c} + \frac{\vec{r}(t)\cdot \hat{\dot{n}}(t_{\rm r})}{c} \right) + f^{\prime }_1 t, \end{eqnarray*}$$(21) where |$\hat{n}_0 \equiv \hat{n}(t_{\rm r})$|⁠, |$\vec{v}(t)$| is the velocity of the detector, and |$\hat{v}_{\rm s}(t)\approx \hat{\dot{n}}(t_{\rm r})$| is the velocity of the source. When the proper motion parameters are not searched the frequency-time pattern is given by $$\begin{eqnarray*} f(t) = f^{\prime }_0 + f^{\prime }_0\frac{\vec{v}(t)\cdot \hat{n}_0}{c} + f^{\prime }_1 t. \end{eqnarray*}$$(22) The difference between the two frequency-time patterns is $$\begin{eqnarray*} |f_{\rm P}(t)- f(t)| &= f^{\prime }_0 \left(\frac{\vec{v}(t)\cdot \hat{\dot{n}}(t_{\rm r}) (t-t_{\rm r})}{c} + \frac{\vec{r}(t)\cdot \hat{\dot{n}}(t_{\rm r})}{c} \right). \end{eqnarray*}$$(23) A quick estimate shows that for signals with f0 = 1000 Hz and an observation time of 1 yr, the change in frequency will be smaller than 10−7 Hz for a total proper motion of 10−14 rad s−1. The coherent time of these searches is usually less than 7200 s, which implies a frequency resolution of df0 ∼ 10−4 Hz. This is shown in Fig. 6, where four different traces are shown, for two different observation times and total proper motion values. Figure 6. Open in new tabDownload slide Maximum difference between the frequency-time patterns given by equation (23). Two different observing times and total proper motion values are shown. Figure 6. Open in new tabDownload slide Maximum difference between the frequency-time patterns given by equation (23). Two different observing times and total proper motion values are shown. This shows that most semicoherent searches are not able to detect the changes produced by proper motion since the calculated frequency evolution does not deviate by more than a frequency bin. This estimate is in agreement with the results shown in Fig. 4, where it can be seen that searches with a large number of segments (as is the case for these frequency tracking methods with short coherent times) do not have a mismatch higher than 10−3. 5 PROPER MOTION PARAMETER ESTIMATION In the previous section we have shown the mismatch present when the proper motion of the neutron star is assumed to be zero. In this section we study the accuracy that can be achieved when these parameters are searched. A lower bound on the best achievable accuracy on the estimation of parameters can be obtained with the Fisher information matrix. As discussed in Prix (2007), the Fisher information matrix is related to the mismatch metric: $$\begin{eqnarray*} \Gamma _{ij} = \rho ^2 g_{ij}. \end{eqnarray*}$$(24) The lower bound on the standard deviation is $$\begin{eqnarray*} \sigma _i \ge \sqrt{\Gamma ^{-1}_{ii}} = \frac{1}{\rho } \sqrt{g^{-1}_{ii}}. \end{eqnarray*}$$(25) From this expression it is clear that the accuracy depends on the inverse of the S/N and the inverse of the metric. For this reason, the accuracy will depend on the frequency, the coherent time, and the S/N of the signal, and it is independent of the absolute value of the proper motion, i.e. we do not get better resolution for higher proper motions. On the other hand, the dependence of the accuracy on the number of detectors is present through the S/N of the signal: Although the metric components are the same when more detectors are added (if the noise floors are the same), the S/N of the signal increments and so does the accuracy. From the expression, it is also clear that a better accuracy will be achieved when only the proper motion parameters are unknown, since a search over other parameters will decrease the accuracy due to covariance. This means that in general the estimation of the proper motion parameters will be better for directed searches as compared to all-sky searches. In order to empirically study the accuracy, we have modified the MCMC routine part of the pyfstat repository (Ashton & Prix 2018; Ashton et al. 2020). This software allows to do an |$\mathcal {F}$|-statistic search using a parallel tempered MCMC follow-up (Foreman-Mackey et al. 2013; Vousden, Farr & Mandel 2015). We add the two proper motion parameters, and we do a coherent search with different coherent times (without data gaps) by searching an interval around the true parameters. We simulate an all-sky search where 4 + 2 parameters are searched, with two detectors of stationary Gaussian noise, and we add signals with a range of S/Ns and with different Doppler parameters. First, we check that we recover the injections with the correct parameters. The pp-plot in Fig. 7 shows that we recover the proper motions of the injections within the expected credible regions (the same happens for all the other parameters describing the signal). This is calculated by taking the output chains of the MCMC procedure and calculating the credible region of each parameter. When the proper motion parameters are not searched, the obtained points do not follow the expected straight line (as expected from the biases found in subSection 4.1), but when these parameters are included this figure shows that we recover the expected behaviour. Figure 7. Open in new tabDownload slide Probability–probability plot showing the quantiles of the posterior distribution and the fraction of recovered signals whose true parameters are located inside the quantile. Figure 7. Open in new tabDownload slide Probability–probability plot showing the quantiles of the posterior distribution and the fraction of recovered signals whose true parameters are located inside the quantile. Secondly, we calculate the uncertainty of the recovered posterior distribution of the proper motion parameters. Fig. 8 shows the relative uncertainty as a function of the proper motion and the squared S/N of the signal, for a 2-yr coherent search. As previously done, we inject signals at a range of S/Ns, with isotropic orientations and sky distribution. It can be seen that, as expected, higher S/Ns produce lower uncertainties, but other parameters not shown in this plot (such as the sky position or the frequency) also contribute to the vertical spread of similar proper motion values. The figure shows that for proper motions smaller than 10−15 rad s−1 a relative uncertainty smaller than 1 cannot be obtained (for a 2-yr search with two detectors). Figure 8. Open in new tabDownload slide Relative uncertainty (1σ divided by the true proper motion value) of the posterior distributions recovered with the MCMC follow-up scheme as a function of the proper motion value, for an all-sky coherent search of 2 yr with two detectors. The colour of each point (circles for the Right Ascension and squares for the Declination) shows the squared S/N of the signal. Figure 8. Open in new tabDownload slide Relative uncertainty (1σ divided by the true proper motion value) of the posterior distributions recovered with the MCMC follow-up scheme as a function of the proper motion value, for an all-sky coherent search of 2 yr with two detectors. The colour of each point (circles for the Right Ascension and squares for the Declination) shows the squared S/N of the signal. We have repeated the same injections for a 1-yr search, and the theoretical increase of the uncertainty region following T3/2 has been verified. Furthermore, we have also done a follow-up search where all the parameters except the proper motion ones were fixed to the true values, and the mean reduction of the relative uncertainty is around a factor of 2, meaning that we can get twice (on average) better uncertainty if the other parameters are known exactly and not searched. This could happen if, for example, a source is also found with an electromagnetic search. From these results, we remark that systematic calibration errors present in the current gravitational-wave data (Sun et al. 2020), which are smaller than |$5{{\ \rm per\ cent}}$|⁠, are narrower than our obtained relative uncertainties, at least for observation times less than a couple of years. The accuracy on the proper motion parameters could be improved for some sky positions if ecliptical coordinates were used instead of equatorial coordinates, as, for example, discussed in Matthews et al. (2016). The metric components in ecliptical coordinates are obtained in the Appendix, and it can be seen that the accuracy that can be obtained at the same sky position is different for both coordinate systems. 6 CONCLUSIONS In this pape,r we have studied the effects that proper motion of neutron stars produce on searches for continuous gravitational waves. All past searches have assumed the effect of proper motion to be negligible, but as we have seen this might be dangerous for coherent times longer than a year at frequencies higher than ∼1000 Hz. Our results indicate that at these regions of parameter space, follow-up efforts should include these two extra parameters in the analysis because otherwise a real signal could be missed. Since only a small number of outliers reach the final follow-up stage (where the number of segments is small) of searches for unknown neutron stars, this strategy should not highly increment the total computational cost of a search. This will be important for the upcoming observing runs that are planned to be longer than a year, or for when data from different observing runs is combined. Besides, the danger of missing a signal, we have also seen that even if this is not the case, the estimated parameters of the signal are biased. If these parameters are not included in a search, the reported uncertainties on the estimated parameters should be bigger in order to accomodate these systematic errors. Finally, we have shown the relative uncertainty that can be achieved with a CW search of 2 yr and two detectors, and how it depends on the S/N of the signal. Relative uncertainties smaller than 1 are only possible for proper motion values higher than 10−15 rad s−1 (for observing times less than a couple of years). A higher number of detectors would improve the relative uncertainty even more, due to the increase of the measured S/N. In this paper, we have not discussed the feasibility of measuring radial motion. As briefly mentioned in the signal model section, the effect of radial motion on the phase of the signal is smaller than the effect of transverse motion. Since we have seen that detecting transverse motion requires long observation times, measuring radial motion with CWs will be even more difficult and require observing times of many years. We have assumed no timing noise or spin-wandering, which if present might bias the estimation of parameters in a similar way as the presence of glitches or proper motion. A detailed study of the size of these biases might be useful to uncover if these biases are bigger than the statistical uncertainties associated with our measurements. Another limitation of this study is that we have studied the bias introduced by the dismissal of proper motion in searches for isolated neutron stars. Biases introduced in searches for neutron stars in binary systems may be different, and the parameters describing the binary orbit (both Keplearian and post-Keplerian) might also be affected, as, for example, discussed in Splaver et al. (2005). ACKNOWLEDGEMENTS The authors want to thank Karl Wette, Reinhard Prix, and David Keitel for many helpful discussions. The authors also want to thank the anonymous referee for helping to improve the quality of this paper. We acknowledge the support of the Spanish Agencia Estatal de Investigación and Ministerio de Ciencia, Innovación y Universidades grants FPA2016-76821-P, FPA2017-90687-REDC, FPA2017-90566-REDC, FPA2015-69815-REDT, FPA2015-68783-REDT, the Vicepresidencia i Conselleria d’Innovació, Recerca i Turisme del Govern de les Illes Balears and the Fons Social Europeu 2014-2020 de les Illes Balears, the European Union FEDER funds, and the EU COST actions CA16104, CA16214, and CA17137. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459. This paper has LIGO document number P2000238. DATA AVAILABILITY The data underlying this paper were either downloaded from the ATNF Pulsar Database (Section 3) or generated by using the LALSuite code repository maintained by the LIGO Scientific Collaboration (LIGO Scientific Collaboration 2018) (Sections 4 and 5). More specifically, the codes lalapps_ComputeFstatistic_v2 and UniversalDopplerMetric were used for the results obtained in Section 4, while the software from Ashton & Prix (2018) was used in Section 5, both for the generation of simulated signals and for the subsequent MCMC search. Due to the high number of simulated signals and the large time span of the different generated data sets, the simulated signals cannot be saved to disc and are eliminated after their analysis, although the parameters of the different simulated signals are saved in a text file and are available upon request. REFERENCES Aasi J. et al. , 2015 , Class. Quantum Gravity , 32 , 7 Crossref Search ADS Abbott B. P. et al. , 2017 , Phys. Rev. D , 96 , 122004 10.1103/PhysRevD.96.122004 Crossref Search ADS Abbott B. P. et al. , 2018 , Living Rev. Relativ. , 21 , 3 10.1007/s41114-018-0012-9 Crossref Search ADS PubMed Abbott B. P. et al. , 2019a , Phys. Rev. D , 100 , 024004 10.1103/PhysRevD.100.024004 Crossref Search ADS Abbott B. P. et al. , 2019b , ApJ , 879 , 1 10.3847/1538-4357/ab20cb Crossref Search ADS Acernese F. et al. , 2014 , Class. Quantum Gravity , 32 , 2 Akutsu T. et al. , 2019 , Nat. Astron. , 3 , 35 10.1038/s41550-018-0658-y Crossref Search ADS Ashton G. , Prix R., 2018 , Phys. Rev. D , 97 , 103020 10.1103/PhysRevD.97.103020 Crossref Search ADS Ashton G. , Jones D. I., Prix R., 2015 , Phys. Rev. D , 91 , 062009 10.1103/PhysRevD.91.062009 Crossref Search ADS Ashton G. , Prix R., Jones D. I., 2017 , Phys. Rev. D , 96 , 063004 10.1103/PhysRevD.96.063004 Crossref Search ADS Ashton G. , , Keitel D., Prix R., Tenorio R., 2020 , PyFstat/PyFstat: v1.9.0 : 10.5281/zenodo.4244503 Crossref Chatterjee S. et al. , 2005 , ApJ , 630 , 1 Crossref Search ADS Cordes J. M. , 1987 , in Helfand D., Huang J.-H.eds, The Origin and Evolution of Neutron Stars , Proc. IAU Symp. 125, Kluwer , Dordrecht , p. 35 Covas P. B. , Sintes A. M., 2019 , Phys. Rev. D , 99 , 124019 10.1103/PhysRevD.99.124019 Crossref Search ADS Covas P. B. , Sintes A. M., 2020 , Phys. Rev. Lett. , 124 , 191102 10.1103/PhysRevLett.124.191102 Crossref Search ADS PubMed Cui X. H. , Wang H. G., Xu R. X., Qiao G. J., 2007 , A&A , 472 , 1 10.1051/0004-6361:20077200 Crossref Search ADS Edwards R. T. , Hobbs G. B., Manchester R. N., 2006 , MNRAS , 372 , 4 Crossref Search ADS Foreman-Mackey D. , Hogg D. W., Lang D., Goodman J., 2013 , PASP , 125 , 925 Crossref Search ADS Harrison E. R. , Tademaru E., 1975 , ApJ , 201 , 447 10.1086/153907 Crossref Search ADS Hobbs G. , Lorimer D. R., Lyne A. G., Kramer M., 2005 , MNRAS , 360 , 3 Crossref Search ADS Hui C. Y. , Becker W., 2006 , A&A , 454 , 2 Crossref Search ADS Jaranowski P. , Królak A., 1999 , Phys. Rev. D , 59 , 063003 10.1103/PhysRevD.59.063003 Crossref Search ADS Jaranowski P. , Królak A., Schutz B. F., 1998 , Phys. Rev. D , 58 , 063001 10.1103/PhysRevD.58.063001 Crossref Search ADS Kaplan D. L. , Chatterjee S., Gaensler B. M., Anderson J., 2008 , ApJ , 677 , 2 Crossref Search ADS Krishnan B. , Sintes A. M., Papa M. A., Schutz B. F., Frasca S., Palomba C., 2004 , Phys. Rev. D , 70 , 082001 10.1103/PhysRevD.70.082001 Crossref Search ADS Lai D. , Goldreich P., 2000 , ApJ , 535 , 1 Crossref Search ADS Lai D. , Qian Y. Z., 1998 , ApJ , 505 , 2 Crossref Search ADS Lai D. , Chernoff D. F., Cordes J. M., 2001 , ApJ , 549 , 2 Crossref Search ADS LIGO Scientific Collaboration , 2018 , LIGO Algorithm Library – LALSuite, Free Software (GPL) . 10.7935/GT1W-FZ16 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref Manchester R. N. , Hobbs G. B., Teoh A., Hobbs M., 2005 , AJ , 129 , 4 Crossref Search ADS Matthews A. M. et al. , 2016 , ApJ , 818 , 1 Crossref Search ADS Mukherjee A. , Messenger C., Riles K., 2018 , Phys. Rev. D , 97 , 043016 10.1103/PhysRevD.97.043016 Crossref Search ADS Noutsos A. , Kramer M., Carr P., Johnston S., 2012 , MNRAS , 423 , 3 Crossref Search ADS Pitkin M. , 2018 , J. Open Source Softw. , 3 , 22 Crossref Search ADS Pletsch H. J. , 2011 , Phys. Rev. D , 83 , 122003 10.1103/PhysRevD.83.122003 Crossref Search ADS Prix R. , 2007 , Phys. Rev. D , 75 , 023004 10.1103/PhysRevD.75.023004 Crossref Search ADS Reardon D. J. , Coles W. A., Hobbs G., Ord S., Kerr M., Bailes M., Bhat N. D. R., Krishnan V. V., 2019 , MNRAS , 485 , 3 Crossref Search ADS Sieniawska M. , Bejger M., 2019 , Universe , 5 , 217 10.3390/universe5110217 Crossref Search ADS Splaver E. M. , Nice D. J., Stairs I. H., Lommen A. N., Backer D. C., 2005 , ApJ , 620 , 1 Crossref Search ADS Sun L. et al. , 2020 , Class. Quantum Gravity , 37 , 45 Vousden W. , Farr W. M., Mandel I., 2015 , MNRAS , 455 , 2 Wette K. , 2014 , Phys. Rev. D , 90 , 122010 10.1103/PhysRevD.90.122010 Crossref Search ADS Wette K. , 2016 , Phys. Rev. D , 94 , 122002 10.1103/PhysRevD.94.122002 Crossref Search ADS APPENDIX A: EXPLICIT DERIVATION OF METRIC COMPONENTS In this Appendix, we explicitly derive the |$g_{\mu _{\alpha } \mu _{\alpha }}$| metric component shown in Section 4 as an example. The other metric components can be obtained in a similar way as the one showed below. The phase in the detector frame is given by (up to first-order in frequency derivatives, and assuming that |$f_k \approx f^{\prime }_k$|⁠) equation (3): $$\begin{eqnarray*} \phi (t) &= \phi _0 + 2\pi f_0 \left[ (t-t_{\rm r}) + \frac{\vec{r}(t)\cdot \hat{n}(t)}{c}\right] + \pi f_1 (t-t_{\rm r})^2. \end{eqnarray*}$$(A1) The first step consists on calculating the phase derivative with respect to the parameter: $$\begin{eqnarray*} \frac{\partial \phi }{\partial \mu _{\alpha }} &=& \frac{2\pi f_0}{c} \vec{r}(t)\cdot \frac{\partial \hat{n}(t)}{\partial \mu _{\alpha }} \nonumber \\ \frac{\partial \hat{n}(t)}{\partial \mu _{\alpha }} &=& (t-t_{\rm r})[-\sin {\alpha _0}\cos {\delta _0}, \cos {\alpha _0}\cos {\delta _0}, 0], \end{eqnarray*}$$(A2) where, for the last equation, we have used the Taylor-expanded sky position given by equation (6). The relative importance of the different terms can be estimated with the following order of magnitude calculation (separating |$\vec{r}(t)$| in its two contributions): $$\begin{eqnarray*} \frac{\partial \phi }{\partial \mu _{\alpha }} \sim &\frac{f_0 T R_{\rm ES}}{c} + \frac{f_0 T R_{\rm E}}{c}, \end{eqnarray*}$$(A3) which shows that we can approximate the derivative by just keeping the terms dependent on RES, thus neglecting the rotation of the Earth. The last step consists on calculating the time integrals: $$\begin{eqnarray*} &&\langle \frac{\partial \phi }{\partial \mu _{\alpha }} \frac{\partial \phi }{\partial \mu _{\alpha }} \rangle = \frac{1}{T} \int _{t_0}^{t_0+T} \left(\frac{2\pi f_0}{c} \vec{r}_{\rm O}(t)\cdot \frac{\partial \hat{n}(t)}{\partial \mu _{\alpha }}\right)^2 {\rm d}t \nonumber \\ &&= \frac{4\pi ^2 f_0^2 R^2_{\rm ES}}{T c^2} \int _{t_0}^{t_0+T} (t-t_{\rm r})^2 (-\sin {\alpha _0}\cos {\delta _0}\cos {[\phi _O + \Omega _O (t-t_{\rm r})]} \nonumber \\ &&+ \cos {\alpha _0}\cos {\delta _0}\cos {\epsilon }\sin {[\phi _O + \Omega _O (t-t_{\rm r})]})^2 {\rm d}t \nonumber \\ &&\approx \frac{4\pi ^2 f_0^2 R^2_{\rm ES}}{T c^2}\left[ \right. (\sin ^2{\alpha _0}\cos ^2{\delta _0}+\cos ^2{\alpha _0}\cos ^2{\delta _0}\cos ^2{\epsilon }) \nonumber \\ &&\left(\frac{(T+t_0)^3}{6} \right. - \frac{t_{\rm r}(T+t_0)^2}{2} + \frac{t_{\rm r}^2(T+t_0)}{2} - \frac{t_{\rm r}^2 t_0}{2} + \frac{t_{\rm r} t_0^2}{2} - \frac{t_0^3}{6} \nonumber \\ &&+ \mathcal {O} \left. \left(\frac{T^2}{\Omega _{\rm O}}\right) + \mathcal {O}\left(\frac{T}{\Omega _{\rm O}^2}\right) + \mathcal {O}\left(\frac{1}{\Omega _{\rm O}^3}\right)\right) \nonumber \\ &&- (\sin {\alpha _0}\cos ^2{\delta _0}\cos {\alpha _0}\cos {\epsilon })\left.\left(\mathcal {O}\left(\frac{T^2}{\Omega _{\rm O}}\right) + \mathcal {O}\left(\frac{T}{\Omega _{\rm O}^2}\right) + \mathcal {O}\left(\frac{1}{\Omega _{\rm O}^3} \right) \right) \right], \nonumber \\ \end{eqnarray*}$$(A4) $$\begin{eqnarray*} \langle \frac{\partial \phi }{\partial \mu _{\alpha }} \rangle \langle \frac{\partial \phi }{\partial \mu _{\alpha }} \rangle = \frac{1}{T^2} \left(\int _{t_0}^{t_0+T} \frac{2\pi f_0}{c} \vec{r}_{\rm O}(t)\cdot \frac{\partial \hat{n}(t)}{\partial \mu _{\alpha }} {\rm d}t \right)^2 \propto \mathcal {O}(\frac{T}{\Omega _{\rm O}}). \nonumber \\ \end{eqnarray*}$$(A5) It can be seen that for integration times longer than a year, terms such as T2/ΩO are smaller than T3. Now, we fix the reference time to two different values. For tr = t0 + T/2, $$\begin{eqnarray*} g_{\mu _{\alpha } \mu _{\alpha }} = \frac{4\pi ^2 f_0^2 R^2_{\rm ES} T^2}{24 c^2} (\sin ^2{\alpha _0}\cos ^2{\delta _0}+\cos ^2{\alpha _0}\cos ^2{\delta _0}\cos ^2{\epsilon }) + \mathcal {O}(T/\Omega _{\rm O}), \nonumber \\ \end{eqnarray*}$$(A6) while for tr = t0 or tr = t0 + T, $$\begin{eqnarray*} g_{\mu _{\alpha } \mu _{\alpha }} = \frac{4\pi ^2 f_0^2 R^2_{\rm ES} T^2}{6c^2} (\sin ^2{\alpha _0}\cos ^2{\delta _0}+\cos ^2{\alpha _0}\cos ^2{\delta _0}\cos ^2{\epsilon }) + \mathcal {O}(T/\Omega _{\rm O}). \nonumber \\ \end{eqnarray*}$$(A7) Reference times selected between the initial and mid-time will produce metric components with values between these two extremes. With this approximation and equation (17), we obtain the metric elements shown in Section 4. In ecliptical coordinates where l is the longitude and b is the latitude, we define the source and Earth positions: $$\begin{eqnarray*} \hat{n}(t) &=& \hat{n}(t_{\rm r}) + \hat{\dot{n}}(t_{\rm r})(t-t_{\rm r}) \nonumber \\ &=& [\cos {l_0}\cos {b_0},\sin {l_0}\cos {b_0},\sin {b_0}] \nonumber \\ &+& (t-t_{\rm r})[-\mu _{l}\sin {l_0}\cos {b_0} - \mu _{b}\cos {l_0}\sin {l_0}, \nonumber \\ &\mu & _{l}\cos {l_0}\cos {b_0} - \mu _{b}\sin {l_0}\sin {b_0}, \mu _{b}\cos {b_0}], \end{eqnarray*}$$(A8) $$\begin{eqnarray*} \vec{r}_O (t) &=& R_{ES}[\cos {\left(\phi _O + \Omega _O (t-t_r)\right)}, \nonumber \\ && \sin {\left(\phi _O + \Omega _O (t-t_r)\right)}, 0] . \end{eqnarray*}$$(A9) When the time integrals are done in these coordinates, the results are (for tr = t0 + T/2) $$\begin{eqnarray*} g_{\mu _{l}\mu _{l}} &=& \frac{4\pi ^2 R_{\rm ES}^2 f_0^2 T^2}{24c^2} \cos ^2{b} + \mathcal {O}(T/\Omega _{\rm O}), \nonumber \\ g_{\mu _{b}\mu _{b}} &=& \frac{4\pi ^2 R_{\rm ES}^2 f_0^2 T^2}{24c^2} \sin ^2{b} + \mathcal {O}(T/\Omega _{\rm O}), \nonumber \\ g_{\mu _{l}\mu _{b}} &\propto & \mathcal {O}(T/\Omega _{\rm O}). \end{eqnarray*}$$(A10) It can be seen that in these coordinates, the covariant component |$g_{\mu _{l}\mu _{b}}$| is reduced as compared to the equatorial coordinates |$g_{\mu _{\alpha }\mu _{\delta }}$| component. © 2020 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Effects of proper motion of neutron stars on continuous gravitational-wave searches JO - Monthly Notices of the Royal Astronomical Society DO - 10.1093/mnras/staa3624 DA - 2020-12-11 UR - https://www.deepdyve.com/lp/oxford-university-press/effects-of-proper-motion-of-neutron-stars-on-continuous-gravitational-m2y0tr2OP0 SP - 5167 EP - 5176 VL - 500 IS - 4 DP - DeepDyve ER -