Constraints on ultra-low-frequency gravitational waves with statistics of pulsar spin-down rates
Hiroki Kumamoto, Yuya Imasato, Naoyuki Yonemaru, Keitaro Takahashi, Sachiko Kuroyanagi
aa r X i v : . [ a s t r o - ph . H E ] S e p MNRAS , 1–7 (2019) Preprint 1 October 2019 Compiled using MNRAS L A TEX style file v3.0
Constraints on ultra-low-frequency gravitational waveswith statistics of pulsar spin-down rates
Hiroki Kumamoto , ⋆ , Yuya Imasato , Naoyuki Yonemaru , , Sachiko Kuroyanagi ,and Keitaro Takahashi , . Kumamoto University, Graduate School of Science and Technology, Japan CSIRO Astronomy and Space Science, PO Box 76, Epping NSW 1710, Australia Nagoya University, Graduate School of Science, Japan International Research Organization for Advanced Science and Technology, Japan
Accepted 2019 August 18. Received 2019 August 18; in original form 2019 March 4
ABSTRACT
We probe ultra-low-frequency gravitational waves (GWs) with statistics of spin-downrates of milli-second pulsars (thereafter MSPs) by a method proposed in our prevouswork (Yonemaru et al. 2016). The considered frequency range is − Hz . f GW . − Hz . The effect of such low-frequency GWs appears as a bias to spin-down rateswhich has a quadrupole pattern in the sky. We use the skewness of the spin-down ratedistribution and the number of MSPs with negative spin-down rates to search for thebias induced by GWs. Applying this method to 149 MSPs selected from the ATNFpulsar catalog, we derive upper bounds on the time derivative of the GW amplitudesof Û h < . × − sec − and Û h < . × − sec − in the directions of the GalacticCenter and M87, respectively. Approximating the GW amplitude as Û h ∼ π f GW h , thebounds translate into h < × − and h < × − , respectively, for f GW = /( ) .Finally, we give the implications to possible super-massive black hole binaries at thesesites. Key words: gravitational waves – methods: data analysis – methods: statistical –pulsars: general.
Laser Interferometer Gravitational Wave Observatory(LIGO) has succeeded in detecting gravitational waves(GWs) with frequencies ∼
100 Hz radiated from blackhole binaries with masses of about M ⊙ (Abbott et al.2016). With KAGRA (KAGRA Collaboration 2019) andVIRGO (Acernese et al. 2015), ground-based interferome-ters will pave the way for gravitational-wave astronomy.In fact, this is the first step toward multi-wavelengthgravitational-wave astronomy and lower-frequency GWs areto be probed: space interferometers such as Laser In-terferometer Space Antenna (LISA) (Amaro-Seoane et al.2017) and DECIGO (Kawamura et al. 2011), Pulsar Tim-ing Arrays (PTAs) such as Parkes Pulsar Timing Ar-ray (PPTA) (Manchester et al. 2012), European PulsarTiming Array (EPTA) (Kramer & Champion 2013), NorthAmerican Nanohertz Observatory for Gravitational Waves(NANOGrav) (Jenet et al. 2009) and International PulsarTiming Array (IPTA) (Verbiest et al. 2016), and observa- ⋆ E-mail: [email protected] tions of B-mode polarization of the cosmic microwave back-ground such as Simons Observatory (Ade et al. 2018) andLiteBIRD (Ishino et al. 2016).GWs radiated from possible super-massive black hole(SMBH) binaries in the Galactic Center (GC) will greatlyimprove our understanding of gravity theories, astrophysicsof SMBHs and environment of the GC region. In the latestage of SMBH binary evolution, GWs with frequencies of − Hz . f GW . − Hz are radiated and the range ismain target of PTAs. The frequency range of PTA is de-termined by the observational time span and cadence, andlower frequencies f GW . − Hz are difficult to probe. Infact, the sensitivety is expected to scale as f − toward lowerfrequencies (Blandford et al. 1984; Moore et al. 2015).In Yonemaru et al. (2016), we proposed a new detec-tion method for GWs with ultra-low-frequencies of f GW . − Hz (for other methods, see Bertotti et al. (1983);Kopeikin (1997); Potapov et al. (2003)). The method uti-lizes the fact that the spin-down rate of milli-second pul-sars (MSPs) is biased by such GWs, since both give thesame quadratic time dependence to the time of arrival ofpulses. This effect depends on the relative direction of the © H., Kumamoto. et al.
GW source and a pulsar. Thus, statistics of the spin-downrate distributions can probe such GWs as we describe later.In Yonemaru et al. (2018), by using simulated 3,000 MSPs,which are expected to be discovered by the Square Kilome-tre Array (SKA) survey, we estimated the sensitivity of thismethod in a simple situation, where we assume that MSPsare located uniformly in the sky and the “pulsar term” isneglected. We concluded that GWs with the derivative ofamplitude as small as × − s − could be detected. Then,in Hisano et al. (2019), we considered a more realistic modelof MSP distribution in the Galaxy and took the pulsar terminto account in order to obtain more accurate estimates ofthe sensitivity by extending the analysis of Yonemaru et al.(2018). We found that the sensitivity depends on the direc-tion, polarization and frequency of GWs and becomes worseat low frequencies ( f GW . − Hz ) because of the pulsarterm.This work is the first attempt to apply the abovemethod to real data. We use MSPs selected from the currentATNF pulsar catalog (Manchester et al. 2005) and probeGWs with a frequency of − Hz . f GW . − Hz . InSection 2, we give a brief summary on our method proposedin Yonemaru et al. (2016). Then, after describing our dataset, upper bounds on the derivative of GW amplitude are de-rived in Section 3. In Section 4, we discuss the implication ofthe upper bounds to possible SMBH binaries at the Galac-tic Center and M87. Finally, our results are summarized inSection 5. Let us first describe the detection method of ultra-low-frequency GWs following Yonemaru et al. (2016). Timingresiduals induced by GW are given by (Detweiler 2010), r GW ( t ) = Õ A =+ , × F A ( ˆ Ω , ˆ p ) ∫ t ∆ h A ( t ′ , ˆ Ω , θ ) dt ′ , (1)where we denote the direction of pulsar as ˆ p , the propagationdirection of GW as ˆ Ω and the GW polarization angle as θ . Here, antenna beam pattern F A ( ˆ Ω , ˆ p , θ ) is the geometricfactor written by (Anholm et al. 2009), F A ( ˆ Ω , ˆ p ) =
12 ˆ p i ˆ p j + ˆ Ω · ˆ p e Aij ( ˆ Ω ) , (2)where e Aij ( ˆ Ω ) ( A = + , ×) are the GW polarization tensor givenby e + ij ( ˆ Ω ) = ˆ m i ˆ m j − ˆ n i ˆ n j , (3) e × ij ( ˆ Ω ) = ˆ m i ˆ n j + ˆ n i ˆ m j , (4)with ˆ m and ˆ n being the polarization basis vectors. In Eq.(1), ∆ h A ( t ′ , ˆ Ω , θ ) is the difference of geometric perturbation be-tween the earth and the pulsar. This is given by, ∆ h A ( t ′ , ˆ Ω , θ ) = h A ( t , ˆ Ω , θ ) − h A ( t p , ˆ Ω , θ ) , (5)where t p = t − τ with τ = L / c ( + ˆ Ω · ˆ p ) being the pulsepropagation time from the pulsar at the distance L to theearth.In the following, we will discuss GWs with periods muchlonger than the observational time span, and in this case, the -1x10 -18 -8x10 -19 -6x10 -19 -4x10 -19 -2x10 -19 -19 -19 -19 -19 -18 Figure 1.
Spatial pattern of the bias factor α ( ˆ Ω , ˆ p , θ ) in the skyfor Û h + = − s − . The GW source position is placed at the centerof the celestial sphere in equatorial coordinates. GW amplitude changes linearly with time. At the same time,we use the assumption that the second term (”pulsar term”)is a random noise with zero average, which is reasonablewhen the GW wavelength is much shorter than the typicaldistance to pulsars. Thus, the GW frequency range we con-sider here is − Hz . f GW . − Hz. For such GWs, wecan rewrite Eq.(5) as ∆ h A ( t ′ , ˆ Ω , θ ) = Û h A ( ˆ Ω , θ ) t . (6)Then, substituting Eq.(5) into Eq.(1), we find the timingresidual induced by ultra-low frequency GWs is describedby r GW ( t ) = Õ A =+ , × F A ( ˆ Ω , ˆ p ) Û h A ( ˆ Ω , θ ) t . (7)This time dependence is the same as timing residual inducedby pulsar spin down, which is given by r Û p ( t ) = Û pp t , (8)where p and Û p are the pulse period and spin-down rate,respectively. Therefore, the influence of ultra-low-frequencyGWs is absorbed into the spin-down rate of the pulsar andcannot be identified in the standard analysis of PTA. Onone hand, in the presence of ultra-low-frequency GWs, spin-down rate is biased as, Û p obs p = Û p p + α ( ˆ Ω , ˆ p , θ ) , (9)where Û p obs and Û p are observed and intrinsic spin-down rates,respectively, and the bias factor α ( ˆ Ω , ˆ p ) is given by, α ( ˆ Ω , ˆ p , θ ) = Õ A =+ , × F A ( ˆ Ω , ˆ p ) Û h A ( ˆ Ω , θ ) . (10)The bias factor depends on the relative direction between theGW propagation and each pulsar and the spatial pattern isplotted in Fig. 1. Here, Û h + ( Ω , θ ) and Û h × ( Ω , θ ) are depend onGW polarization θ and given by Û h + ( ˆ Ω , θ ) = Û h ( ˆ Ω ) cos 2 θ, (11) Û h × ( ˆ Ω , θ ) = Û h ( ˆ Ω ) sin 2 θ. (12)In Yonemaru et al. (2016), we proposed a method to MNRAS , 1–7 (2019) onstrains on ultra-low-frequency GWs N u m b er log (p ·/ p) Figure 2.
Histogram of spin-down rates of 149 MSPs. utilize the spatial pattern of α ( ˆ Ω , ˆ p , θ ) to probe ultra-low-frequency GWs. First, by assuming GW source position andpolarization, we divide pulsars into two groups according tothe sign of the bias factor, depending on the location ofeach pulsar. Although GW signals cannot be extracted fromindividual pulsars, since spin-down rates are biased to pos-itive and negative values in the two groups respectively, itis possible to detect GWs by measuring the systematic dif-ference in the spin-down rate distribution between the twogroups. We use the skewness of the spin-down rate distribu-tion to characterize the bias induced by GWs, and convertthe difference in the skewness of the two groups to the valueof Û h . Below, we apply this method to real data and deriveconstraints on Û h by analyzing the skewness difference. In ad-dition, if the amplitude of GW is too strong, some of pulsarsin the region with negative α ( ˆ Ω , ˆ p , θ ) get negative spin-downrate. We derive constraints on the GW amplitude by usingthe number of pulsars with negative spin-down rates in thereal pulsar observation. We use data of observed MSPs in ATNF pulsar catalogver. 1.59. The data set includes 181 MSPs with the mea-sured periods shorter than 30 msec and the time derivatives.Among 181 MSPs, we exclude 30 MSPs in globular clusterssince they would be biased significantly by the gravitationalpotential and complicated dynamics inside the cluster. Inaddition, two MSPs are removed as outliers: one with anegative spin-down rate ( Û p obs / p = − − . [ sec − ] , J1801-3210) and one with an exceptionally large spin-down rate ( Û p obs / p = − . [ sec − ] , J0537-6910). Thus, 149 MSPs areused for our analysis below.Fig. 2 shows the histogram of Û p obs / p of 149 MSPs. Mean,standard deviation, skewness and kurtosis of the distributionare − . , . , . and . , respectively, and the deviationfrom Gaussian distribution was shown to be statistically sig-nificant by the Jarque-Bera test (Yonemaru et al. 2018). Given the propagation direction and polarization of GWs,the sky is divided into two areas according to the sign of -0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 0 20 40 60 80 100 120 140 160 180 Sk e w n e ss d i ff ere n ce Polarization
GCM87
Figure 3.
Skewness difference as a function of GW polarizationangle for the cases where the location of the GW source is assumedin the direction of the Galactic Center (red) and M87 (blue). the bias factor α ( ˆ Ω , ˆ p , θ ) . Then MSPs are classified into twogroups and skewness of the Û p obs / p distribution is calcu-lated for each group. The skewness in the positive (negative) α ( ˆ Ω , ˆ p , θ ) region is given by S α + (−) = σ + (−) N + (−) N + (−) Õ i (cid:18) log (cid:18) Û p obs p (cid:19) i − µ + (−) (cid:19) . (13)where i = , · · · , N + (−) is the number of MSP in the positive(negative) α ( ˆ Ω , ˆ p , θ ) region, and µ + (−) and σ + (−) are the meanvalue and variance of the log Û p / p distribution, respectively, µ + (−) = N + (−) N + (−) Õ i log (cid:18) Û p obs p (cid:19) i , (14) σ + (−) = N + (−) N + (−) Õ i (cid:18) log (cid:18) Û p obs p (cid:19) i − µ + (−) (cid:19) . (15)Then the skewness difference is given by ∆ S = S α + − S α − . (16)Fig. 3 shows the skewness difference as a function of GWpolarization angle for the cases where we assume that theGW source is located in the direction of the Galactic Centerand M87. It should be noted that polarization angles of and degree correspond to the same polarization, but thesign of α is inverted so that the sign of skewness difference isalso inverted. Maximum values for the case of the GalacticCenter and M87 are 0.672 and 0.676, respectively.In the same way, we calculate the skewness differencefor all directions of GW source and all angles of GW polar-ization. In Fig. 4, at each point of the sky where GW sourceis assumed to be, we have searched for the largest skewnessdifference by changing the GW polarization. The skewnessdifference is mostly smaller than unity and the largest valuein the whole sky is . . The two red regions are antipodes.In order to study the statistical significance of this value,we perform a series of simulations. First, we make mockdata of 149 MSPs in absence of GWs. MSPs are locatedat the same position as that of the real data. The valueof Û p obs / p is randomly allocated to each MSP according toa generalized Gaussian distribution with the same values of MNRAS , 1–7 (2019)
H., Kumamoto. et al.
Figure 4.
Skewness difference distribution in the sky. At eachpoint of the sky, where the GW source is assumed to be, wecompute the skewness difference by varying GW polarization an-gle and the colors represent the maximum skewness difference ateach point. Black points represent the position of MSPs. Red ”+”and blue ” × ” represent the position of the Galactic Center andM87, respectively. mean, standard deviation and skewness as the ones obtainedfrom the catalog. The generalized Gaussian distribution isgiven by f ( x ) = φ ( y ) α − κ ( x − ξ ) , (17)where φ ( y ) is the standard Gaussian distribution and y isgiven by y = ( − κ log [ − κ ( x − ξ ) α ] ( κ , ) κ ( x − ξ ) α ( κ = ) , (18)Here, ξ , α and κ are the location, scale and shape parame-ters, respectively, and the mean µ , standard deviation σ andskewness S are expressed by these parameters. µ = ξ − ακ (cid:16) e κ / − (cid:17) , (19) σ = s α κ (cid:16) e κ − (cid:17) , (20) S = e κ − e κ − ( e κ − ) / sgn ( κ ) . (21)For each realization of mock data, we obtain the skewnessdifference in the same way as above and search for the max-imum varying the position and polarization.Fig. 5 shows the probability distribution function of themaximum skewness difference in the sky obtained through10,000 realizations of mock MSP data without GW injection.We find that the distribution extends from . to . and,as a result, the obtained value . is fairly consistent withthe statistical fluctuations without GWs. In the previous subsection, we have shown that the cur-rent pulsar data is consistent with the non-existence of GWswithin the statistical error. In this subsection, we derive up-per bounds on the derivative of the GW amplitude Û h . Due P r o b a b ili t y Maximum skewness difference
Figure 5.
Simulated probability distribution function of the max-imum skewness difference in the sky, which is obtained through10,000 realizations of mock MSP data. The vertical line shows thevalue . obtained from the real MSP data. to the limited computational power, we focus on two as-trophysically important directions of Galactic Center andM87 where the existence of supermassive black hole binarieshas been suggested (Yu & Tremaine 2003; Oka et al. 2016;Yonemaru et al. 2016).As we saw in the previous subsection, the maximumskewness differences in the direction of the Galactic Cen-ter and M87 are 0.672 (at polarization angle of 25 deg) and0.676 (at polarization angle of 108 deg), respectively. In thepresence of GWs with large enough value of Û h , the probabil-ity of obtaining such small values is low. We place the upperbound on Û h by using the threshold where the probability ofhaving skewness difference less than 0.672 (0.676) is forthe direction of the Galactic Center (M87).In order to evaluate the upper bounds, we make mockdata of 149 MSPs in the same way as the previous subsec-tion and inject the bias due to GWs to the simulated in-trinsic spin-down rates Û p / p . The polarization angle of GWsis set to be the same as the above: 25 deg for the GalacticCenter and 108 deg for M87. In this way, a simulated valueof Û p obs / p is allocated to each MSP and then we computethe skewness difference for the data set. It should be notedthat MSPs with very small values of intrinsic spin-down rate( ∼ − sec ) can have negative values of observed spin-downrate if they are located at an area with negative bias. Theyare removed from the computation of skewness differenceand the total number of used MSPs is slightly smaller than149. Fig. 6 shows the probability distribution function ofskewness difference in the directions of the Galactic Cen-ter and M87 estimated from 10,000 realizations of simula-tions with and without GWs. Here, we fix GW polarizationangles at
25 deg (GC) and
108 deg (M87), which give thelargest skewness difference in Fig. 3. We see that, as the valueof Û h increases, the probability of having a larger value ofskewness difference becomes higher. Comparing them withthe observed values (0.672 and 0.676), we obtain an upperbound of Û h < . × − sec − for the Galactic Center and Û h < . × − sec − for M87. The implication of the upperbounds will be discussed in Section. 4. MNRAS , 1–7 (2019) onstrains on ultra-low-frequency GWs P r o b a b ili t y Skewness difference no GWGW 1e-18GW 6.2e-18 0 0.5 1 1.5 2 0 1 2 3 4 5 P r o b a b ili t y Skewness difference no GWGW 1e-18GW 8.1e-18
Figure 6.
Simulated probability distribution function of skew-ness difference in the directions of the Galactic Center (top) andM87 (bottom) estimated from 10,000 realizations of simulationswithout GWs (red) and with GWs (green and blue). The value of Û h is − sec − (green) and . × − sec − (blue) for the Galac-tic Center, and − sec − (green) and . × − sec − (blue) forM87. As we mentioned in the previous subsection, when MSPswith very small intrinsic spin-down rates are biased nega-tively, they can have negative values of observed spin-downrate. The number of such MSPs will increase for strongerGWs. Therefore, the number of MSPs with negative ob-served spin-down rates could be used for another measureto probe ultra-low-frequency GWs. In the current data set,there is only one MSP with a negative Û p obs / p except onesin globular clusters. Here, we set upper bounds on Û h usingthe threshold where the probability of having two or moreMSPs have negative Û p obs / p is .Fig. 7 shows upper bounds on Û h as a function of GWpolarization angle for the Galactic Center and M87. The up-per bounds are of order − − − . sec − and comparableto those from skewness difference obtained in the previoussubsection. GWs from the Galactic Center are slightly wellconstrained than those from M87 and the dependence on thepolarization angle is very weak. -17.4-17.35-17.3-17.25-17.2-17.15-17.1-17.05-17 0 20 40 60 80 100 120 140 160 180 h · Polarization
GCM87
Figure 7.
Upper bounds on Û h from the number of MSPs with Û p obs / p as a function of GW polarization angle for the GalacticCenter (red) and M87 (green). In the previous section, we have obtained upper bounds onthe time derivative of GW amplitudes, rather than the am-plitudes themselves. Using an approximation Û h ∼ π f GW h where f GW is the frequency of GW, which is reasonable formost of the periods, our constraints for the Galactic Centerand M87 can be approximated as, h GC . × − (cid:18) / f GW (cid:19) , (22) h M87 . × − (cid:18) / f GW (cid:19) , (23)respectively. On the other hand, the recent PTA analysesput upper bounds of h GW ∼ . × − at the frequency of8 nHz and the bound scales as f − (Aggarwal et al 2018;Babak et al. 2015). Thus, our constraints are comparable tothose of standard PTAs at f GW ∼ / ∼ − Hz and are better at even lower frequencies. Although our con-straints are weaker for f GW & − Hz , they are still valu-able as independent constraints. It should be noted that atfrequencies lower than − Hz , the pulsar term cannotbe treated as random noise and pulsar distances are nec-essary to account for it (Yonemaru et al. 2018; Hisano et al.2019). Currently, the distance is not available for most pul-sars and we do not consider this frequency range here.Thus, the constraints Eqs. (22) and (23) are applicable for f GW & − Hz .Let us consider possible SMBH binaries at these cites.A SMBH with mass of . × M ⊙ is known to reside inthe Galactic Center and the possibility of the existence ofanother SMBH has been discussed (e.g. Oka et al. (2016)).If there exists a SMBH orbiting around the known SMBH,it could be a source of GWs. Assuming the period of thebinary motion to be
100 years , the upper bound of Eq.(22)translates into a upper bound on the companion mass of × M ⊙ .Concerning the M87, the mass of a SMBH is estimatedto be . × M ⊙ at the center of M87. it has been indicatedto have secondary SMBH and would be binary and GWsfrom such a potential pc-scale SMBH binary has been dis-cussed (Betcheldor et al. 2010; Yonemaru et al. 2016), while MNRAS , 1–7 (2019)
H., Kumamoto. et al. [t] -18.5-18-17.5-17-16.5-16-15.5 N u m b er log (p ·/ p) |b| < 10 ° |b| > 10 ° Figure 8. (Top) Position of 149 MSPs in the galactic coordi-nates with the spin-down rates represented by the color. (Bot-tom) Histogram of spin-down rates for MSPs in the Galactic plane( | b | <
10 deg , red) and outside ( | b | >
10 deg , blue). constraints on the amplitude of GWs emitted by a milli-pcscale SMBH binary in the PTA frequency bands has beenalready studied (Schutz & Ma 2016). Assuming the orbitalperiod to be
100 years , the upper bound of Eq.(23) resultsin a upper bound on the companion mass of × M ⊙ .In the above considerations, binaries are assumedto have circular orbits and inclination is zero degree(face-on). In our previous works (Yonemaru et al. 2018;Hisano et al. 2019), we estimated future constraints onultra-low-frequency GWs with 3,000 MSPs, which are ex-pected to be found by the SKA2. There, we found that GWswith Û h as small as ∼ × − sec − can be detected and thesecond SMBH mass as small as × M ⊙ could be probedin the case of circular orbits and zero inclination. In fact,the GW amplitude is sensitive to the eccentricity and thephase of the binary motion.Our simulation is based on the assumption that the dis-tribution of spin-down rate in the sky is isotropic in the ab-sense of GWs. However, considering the evolution of MSPs,it may not be the case. Neutron stars are mostly producedin Galactic plane and often have large peculiar velocitiesdue to the kick at supernovae (Hansen & Phinney 1997). Since MSPs with small values of Û p / p have long characteris-tic age ( p / Û p ), they may tend to be located far from the birthplace outside the Galactic plane. Thus, there is a possibil-ity that MSPs with large (small) Û p / p are populated outside(inside) of the Galactic plane, which induces the anisotropyof Û p / p distribution in the sky and results in systematics inour method.In the top panel of Fig. 8, the position of 149 MSPsin the galactic coordinates is shown with the indication ofspin-down rates. In the bottom panel, we show the his-togram of spin-down rates of MSPs within Galactic plane( | b | <
10 deg ) and outside ( | b | >
10 deg ), separately. Themean and standard deviation of the histograms are -17.4 and0.5 for | b | <
10 deg , and -17.5 and 0.4 for | b | >
10 deg , respec-tively. Thus, significant difference is not found between thetwo histograms and the systematics is considered to be neg-ligible. Here it should be noted that MSPs with | b | >
10 deg may reside within the Galactic disk but it cannot be knownwithout information on the distance.Finally, we note that the observed spin-down rates couldbe biased by other factors than GWs such as the Shklovskiieffect (Shklovskii 1970), the Galactic differential rotation(Damour & Taylor 1991; Rong & Tan 1999) and accelera-tion toward the Galactic disk (Nice & Taylor 1995). Al-though the biases from the Galactic differential rotation andacceleration toward the disk would have spatial correlationsin the sky, these effects are less significant ( ∆ ( Û p / p ) − for MSPs at
10 kpc) (Nice & Taylor 1995) and will be re-moved if the distance to MSPs is measured precisely in thefuture.
In this paper, we have placed constraints on GWs froma single source with ultra-low frequencies ( − Hz . f GW . − Hz ) by applying a method proposed inYonemaru et al. (2016) to observed milli-second pulsars(MSPs). This method is based on the statistics of spin-downrate distribution, where the skewness difference between twoMSP groups divided according to the position in the sky, isused to find a bias induced by GWs. We selected 149 MSPsfrom the ATNF pulsar catalog and calculated the skewnessdifference. By comparing with mock MSP data, we haveshown that the current MSP data is consistent with no GWsfrom any direction of the sky. Furthermore, we have derivedupper bounds on the time derivative of the GW amplitude Û h < . × − sec − and Û h < . × − sec − for the Galac-tic Center and M87, respectively. Consistent bounds werederived from the number of MSPs with negative spin-downrates. Approximating the GW amplitude as Û h ∼ π f GW h ,the bounds respectively translate into h < × − and h < × − for f GW = /(
100 yr ) . The constraints will im-prove by more than one order of magnitude with 3,000 MSPsin the SKA era (Yonemaru et al. 2018; Hisano et al. 2019). ACKNOWLEDGEMENTS
We thank George Hobbs for useful discussion. NY was fi-nancially supported by the Grant-in-Aid from the OverseasChallenge Program for Young Researchers of JSPS. SK is
MNRAS , 1–7 (2019) onstrains on ultra-low-frequency GWs REFERENCES
Abbott, B. P., et al. (the LIGO Scientific Collaboration, TheVirgo Collaboration), 2016, PhysRevLett. 116. 131103Acernese, F., et al. (VIRGO Collaboration), 2015, Classical Quan-tum Gravity 32, 024001Ade, P., et al. (The Simons Observatory Collaboration), 2018,arXiv: 1808.07445Aggarwal, K., et al. (the NANOGrav collaboration),2018, ac-cepted to Astrophysical Journal, arXiv:1812.11585Anholm, M., Ballmer, S., Creighton, J. D. E., Price, L. R. &Xavier, S., 2009, Phys. Rev. D, 79, 084030Babak, S., et al. (the EPTA collaboration), 2015, MNRAS, 455,2Batcheldor, D., Robinson, A., Axon, D. J., Perlman, E. S. & Mer-ritt, D., 2010, ApJ, 717, L6Bertotti, B., Carr B. J. & Rees M. J., 1983, MNRAS, 243, 945Blandford, R., Narayan, R., Romani & R. W., 1990, Journal ofAstrophysics and Astronomy, vol. 5, p. 369-388Damour, T. & Taylor, J. H., 1991, ApJ, 366, 501Detweiler, S., 2010, ApJ, 234, 1100Amaro-Seoane, P., et al., 2017, arXiv:1702.00786Hansen, B. M. S., & Phinney, E. S., 1997,Mon.Not.Roy.Astron.Soc. 291Hisano, S., Yonemaru, N., Kumamoto, H., & Takahashi, K., ac-cepted to MNRAS, arXiv:1902.04787KAGRA Collaboration, 2019, Nat. As., 3, 35Kawamura, S., et al., 2011, Classical and Quantum Gravity,28.9.094011Ishino, H., et al. 2016, Proceedings of the SPIE, 9904Jenet, F., et al. 2009, arXiv:0909.1058Kopeikin, S. M., 1997, Phys. Rev. D, 56, 4455Kramer, M., & Champion, D. J., 2013, Classical and QuantumGravity, 30, 22Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M., As-tronomical Journal, 129, 1993-2006, ”The ATNF Pulsar Cat-alogue”Manchester, R. N., et al. 2012, PASA, 30, 17Moore, C. D., Taylor, S., R. & Gair, J.,R., 2015, Classical andQuantum Gravity, 30, 32, 5Nice, D. J. & Taylor, J. H.,1995, ApJ, 441, 429Oka, T., Mizuno, R., Miura, K. & Tatekawa, S., 2016, ApJ, 816,L7Potapov, V. A., Ilyasov, Yu. P., Oreshko, V. V. & Rodin, A. E.,2003, Astron. Lett., 29, 241Rong, J., Xiao, N. & Tan, L., 1999, Science in China Series, A-Math, 42, 444Schutz, K. & Ma, C., 2016, MNRAS, 459, 1737Shklovskii, I. S., 1970, Soviet Astronomy, 13, 562Verbiest, J. P. W., et al. 2016, MNRAS, 458, 2, 1267Yonemaru, N., Kumamoto, H., Kuroyanagi, S., Takahashi, K. &Silk, J., 2016, PASJ, 68, 106 Yonemaru, N., Kumamoto, H., Kuroyanagi, S. & Takahashi, K.,2018, MNRAS, 478.1670Yu, Q. & Tremaine, S., 2003, Astrophys.J.599,1129-1138This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000