Search for ultralight scalar dark matter with NANOGrav pulsar timing arrays
KKOBE-COSMO-19-07
Search for ultralight scalar dark matter with NANOGrav pulsar timing arrays
Ryo Kato ∗ and Jiro Soda
1, † Department of Physics, Kobe University, Rokkodai 1-1, Kobe 657-8501, Japan (Dated: March 18, 2020)An ultralight scalar field is a candidate for the dark matter. The ultralight scalar dark matter with mass around10 − eV induces oscillations of the pulse arrival time in the sensitive frequency range of the pulsar timingarrays. We search for the ultralight scalar dark matter using the North American Nanohertz Observatory forGravitational Waves 11-year Data Set. We give the 95% confidence upper limit for the signal induced by theultralight scalar dark matter. In comparison with the published Bayesian upper limits on the amplitude of theultralight scalar dark matter obtained by Bayesian analysis using the Parkes Pulsar Timing Array 12-year data set(Porayko et al. 2018), we find three times stronger upper limit in the frequency range from 10 − . to 10 − . Hzwhich corresponds to the mass range from 9 . × − to 1 . × − eV. In terms of the energy density of thedark matter, we find that the energy density near the Earth is less than 7GeV / cm in the range from 10 − . to10 − . Hz (from 5 . × − to 2 . × − eV). The strongest upper limit on the the energy density is givenby 2GeV / cm at a frequency 10 − . Hz (corresponding to a mass 1 . × − eV). We find that the signal ofthe ultralight scalar dark matter can be explained by the solar system ephemeris effect. Also, we reveal that themodel of the solar system ephemeris effect prefers parameters which are contrary to the expectation that noisewill be reduced on all pulsars. PACS numbers:
I. INTRODUCTION
The dark matter problem is clearly one of the most important issues in modern cosmology. Recently, motivated by stringtheory, an ultralight scalar dark matter has been intensively studied [1, 2]. In particular, an ultralight scalar field with mass10 − eV can behave like the cold dark matter (CDM) on cosmological scales and resolve a cusp problem [3, 4]. In this article,we call it simply the fuzzy dark matter (FDM). The FDM can be treated as a classical scalar field because the occupationnumber of the FDM accounting for the energy density of the dark matter is very large. The main difference between FDM andCDM is that the pressure of the FDM is coherently oscillating, while that of CDM almost vanishes. Khmelnitsky and Rubakovhave pointed out that the effect of oscillating pressure might be detected with the pulsar timing arrays (PTAs) [5]. Indeed, theoscillation of the pressure induces the oscillation of the gravitational potential, and as a result, it induces the oscillation of thearrival time of the pulse passing through the gravitational potential.It would be worth noting that there exists experimental duality between gravitational wave and scalar dark matter detections.More precisely, the detection method for gravitational waves is useful for scalar field dark matter, and vice versa. The idea ofKhmelnitsky and Rubakov inspired us to use the gravitational wave interferometers for detecting scalar dark matter [6]. Recently,the importance of the reverse direction has been promoted and a novel constraint on GHz gravitational waves was obtained [7].Hence, it is importnt to investigate thoroughly the duality.An attempt to search for long wavelength gravitational waves with the PTAs composed of long-term observation of manypulsars was proposed in the articles [8–10]. Nowadays, the PTAs are most sensitive to the gravitational waves with a fewnanohertz frequency. There are three major pulsar timing projects aimed at observing the pulsars and searching for the gravi-tational waves: the European Pulsar Timing Array (EPTA) [11], the North American Nanohertz Observatory for GravitationalWaves (NANOGrav) [12], and the Parkes Pulsar Timing Array (PPTA) [13]. The collaboration of the three projects is called theInternational Pulsar Timing Array (IPTA) [14]. The point is that we can utilize the PTAs for searching the FDM.Porayko and Postnov [15] gave upper limits for the FDM with the Bayesian analysis using the NANOGrav 5-year Data Set.Moreover, Porayko et al. [16] gave upper limits for the FDM with the Bayesian and the Frequentist analyses using the PPTA 12-year Data Set. In this article, following the previous articles, we search for the FDM by the Bayesian analysis in the time domainusing the NANOGrav 11-year Data Set. We quantitatively investigate whether the ultralight scalar dark matter is detectable ornot using the Bayesian model selection approach. We clarify the prior dependence of constraints on the amplitude of the FDMand obtain three times stronger constraints on the amplitude of the FDM in the frequency range from 10 − . to 10 − . Hz. We ∗ Electronic address: [email protected] † Electronic address: [email protected] a r X i v : . [ a s t r o - ph . H E ] M a r also discuss how the results of Bayesian analysis depend on the the solar system ephemeris noise in the model describing theobservation data.This article is organized as follows. In Section II we describe a model of FDM signal. In Section III we briefly review theBayesian statistics, and explain how to use it for our analysis. In Section IV we describe the model, the data, and the functionused in the Bayesian analysis. In Section V we briefly review the MCMC simulation. In Section VI we describe the analysis ofthe white noise that is performed before the main analysis. In Section VII we summarize the results of Bayesian analysis usingNANOGrav 11-year data set. The last Section is devoted to conclusion. II. FDM SIGNAL
As we mentioned in the introduction, the oscillation of the scalar field with the mass m induces the oscillation of the gravita-tional potential and the oscillation of the arrival time of the pulse passing through the gravitational potential. The oscillation ofthe arrival time of the pulse induced by the FDM is given by [17] s ( t ) = − π f [ Ψ ( x e ) sin ( π f t + α ( x e )) − Ψ ( x p ) sin ( π f ( t − D ) + α ( x p ))] , (2.1)where f = m / π is the frequency, m is the mass, D is the distance between the pulsar and the earth, each x e and x p are the positionof the Earth and the pulsar, and α denotes the phase. Here, we used the gravitational potential Ψ ( x ) ≡ π G ρ ( x ) m , (2.2)where ρ ( x ) is the energy density of the dark matter. We will refer to this oscillation as the FDM signal.The parameters used in the Bayesian estimation are defined as follows: s ( t ) = − Ψ π f [ sin ( π f t + α e ) − sin ( π f t + α p )] , (2.3)where we assumed Ψ ≡ Ψ ( x e ) = Ψ ( x p ) and defined α e ≡ α ( x e ) , α p ≡ α ( x p ) − π f L . (2.4)Here, since we do not aim to estimate the distance, we put together the phase α ( x p ) and the distance L . In fact, the distance hasan uncertainty of tens to hundreds of parsec currently that is too large to determine the phase [18]. Since, as is mentioned in thearticle [16], the distance between the Earth and the pulsar L is not so large, it is reasonable to assume that the amplitudes at theearth Ψ ( x e ) and the pulsar Ψ ( x p ) are equal. The FDM signal is a simple sine wave superposition, and so some noise in the datamay have a similar waveform. Therefore, it is necessary to consider whether the FDM signal obtained by data analysis is theactual FDM signal or some noise.If all dark matter is occupied by the FDM, the predicted amplitude at the position of the Earth is Ψ ≃ . × − ( ρ . / cm )( − eV m ) , ≃ . × − ( ρ . / cm )( × − Hz f ) , (2.5)and f ≃ . × − Hz ( m − eV ) , (2.6)where ρ = . / cm is the estimated energy density of the dark matter at the position of the Earth [19]. III. BAYESIAN PARAMETER ESTIMATION AND MODEL COMPARISON
In this section we explain the Bayesian parameter estimation and the model comparison. For further details about the Bayesiandata analysis, see for example [20–22].
TABLE I: Interpretation of the Bayes factor B Evidence in favor of M against M >
150 Very strong
The purpose of the Bayesian parameter estimation is to estimate the posterior probability distribution p ( θ ∣ D ) of the parameters θ given the data D . Having the data, we can update our belief about the parameters using Bayes’ rule, namely p ( θ ∣ D ) = p ( D ∣ θ ) p ( θ ) p ( D ) , (3.1)where p ( θ ) is the prior probability distribution, p ( D ∣ θ ) is the likelihood function, and p ( D ) is the evidence. In the parameterestimation, the evidence can be regarded as a normalization constant which can be omitted, because it does not involve theparameter. We will define the data D , the parameter θ , the likelihood function p ( D ∣ θ ) , and the prior probability distribution p ( θ ) in Section IV.Using the posterior probability distribution for the amplitude of the FDM signal Ψ , we define the 95% upper limit R by ∫ R p ( Ψ ∣ D ) d Ψ = . . (3.2)We will calculate this 95% upper limit with the samples of the posterior probability distribution which is generated by theMCMC in Section VII A.For the Bayesian model comparison, it is often considered the Bayes factor: B = p ( D ∣M ) p ( D ∣M ) , (3.3)where M and M are competing models which assign a meaning to the parameters. The Table I gives the interpretation of theBayes factor in terms of the strength of the evidence [23–25].If M and M are nested models, we can use the Savege-Dickey density ratio to calculate the Bayes factor [24, 26–28]. Weassume that M is a model described in Section IV B, and M is a model in which the amplitude Ψ of the model M is fixed toa very small value 10 − . In this case, the model M can be regarded as a model with no FDM signal. Then the Savege-Dickeydensity ratio is B = p ( Ψ = − ∣M ) p ( Ψ = − ∣ D , M ) . (3.4)We will calculate this Bayes factor with the samples of the posterior probability distribution which is generated by the MCMCin Section IV. Specifically, we calculate the Bayes factor in multiple small bins around a fixed value Ψ , then derive the meanand unbiased standard deviation. IV. BAYESIAN ANALYSIS IN THE TIME DOMAIN
In this section, we explain the data D , the model M , and the parameter θ , and define the likelihood function p ( D ∣ θ ) , and theprior probability distribution p ( θ ) . A. Data
We used the NANOGrav 11-year data set [18] and chose six pulsars: PSRs J0613-0200, J1012+5307, J1600-3053,J1713+0747, J1744-1134, and J1909-3744. In this dataset, these pulsars have relatively good time-of-arrival (TOA) preci-sion and long observation period, which would be suitable for detecting the FDM signal which becomes larger as the frequencybecomes lower.The data D for the Bayesian analysis are timing residuals which are calculated by subtracting the timing model from the TOAs[29]. The fitting to make timing residual is called timing fit, which removes the currently well-known effects. The parametersincluded in the timing model are spin parameters, astrometry parameters, binary parameters (if pulser is binary), dispersionmeasure parameters, frequency dependency parameters, jump parameters, see the articles [18, 30] for details. Note that, byfitting the spin parameters, some of the low frequency signal that we are searching for in this article will be absorbed [31, 32].In the timing fit, TT (BIPM2015) is used for the Terrestrial Time, and JPL DE436 [33] is used for the planetary ephemeris. Itis known that the timing residuals change greatly depending on which planetary ephemeris is selected. To account for this error,the SSE noise described in Section IV B was first introduced to the model by the article [34].In order to obtain the timing residuals, we use the libstempo which is the PYTHON interface to TEMPO2 [35] timingpackage. We used the identical data set except for the parameter file of PSR J1713+0747. In the parameter file of thePSR J1713+0747, we changed only a parameter EPHEM from DE430 [36] to DE436, where this parameter specifies whichephemerides to be used. Then we used libstempo to fit the timing parameters of the PSR J1713+0747 and created a new param-eter file. We iterated the timing fitting five times, which would be sufficient for parameters to converge to certain values. Wefound that the LAMBDA and BETA parameters of the astrometry parameters shift by more than one sigma away from the oldervalue. This result can be inferred from the deviation of the position parameter due to the difference in ephemeris as pointed outin the article [37].
B. Model
Following the paper [30, 34], the model M of the timing residuals δ ttt for each pulsar includes the FDM signal, the linearizedtiming model noise, the red noise, the Solar System ephemeris noise, and the white noise. Each noise is described below.First, the linearized timing model noise characterize the inaccuracies of the timing fit. The linearized timing model noiseis represented by the product of a design matrix and small offsets. The design matrix describes the dependence of the timingresiduals on respective timing model parameters. We obtain the design matrix using the TEMPO2 via libstempo. The timingmodel parameters are listed in [18].Second, the red noise has most of their power at low frequencies in a given data set. The red noise is known to have achromatic(observing-frequency-independent) and chromatic (observing-frequency-dependent) components [38]. The achromatic compo-nents are thought to be caused by a random walk in one of the pulsar spin parameters [39–42] and contributions to TOAs by anasteroid belt around the pulsar [43]. The chromatic components are thought to be caused by the pulse propagating through theionized interstellar medium [38, 44]. Although the origins of red noise are various, simple power-law spectrum form is oftenused as the power spectral density. Specifically, the power spectral density of the red noise P ( f ) is defined as P ( f ) = A π ( ff yr ) − γ f − , (4.1)where f is a red noise frequency, f yr is 1yr − , A is a dimensionless amplitude of the red noise, and γ is a spectral index of thered noise.Third, the Solar System ephemeris noise characterizes the inaccuracy in the vector between the geocenter and the Solar Systembarycenter caused by the inaccuracies of the Solar System ephemeris. We will refer to this noise as the SSE noise. It is knownthat SSE noise affect upper limits and Bayes factors for amplitudes of the stochastic gravitational wave background [34]. Thecomponents of the SSE noise are considered to be the uncertainties of the planet mass and the planet orbit, and the uncertainty ofthe rotation rate around the ecliptic pole. The uncertainties of the planet mass are taken into account in Jupiter, Saturn, Uranusand Neptune. The uncertainties of the planet orbit are taken into account in Jupiter. We used the values and the data implementedin ENTERPRISE (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE) which is a pulsar timing analysiscode . Thus, regarding the SSE noise caused by the uncertainties of the planet orbit, the Jupiter mass is the value of the IAU2009 system of astronomical constants [45]. Also, the data describing the dependence of the Jupiter orbit on respective set-IIIparameters [46] is the same as ENTERPRISE. Regarding the SSE noise caused by the uncertainty of the rotation rate around theecliptic pole, the reference date corresponds to MJD 55197. http://vallis.github.io/libstempo https://bitbucket.org/psrsoft/tempo2.git We confirmed that each pulsar’s value of the chi-square and the degrees of freedom which can be derived by TEMPO2 are consistent with values listed in thefile ‘stats 11y 20180226.dat’, where this file is included in the data set and can be used as to see if TEMPO2 is properly constructed. Therefore, TEMPO2was installed as expected. https://github.com/nanograv/enterprise Finally, the white noise has the same power in all frequency region. The components of the white noise are EFAC, EQUAD,and ECORR. The EFAC characterize systematic errors of the TOA measurement uncertainties. The EQUAD is an additionalwhite noise. The EFAC and the EQUAD are used independently for each combination of receivers and backend systems. TheECORR differs from the EQUAD in that there is a correlation between the TOAs obtained during a single observation. ThisECORR characterizes pulse jitter caused by stochastic amplitude and phase variations in pulse, which correlates in a certainfrequency band and doesn’t correlate in time [44].
C. Likelihood Function
The likelihood for each pulsar can be written as [30, 47]: p ( δ ttt ∣ θθθ ) = √( π ) N TOA det ( CCC white ) det ( BBB ) det ( TTT T CCC − TTT + BBB − )× exp (− [ δ rrr T CCC − δ rrr − ( TTT T CCC − δ rrr ) T ( TTT T CCC − TTT + BBB − ) − TTT T CCC − δ rrr ]) , (4.2)where δ rrr ≡ δ ttt − sss − nnn SSE , TTT ≡ (
MMM FFF ) , BBB ≡ diag ( , , ⋯ ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ N TM , Ξ , Ξ , ⋯ Ξ N red ) , (4.3) δ ttt is the N TOA dimensional timing residuals, N TOA denotes the number of TOAs of the pulsar, sss is the N TOA dimensional FDMsignal, nnn
SSE is the N TOA dimensional linearized timing model noise,
MMM is a N TOA × N TM design matrix, N TM denotes the numberof the timing model parameters, FFF is a N TOA × N red matrix which has columns of alternating cosine and sine functions for thered noise, N red denotes the number of frequencies, 10 are the variance of the TM noise parameters, Ξ is the power spectrum ofthe corresponding frequency times the inverse of the observation period of the pulsar, and CCC white is the N TOA × N TOA covariancematrix which is composed by the EFAC, EQUAD, and ECORR. The likelihood function using multiple pulsars is representedby the product of the likelihood function of each pulsar.
D. Prior Probability Distribution
We use specific knowledge only for the mass errors of each planet as the prior information. Using the propagation of uncer-tainty law, the variances of Jupiter mass, Saturn mass, Neptune mass are calculated from the IAU 2009 system of astronomicalconstants. Also, the variance of Uranus mass is calculated from the values in the article [48] which is newer than the IAU 2009system of astronomical constants. Then we assume a normal distribution for the uncertainty of each planet mass and apply theobtained variances. For parameters without specific knowledge, we use a log-uniform distribution for parameters which areneed to be searched over several orders of magnitude with only positive values, and we use a uniform distribution for the otherparameters. The range of the log-uniform distribution and the uniform distribution is taken sufficiently wider than the value thatthe parameter would take. Regarding the amplitude of the FDM signal, we especially consider both cases of uniform distributionand log-uniform distribution as in the article [34]. The uniform distribution is used for the upper limit, and the log-uniformdistribution is used for the model comparison. The parameters and their prior probability distribution are given in Table II.As is done in [34, 49, 50], we analyze the white noises in advance before the main analysis. The resulting MCMC chains areused to calculate the value that maximizes the one-dimensional posterior probability distribution corresponding to each whitenoise, where this value is called the maximum a posteriori (MAP) value. The main analysis is performed by fixing the possiblevalues of the white noise parameter to the MAP value. See the section VI for the pre-analysis.
V. MARKOV CHAIN MONTE CARLO SIMULATION
The MCMC simulation can be used to generate samples from the posterior probability distribution. The MCMC method weused is called parallel tempering. In the parallel tempering method, a concept of temperature is introduced, and the MCMCsimulations of different temperatures are executed in parallel. The advantage of parallel tempering is that it is possible to
TABLE II: Prior Probability Distributionparameter description prior probability distributionFDM signal a Ψ amplitude Uniform[10 − , 10 − ] (for upper limit)logUniform[ − −
11] (for model comparison) f [ Hz ] frequency logUniform[log f − . f + . α e [ rad ] phase at Earth Uniform[0, 2 π ] α p [ rad ] phase at pulsar Uniform[0, 2 π ]red noise A amplitude logUniform[ − − γ spectral index Uniform[0 .
02, 6 . b δ M J [ M ⊙ ] mass error of Jupiter N (0, 1 . × − ) δ M S [ M ⊙ ] mass error of Saturn N (0, 8 . × − ) δ M U [ M ⊙ ] mass error of Uranus N (0, 5 . × − ) δ M N [ M ⊙ ] mass error of Neptune N (0, 7 . × − ) δ a J µ small offsets of set-III parameters Uniform[ − .
05, 0 . δ z [ rad / year ] rotation rate around ecliptic pole Uniform[ − − , 10 − ]white noiseEFAC EFAC parameter Uniform[0 . [ s ] EQUAD parameter logUniform[ − −
4] (for pre-analysis)ECORR [ s ] ECORR parameter logUniform[ − . −
4] (for pre-analysis) a f is the center frequency of the bin: f = { − , − . , ⋯ , − . } b µ is the subscript for the set-III parameters : µ = { , , ⋯ , } reduce the tendency of the samples of the posterior distribution to be trapped in a local minimum, compared to the Metropolis-Hastings method which is the one of the most famous MCMC methods [20]. We carry out the analysis using four temperatures T = . , . , . , [51] via PAL2 [52] whichis a Bayesian inference package for PTA and can include the PTMCMCSampler. Regarding models not implemented in thePAL2, the FDM signal is implemented like the continuous gravitational waves and the SSE noise is implemented like any othernoises. Following the article [49], we use adaptive Metropolis [53], single component adaptive Metropolis [54], and differentialevolution [55], as a proposal algorithm which is used to generate next samples using past samples. Furthermore, we also use asimple proposal algorithm to generate the next sample of each parameter by proposal distribution which is the same distributionas the probability distribution. All of these proposal algorithms are used in a single MCMC simulation and which one is usedis chosen randomly for each proposal in the MCMC simulation. In this article, we use the value written in the PAL2 for eachvariable used in PTMCMCSampler, unless specifically mentioned. VI. PRE-ANALYSIS
As is usual [34, 49, 50], in order to obtain the MAP values of the parameters of the white noise, we analyze the white noisefirst before the main analysis. By doing this, the number of the free parameters can be reduced in the main analysis. In thepre-analysis, we performed independent analysis for each pulsar, and we used the model which contains the red noise in additionto white noise. We ran the MCMC simulation with 10 iterations and removed the first 25% as a burn-in period, where theburn-in period is the period during which samples have not yet been obtained from the target distribution.Table III shows the results of the white noise and the red noise obtained by the pre-analysis. The obtained posterior distribution https://github.com/jellis18/PTMCMCSampler https://github.com/jellis18/PAL2 TABLE III: Result of pre-analysisparameter J0613-0200 J1012+5307 J1600-3053 J1713+0747 J1744-1134 J1909-3744EFAC noiseRcvr 800 GASP 1 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . Rcvr 800 GUPPI 1 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . Rcvr1 2 GASP 1 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . Rcvr1 2 GUPPI 1 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . L-wide ASP 1 . + . − . L-wide PUPPI 1 . + . − . S-wide ASP 1 . + . − . S-wide PUPPI 1 . + . − . EQUAD noise a Rcvr 800 GASP − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . Rcvr 800 GUPPI − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . Rcvr1 2 GASP − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . Rcvr1 2 GUPPI − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . L-wide ASP − . + . − . L-wide PUPPI − . + . − . S-wide ASP − . + . − . S-wide PUPPI − . + . − . ECORR noise a Rcvr 800 GASP − . + . − . − . + . − . − . + . + . − . + . − . − . + . − . − . + . − . Rcvr 800 GUPPI − . + . − . − . + . . − . + . − . − . + . − . − . + . − . − . + . − . Rcvr1 2 GASP − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . Rcvr1 2 GUPPI − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . L-wide ASP − . + . − . L-wide PUPPI − . + . − . S-wide ASP − . + . − . S-wide PUPPI − . + . − . red noise A a − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . − . + . − . γ . + . − . . + . − . . + . − . . + . − . . + . − . . − . − . a The values are expressed in logarithm with base 10. is expressed by the MAP value and the 95% confidence interval. At the Green Bank Observatory, the receivers are Rcvr 800and Rcvr1 2, and the backend system is GASP in early observations and GUPPI in later observations. Similarly, at the AreciboObservatory, the receivers are L-wide and S-wide, and the backend system is ASP in early observations and PUPPI in laterobservations. Please see the article [30] for details on the receiver and the backend system.
VII. RESULT
In this section we describe the upper limits on the amplitude of the FDM signal and how much the FDM signal is absorbedby other noises. All our result was calculated using six pulsars: PSRs J0613-0200, J1012+5307, J1600-3053, J1713+0747,J1744-1134, and J1909-3744 in the NANOGrav 11-year data set. ⇥ ⇥ ⇥ Frequency f [Hz] A m p li t ud e FDM+Red+SSE uniformFDM+Red+SSE log-uniformFDM+Red uniformFDM+Red log-uniformpredicted amplitude ⇥ ⇥ FDM mass m [eV] 10 ⇥ ⇥ ⇥ Frequency f [Hz]510152013 B a y e s F a c t o r B ± ± FDM+Red+SSE log-uniformFDM+Red log-uniform
PPTA2018 ⇥ ⇥ ⇥ A m p li t ud e FDM+Red+SSE uniformFDM+Red+SSE log-uniformFDM+Red uniformFDM+Red log-uniformpredicted amplitude ⇥ ⇥ FDM mass m [eV] 10 ⇥ ⇥ ⇥ Frequency f [Hz]510152013 B a y e s F a c t o r B ± ± FDM+Red+SSE log-uniformFDM+Red log-uniform
FIG. 1: Top: The 95% upper limits on the amplitude of the FDM signal Ψ using the NANOGrav 11-year data set. As a prior probabilitydistribution of the amplitude, the uniform prior was used for the black solid lines and the log-uniform prior was used for black dashed lines.The red lines are the upper limit obtained when the SSE noise is not included in the model describing the observed data, and the solid anddashed lines indicate that unifrom and loguniform were used, respectively. The bold black line is the upper limit obtained by the Bayesiananalysis of the PPTA data ( taken from Figure 3 in [16]). The green line is the predicted amplitude. The purple vertical line is the inverseof the observation period of the pulsars. Bottom: The values of the Bayes factor obtained when using log-uniform prior. The black and redindicate when the SSE noise is included in the model or not, respectively. To improve the visibility of the plot, when the value of the Bayesfactor exceeds 20, the upper triangle is used. A. Upper limits
We calculated the 95% confidence upper limits on the amplitude of the FDM signal by the Bayesian analysis. We ran all theMCMC simulation with 10 iterations and removed the first 25% as a burn-in period. As the prior probability distribution of theamplitude of the FDM signal, we considered two cases: the uniform distribution and the log-uniform distribution. The uniformprior was used to place the conservative upper limits, the log-uniform prior was used to calculate the Bayes factors, where theupper limits were calculated using the Eq. (3.2) and the Bayes factors were calculated using Eq. (3.4). In order to see the effectof including the SSE noise in the model on the results, we also calculated the upper limits and the Bayes factors when the SSEnoise is not included in the model.In Figure 1, we show the upper limits and the Bayes factors of the amplitude Ψ as a function of the frequency f and themass m . The relation between the frequency f and the mass m is given by Eq. (2.6). First, in the above plot, the black solidand dashed lines denote the upper limit using uniform prior and log-uniform prior, respectively. Here, we plotted the resultsobtained using log-uniform prior, but as we mentioned in the section IV D, we regard the results obtained by the uniform prior asconservative upper limits. The red solid and dashed lines denote the upper limit obtained when the SSE noise was not includedin the model, and using uniform-prior and log-uniform prior, respectively. The bold black line denotes the upper limit of theBayesian analysis obtained in [16] ( taken from Figure 3). The green line denotes the predicted amplitude of the FDM signalgiven by Eq. (2.5) with ρ = . / cm . The purple vertical lines denote the inverse of the observation periods of pulsars andcorresponds to PSRs J1744-1134, J1012+5307, J1909-3744, J1713+0747, J0613-0200, and J1600-3053 in order from the left.Next, in the bottom plot, the black and red dots denote the mean value of the Bayes factor using the model with and without theSSE noise, respectively. The unbiased standard deviation is used for error bars. Only to make the plot easier to see, when theBayes factor exceeds 20, it is represented by the upper triangle and the mean value and the unbiased standard deviations of theBayes factor is written above it.First we consider the red results obtained when the SSE noise is not included in the model. It turns out that the upper limitsfor the log-uniform prior gives stronger limits than for the uniform prior, but the difference is small. The reason the difference ⇥ ⇥ ⇥ Frequency f [Hz] E n e r g y d e n s i t y ⇢ [ G e V / c m ] FDM+Red+SSE uniformFDM+Red+SSE log-uniformFDM+Red uniformFDM+Red log-uniformpredicted amplitude ⇥ ⇥ FDM mass m [eV] 10 ⇥ ⇥ ⇥ Frequency f [Hz]510152013 B a y e s F a c t o r B ± ± FDM+Red+SSE log-uniformFDM+Red log-uniform
PPTA2018 ⇥ ⇥ ⇥ Frequency f [Hz] E n e r g y d e n s i t y ⇢ [ G e V / c m ] FDM+Red+SSE uniformFDM+Red+SSE log-uniformFDM+Red uniformFDM+Red log-uniformpredicted amplitude ⇥ ⇥ FDM mass m [eV] 10 ⇥ ⇥ ⇥ Frequency f [Hz]510152013 B a y e s F a c t o r B ± ± FDM+Red+SSE log-uniformFDM+Red log-uniform
FIG. 2: Top: The 95% upper limits on the energy density of the FDM ρ using the NANOGrav 11-year data set. This plot is the same as inFigure 1 except that the amplitude is converted to energy density. The bold black line denotes the upper limit obtained by the Bayesian analysisof the PPTA data ( taken from Figure 4 in [16]). Bottom: The values of the Bayes factor obtained when using log-uniform prior. This plot isthe same as in Figure 1. is small is that the Bayes factor exceeds 3 when the frequency becomes 10 − . Hz (1 . × − eV) or lower. According to theTable I, the Bayes factor exceeds 3 means that there is a signal that is somewhat similar to the FDM signal. Therefore, whicheverprior probability distribution is used, the value of the posterior probability distribution tends to be large at the parameter valuesof that signal, and as a result the upper limits does not change so much. From the above, it was found that the FDM signal ispositively supported in the wide frequency range. However, from a physical point of view, it is hard to think that FDM signalshas been found, because the upper limit obtained is about an order of magnitude greater than the expected amplitude and alsoFDM signals would have only a certain frequency.Next, we consider the black result obtained when the SSE noise is included in the model. This result is our main result.Compared with the case where SSE noise is included in the model, it can be seen that the upper limits of the amplitude ofthe FDM signal obtained using uniform distribution does not change much. However, the difference between the upper limitsobtained from the different prior probability distributions is large. The reason for this is that all Bayes factors are smaller than 3and in most cases they do not exceed 1. Therefore, in this case we conclude that no FDM signal has been detected. The red andblack results show that the Bayes factor is smaller when the SSE noise is considered. Thus, we conclude that the SSE noise canmimic the FDM signal. This result would be inferred from the result that the SSE noise can mimic the stochastic gravitationalwave background [34]. In comparison with the published Bayesian upper limits of the amplitude of the FDM signal using thePPTA 12-year data set [16], i.e. comparing the black and the bold black lines, we found that stronger upper limits were obtainedwhen the frequency was in the range from 10 − . to 10 − . Hz (from 9 . × − to 1 . × − eV). In this range, up to threetimes stronger upper limits were obtained, and in other region, about the same upper limits were obtained.It is also important to see the upper limit on the energy density of the dark matter near the Earth rather than the amplitudeof the FDM signal. Thus, we convert the amplitude of the FDM signal into the energy density using Eq. (2.5), and the resultis plotted in Figure 2. Note that the bold black line denotes the upper limit on the energy density with the Bayesian analysis in[16] (taken from Figure 4). As we can see from Figure 2, our main upper limit represented by the black line is 7 or less in therange from 10 − . to 10 − . Hz (from 5 . × − to 2 . × − eV) where we analyzed. The strongest upper limit on the theenergy density is 2GeV / cm at the frequency 10 − . Hz (1 . × − eV).0 TABLE IV: Red noise and SSE noise used for fixed noise analysisparameter J0613-0200 J1012+5307 J1600-3053 J1713+0747 J1744-1134 J1909-3744red noise A red a − . + . − . − . + . − . − . + . − . − . + . . − . + . − . − . + . − . γ . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . common to all pulsarsSSE noise δ M J − × − + . × − − . × − δ M S × − + . × − − . × − δ M U × − + . × − − . × − δ M N × − + . × − − . × − δ a J1 − . + . − . δ a J2 . + . − . δ a J3 − . + . − . δ a J4 − . + . − . δ a J5 . + . − . δ a J6 . + . − . δ z − . × − + . × − − . × − a The values are expressed in logarithm with base 10.
B. Fixed noise analysis
We analyzed the red noise and the SSE noise first and calculated the upper limits of the amplitude of the FDM signal usingthe obtained MAP values of the parameters. We ran the MCMC simulation with 10 iterations for analysis of red noise and SSEnoise and with 10 iterations for analysis of the FDM signal, and in both cases we removed the first 25% as a burn-in period. Theresult of the red noise and the SSE noise used for the fixed noise analysis are summarized in Table IV. The obtained posteriordistribution is expressed by the MAP value and the 95% confidence interval.As in the previous section, we calculated two cases of uniform and log-unifrom distributions as the prior probability distribu-tion of the amplitude of the FDM signal. The results are plotted in Figure 3 which is the similar plot as Figure 1. The solid anddashed lines indicate that the unifrom and the log-uniform prior were used, respectively.The reason for doing this analysis is to know how much the red noise and the SSE noise can absorb the signal of the FDMsignal. Note that this analysis is not intended to give the upper limits on the amplitude of the FDM signal. The FDM signal,the red noise and the noise induced by Jupiter in the SSE noise have similar waveforms and any of these can not be analyzedin advance. We consider that this analysis is not suitable for giving an upper limit to the amplitude of the FDM signal, and theobtained results are not regarded as the upper limits of the amplitude of the FDM signal.It can be seen from Figure 3 that the values of the upper limits are drastically smaller than the result obtained in the previoussection. In particular, when using the log-uniform prior, surprisingly, the upper limits are smaller than the predicted amplitudein some range. As for Bayesian factors, they are all smaller than 1, which is consistent with the fact that the upper limits arestrongly influenced by the prior probability distribution. From this result, it is inferred that the FDM signal is well absorbed bythe red noise and the SSE noise.In order to investigate the impact of analyzing the noise first, we made simulated noise using MAP values of the SSE noiseThen we calculated the Lomb-Scargle periodogram of the timing residual by subtracting it, where the Lomb-Scargle periodogramcan be used to search for periodic signals in non-uniformly spaced time series data [56, 57]. The reason for not subtracting thered noise from the original timing residual is that it is difficult to create the noise included in the actual data because the rednoise that we can create is only one realization of the stochastic process. For comparison, we also calculated the Lomb-Scargleperiodogram of the original timing residual and the simulated timing residual induced by the red noise only. In order to calculatethe Lomb-Scargle periodogram, we used Astropy [58, 59] which is a Python package for the astronomy. For the purpose of − × − × − × − Frequency f [Hz]10 − − − − − A m p li t ud e Ψ Red+SSE then FDM uniformRed+SSE then FDM log-uniformpredicted amplitude − × − × − FDM mass m [eV] − × − × − × − Frequency f [Hz]0 . . . . . B a y e s F a c t o r B Red+SSE then FDM log-uniform
FIG. 3: Similar plot to Figure 1. The FDM signal, the red noise and the noise induced by Jupiter in the SSE noise have similar waveforms andany of these can not be analyzed in advance. Therefore we do not regard this plot from such an analysis as the upper limits of the amplitude ofthe FDM signal. The purpose of this plot is to know how much the red noise and the SSE noise can absorb the signal of the FDM signal. expressing red noise, the original residual has short observation period and lacks frequency resolution, so that we made red noiseusing simulated observation data that the observation period is 10 days and the data points are every day.In Figure 4, we show the Lomb-Scargle periodogram, where the horizontal axis is the frequency and the vertical axis is theLomb-Scargle power. The black, red, and blue lines represent the periodogram of the original timing residual, the red noise only,and the timing residual subtracted by the SSE noise, respectively. As for PSR J1744-1134, a plot focused on other than the rednoise is also displayed below it, because the red noise is large and the other lines are difficult to see. The purple vertical linerepresents the inverse of the observation period.First, with regard to the black line representing the periodogram of the original timing residual, it can be seen that the power inthe frequency smaller than 10 − Hz is larger than the power in the frequency larger than 10 − Hz. Therefore, it can be understoodthat there is noise in the low frequency region, and the analysis of this article is considered to be meaningful.Second, with regard to the red line representing the periodogram of the red noise only, the red noise can be considered tocharacterize low frequency noise in PSRs J0613-0200 and J1012+5307, but it can not be considered so in the other pulsars. Inparticular, the amplitude of red noise is large in PSR J1744-1134, so that in the analysis in which the MAP value of the red noiseis fixed, it can be seen that this pulsar will contribute little to the analysis result of the amplitude of the FDM signal. The reasonwhy the red noise can be detected properly in PSRs J0613-0200 and J1012+5307 is that, according to the result of analyzingonly the red noise and the SSE noise, the posterior distribution of the red noise is obtained with a sharp peak. On the other hand,in the other pulsars, the posterior probability distribution of the red noise has a peak which spreads over a wide parameter space,which indicates that the MAP value has little meaning.Finally, with regard to the blue line representing the original timing residual after subtracting the SSE noise, it is found thatthe low frequency noise are reduced in PSRs J0613-0200 and J1909-3744. On the other hand, this is not the case with the otherpulsars, and it can be seen that the noise increases rather than decreases. We conclude that the reason why the upper limits aresmaller in fixed noise analysis than in the main analysis is mainly due to PSR J1909-3744, because only this pulsar has smallred noise and the low frequency noise is removed by the SSE noise. Conversely, the other pulsars have not improved much,which may suggest that fixed noise analysis is not successful. In particular, since the SSE noise can mimic the FDM signal, it isa problem that the model of the SSE noise may favor parameters which are contrary to the expectation that noise will be reducedon all pulsars. More precise measurements of the SSE would be needed to detect or exclude the FDM signal. However, theseissues are beyond the scope of this article.2 Frequency f [Hz]0.00.51.01.52.0 L o m b - S c a r g l e P o w e r [ s ]
1e 9 PSR J0613-0200ResidualRedResidual-SSE 10 Frequency f [Hz]0.00.20.40.60.81.0 L o m b - S c a r g l e P o w e r [ s ]
1e 8 PSR J1012+5307PSR J1012+5307 10 Frequency f [Hz]0123 L o m b - S c a r g l e P o w e r [ s ]
1e 10 PSR J1600-3053PSR J1600-305310 Frequency f [Hz]0246 L o m b - S c a r g l e P o w e r [ s ]
1e 10 PSR J1713+0747PSR J1713+0747 10 Frequency f [Hz]0246 L o m b - S c a r g l e P o w e r [ s ]
1e 8 PSR J1744-1134PSR J1744-1134 10 Frequency f [Hz]01234 L o m b - S c a r g l e P o w e r [ s ]
1e 10 PSR J1909-3744PSR J1909-374410 Frequency f [Hz]0.00.51.01.52.0 L o m b - S c a r g l e P o w e r [ s ]
1e 10 PSR J1744-1134PSR J1744-1134
FIG. 4: The Lomb-Scargle periodogram for PSRs J0613-0200, J1012+5307, J1600-3053, J1713+0747, J1744-1134, and J1909-3744. Theblack, red, and blue lines denote the Lomb-Scargle periodogram of the original timing residual, only the red noise, and the original timingresidual after subtracting the SSE noise. As for PSR J1744-1134, the red noise is large, hence a figure with a different scale is depicted belowthat figure to see the other lines.
VIII. CONCLUSION
We searched for the FDM by performing the Bayesian analysis in the time domain using the NANOGrav 11-year Data Set. InSection VII A, we gave the 95% confidence upper limit on the amplitude of the FDM signal. we found that probability that theFDM signal should be included in the model was less than 75% in all frequency region. Compared with the published Bayesianupper limit of the FDM signal using the PPTA 12-year data set [16], we found that our upper limit was up to 3 times stronger thanthe previous study when the frequency was in the range from 10 − . to 10 − . Hz (9 . × − to 1 . × − eV in terms of the3FDM mass). In other region, we also obtained the similar upper limit on the amplitude of the FDM signal. Since the amplitudeof FDM signal can be converted to the energy density of the dark matter near the Earth, it is easy to obtain the upper limit ofthe energy density. The upper limit on the energy density was lower than 7GeV / cm in the range from 10 − . to 10 − . Hz(from 5 . × − to 2 . × − eV) where we analyze. In particular, at a frequency of 10 − . Hz (a mass of 1 . × − eV),we obtained the strongest upper limit 2GeV / cm . In addition to the main analysis, we also investigated the case where the SSEnoise was not included in the model. In this case, we showed that we can not exclude the existence of the FDM signal, becausethe probability that the FDM signal should be included in the model was more than 75% in the frequency region 10 − . Hz orless. This results show that the Bayes factor is smaller when the SSE noise is considered. Thus, we conclude that the SSE noisecan mimic the FDM signal, which would be inferred from the result that the SSE noise can mimic the stochastic gravitationalwave background [34].In Section VII B, by analyzing the noise in advance, we examined how much the FDM signal was absorbed. In this case, weclarified that the probability that the FDM signal should be included in the model was much lower than 50% in all frequencyregion. Compared to our main analysis, we found that the upper limit on the amplitude of the FDM signal became very small.From this, it is expected that the FDM signal will be absorbed very well by analyzing the noise in advance. Then, we made asimulated noise from the parameters obtained by the analysis earlier, and investigated whether removing the SSE noise from theobserved data would reduce the power of the low frequency region of the data. We found that the power of the PSR J1909-3744has become smaller, and we conclude that this pulsar contributed to the improvement of the sensitivity. With other pulsars, wealso found that the power increased on the contrary. If the SSE noise is properly found, it is thought that the noise will be reducedfor all pulsars. Hence, this result seems strange. In order to solve this problem, it is necessary to reduce the uncertainty of theSSE by more accurate observation of the SSE.
Acknowledgments
R.K. was supported by Grant-in-Aid for JSPS Research Fellow and JSPS KAKENHI Grant Numbers 17J00496. J. S. was inpart supported by JSPS KAKENHI Grant Numbers JP17H02894, JP17K18778, JP15H05895, JP17H06359, JP18H04589. J. Sis also supported by JSPS Bilateral Joint Research Projects (JSPS-NRF collaboration) String Axion Cosmology. [1] Peter Svrcek and Edward Witten. Axions In String Theory.
JHEP , 06:051, 2006.[2] Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell. String Axiverse.
Phys. Rev. ,D81:123530, 2010.[3] Wayne Hu, Rennan Barkana, and Andrei Gruzinov. Fuzzy cold dark matter: The wave properties of ultralight particles.
Phys. Rev. Lett. ,85:1158–1161, Aug 2000.[4] David J. E. Marsh. Axion Cosmology.
Phys. Rept. , 643:1–79, 2016.[5] Andrei Khmelnitsky and Valery Rubakov. Pulsar timing signal from ultralight scalar dark matter.
Journal of Cosmology and AstroparticlePhysics , 2014(02):019–019, feb 2014.[6] Arata Aoki and Jiro Soda. Detecting ultralight axion dark matter wind with laser interferometers.
Int. J. Mod. Phys. , D26(07):1750063,2016.[7] Asuka Ito, Tomonori Ikeda, Kentaro Miuchi, and Jiro Soda. Probing GHz Gravitational Waves with Graviton-magnon Resonance. , 2019.[8] S. Detweiler. Pulsar timing measurements and the search for gravitational waves.
Astrophys. J. , 234:1100–1104, December 1979.[9] Roger W. Romani.
Timing a Millisecond Pulsar Array , pages 113–117. Springer Netherlands, 1989.[10] R. S. Foster and D. C. Backer. Constructing a pulsar timing array.
Astrophys. J. , 361:300–308, September 1990.[11] Michael Kramer and David J Champion. The european pulsar timing array and the large european array for pulsars.
Classical andQuantum Gravity , 30(22):224009, nov 2013.[12] M A McLaughlin. The north american nanohertz observatory for gravitational waves.
Classical and Quantum Gravity , 30(22):224008,nov 2013.[13] R. N. Manchester, G. Hobbs, M. Bailes, W. A. Coles, W. van Straten, M. J. Keith, R. M. Shannon, N. D. R. Bhat, A. Brown, S. G.Burke-Spolaor, and et al. The parkes pulsar timing array project.
Publications of the Astronomical Society of Australia , 30:e017, 2013.[14] G Hobbs, A Archibald, Z Arzoumanian, D Backer, M Bailes, N D R Bhat, M Burgay, S Burke-Spolaor, D Champion, I Cognard, W Coles,J Cordes, P Demorest, G Desvignes, R D Ferdman, L Finn, P Freire, M Gonzalez, J Hessels, A Hotan, G Janssen, F Jenet, A Jessner,C Jordan, V Kaspi, M Kramer, V Kondratiev, J Lazio, K Lazaridis, K J Lee, Y Levin, A Lommen, D Lorimer, R Lynch, A Lyne,R Manchester, M McLaughlin, D Nice, S Oslowski, M Pilia, A Possenti, M Purver, S Ransom, J Reynolds, S Sanidas, J Sarkissian,A Sesana, R Shannon, X Siemens, I Stairs, B Stappers, D Stinebring, G Theureau, R van Haasteren, W van Straten, J P W Verbiest,D R B Yardley, and X P You. The international pulsar timing array project: using pulsars as a gravitational wave detector.
Classical andQuantum Gravity , 27(8):084013, apr 2010.[15] N. K. Porayko and K. A. Postnov. Constraints on ultralight scalar dark matter from pulsar timing.
Phys. Rev. D , 90:062008, Sep 2014. [16] Nataliya K. Porayko et al. Parkes Pulsar Timing Array constraints on ultralight scalar-field dark matter. Phys. Rev. , D98(10):102002,2018.[17] Andrei Khmelnitsky and Valery Rubakov. Pulsar timing signal from ultralight scalar dark matter.
JCAP , 1402:019, 2014.[18] Zaven Arzoumanian et al. The NANOGrav 11-year Data Set: High-precision timing of 45 Millisecond Pulsars.
Astrophys. J. Suppl. ,235(2):37, 2018.[19] Fabrizio Nesti and Paolo Salucci. The dark matter halo of the milky way, AD 2013.
Journal of Cosmology and Astroparticle Physics ,2013(07):016–016, jul 2013.[20] Phil Gregory.
Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with Mathematica Support . Cam-bridge University Press, 2005.[21] J. Kruschke.
Doing Bayesian Data Analysis: A Tutorial Introduction with R . Elsevier Science, 2010.[22] A. Gelman, J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari, and D.B. Rubin.
Bayesian Data Analysis, Third Edition . Chapman &Hall/CRC Texts in Statistical Science. Taylor & Francis, 2013.[23] H. Jeffreys.
The Theory of Probability . Clarenson Press, 2nd edition, 1948.[24] Robert E. Kass and Adrian E. Raftery. Bayes factors.
Journal of the American Statistical Association , 90(430):773–795, 1995.[25] Adrian E. Raftery. Bayesian model selection in social research.
Sociological Methodology , 25:111–163, 1995.[26] A. O’Hagan and M. Kendall.
Kendall’s advanced theory of statistics: Vol. IIB . Edward Arnold, 1994.[27] James M. Dickey. The weighted likelihood ratio, linear hypotheses on normal location parameters.
The Annals of Mathematical Statistics ,42(1):204–223, 1971.[28] Isabella Verdinelli and Larry Wasserman. Computing bayes factors using a generalization of the savage-dickey density ratio.
Journal ofthe American Statistical Association , 90(430):614–618, 1995.[29] Russell T. Edwards, G. B. Hobbs, and R. N. Manchester. Tempo2, a new pulsar timing package. 2. The timing model and precisionestimates.
Mon. Not. Roy. Astron. Soc. , 372:1549–1574, 2006.[30] The NANOGrav Collaboration, Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, S. Chamberlin, S. Chatterjee, B. Christy, J. M. Cordes,N. Cornish, K. Crowter, P. B. Demorest, T. Dolch, J. A. Ellis, R. D. Ferdman, E. Fonseca, N. Garver-Daniels, M. E. Gonzalez, F. A.Jenet, G. Jones, M. L. Jones, V. M. Kaspi, M. Koop, M. T. Lam, T. J. W. Lazio, L. Levin, A. N. Lommen, D. R. Lorimer, J. Luo, R. S.Lynch, D. Madison, M. A. McLaughlin, S. T. McWilliams, D. J. Nice, N. Palliyaguru, T. T. Pennucci, S. M. Ransom, X. Siemens, I. H.Stairs, D. R. Stinebring, K. Stovall, J. K. Swiggum, M. Vallisneri, R. van Haasteren, Y. Wang, and W. Zhu. The NANOGrav Nine-yearData Set: Observations, Arrival Time Measurements, and Analysis of 37 Millisecond Pulsars.
Astrophys. J. , 813:65, November 2015.[31] R. Blandford, R. Narayan, and R. W. Romani. Arrival-time analysis for a millisecond pulsar.
Journal of Astrophysics and Astronomy ,5:369–388, December 1984.[32] Curt Cutler, Sarah Burke-Spolaor, Michele Vallisneri, Joseph Lazio, and Walid Majid. The Gravitational-Wave Discovery Space of PulsarTiming Arrays.
Phys. Rev. , D89(4):042003, 2014.[33] W. M. Folkner and R. S. Park. JPL planetary and Lunar ephemeris DE436.
Jet Propulsion Laboratory , September 2016.[34] Z. Arzoumanian et al. The NANOGrav 11-year Data Set: Pulsar-timing Constraints On The Stochastic Gravitational-wave Background.
Astrophys. J. , 859(1):47, 2018.[35] G. B. Hobbs, R. T. Edwards, and R. N. Manchester. tempo2, a new pulsar-timing package i. an overview.
Monthly Notices of the RoyalAstronomical Society , 369(2):655–672, 2006.[36] W. M. Folkner, J. G. Williams, D. H. Boggs, R. S. Park, and P. Kuchynka. The Planetary and Lunar Ephemerides DE430 and DE431.
Interplanetary Network Progress Report , 196:1–81, February 2014.[37] J. B. Wang et al. Comparison of Pulsar Positions from Timing and Very Long Baseline Astrometry.
Mon. Not. Roy. Astron. Soc. ,469(1):425–434, 2017.[38] M. T. Lam, J. M. Cordes, S. Chatterjee, Z. Arzoumanian, K. Crowter, P. B. Demorest, T. Dolch, J. A. Ellis, R. D. Ferdman, E. Fonseca,M. E. Gonzalez, G. Jones, M. L. Jones, L. Levin, D. R. Madison, M. A. McLaughlin, D. J. Nice, T. T. Pennucci, S. M. Ransom, R. M.Shannon, X. Siemens, I. H. Stairs, K. Stovall, J. K. Swiggum, and W. W. Zhu. The NANOGrav Nine-year Data Set: Excess Noise inMillisecond Pulsar Arrival Times.
Astrophys. J. , 834:35, January 2017.[39] P. E. Boynton, E. J. Groth, D. P. Hutchinson, G. P. Nanos, Jr., R. B. Partridge, and D. T. Wilkinson. Optical Timing of the Crab Pulsar,NP 0532.
Astrophys. J. , 175:217, July 1972.[40] E. J. Groth. Timing of the Crab Pulsar III. The Slowing Down and the Nature of the Random Process.
The Astrophysical JournalSupplement Series , 29, November 1975.[41] Ryan M. Shannon and James M. Cordes. Assessing the role of spin noise in the precision timing of millisecond pulsars.
The AstrophysicalJournal , 725(2):1607, 2010.[42] Andrew Lyne, George Hobbs, Michael Kramer, Ingrid Stairs, and Ben Stappers. Switched magnetospheric regulation of pulsar spin-down.
Science , 329(5990):408–412, 2010.[43] R. M. Shannon, J. M. Cordes, T. S. Metcalfe, T. J. W. Lazio, I. Cognard, G. Desvignes, G. H. Janssen, A. Jessner, M. Kramer, K. Lazaridis,M. B. Purver, B. W. Stappers, and G. Theureau. An Asteroid Belt Interpretation for the Timing Variations of the Millisecond PulsarB1937+21.
Astrophys. J. , 766:5, March 2013.[44] M. T. Lam, J. M. Cordes, S. Chatterjee, Z. Arzoumanian, K. Crowter, P. B. Demorest, T. Dolch, J. A. Ellis, R. D. Ferdman, E. F.Fonseca, M. E. Gonzalez, G. Jones, M. L. Jones, L. Levin, D. R. Madison, M. A. McLaughlin, D. J. Nice, T. T. Pennucci, S. M. Ransom,X. Siemens, I. H. Stairs, K. Stovall, J. K. Swiggum, and W. W. Zhu. The NANOGrav Nine-year Data Set: Noise Budget for PulsarArrival Times on Intraday Timescales.
Astrophys. J. , 819:155, March 2016.[45] Brian Luzum, Nicole Capitaine, Agn`es Fienga, William Folkner, Toshio Fukushima, James Hilton, Catherine Hohenkerk, George Krasin-sky, G´erard Petit, Elena Pitjeva, Michael Soffel, and Patrick Wallace. The iau 2009 system of astronomical constants: the report of theiau working group on numerical standards for fundamental astronomy.
Celestial Mechanics and Dynamical Astronomy , 110(4):293, Jul2011. [46] D. Brouwer and G.M. Clemence. Methods of celestial mechanics . Academic Press, 1961.[47] S. R. Taylor, L. Lentati, S. Babak, P. Brem, J. R. Gair, A. Sesana, and A. Vecchio. All correlations must die: Assessing the significanceof a stochastic gravitational-wave background in pulsar-timing arrays.
Phys. Rev. , D95(4):042002, 2017.[48] R. A. Jacobson. The orbits of the uranian satellites and rings, the gravity field of the uranian system, and the orientation of the pole ofuranus.
The Astronomical Journal , 148(5):76, 2014.[49] Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, S. J. Chamberlin, S. Chatterjee, J. M. Cordes, P. B. Demorest, X. Deng, T. Dolch, J. A.Ellis, R. D. Ferdman, N. Garver-Daniels, F. Jenet, G. Jones, V. M. Kaspi, M. Koop, M. T. Lam, T. J. W. Lazio, A. N. Lommen, D. R.Lorimer, J. Luo, R. S. Lynch, D. R. Madison, M. A. McLaughlin, S. T. McWilliams, D. J. Nice, N. Palliyaguru, T. T. Pennucci, S. M.Ransom, A. Sesana, X. Siemens, I. H. Stairs, D. R. Stinebring, K. Stovall, J. Swiggum, M. Vallisneri, R. van Haasteren, Y. Wang, , W. W.Zhu, and NANOGrav Collaboration. Gravitational waves from individual supermassive black hole binaries in circular orbits: Limits fromthe north american nanohertz observatory for gravitational waves.
The Astrophysical Journal , 794(2):141, 2014.[50] Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, S. J. Chamberlin, S. Chatterjee, B. Christy, J. M. Cordes, N. J. Cornish, K. Crowter, P. B.Demorest, X. Deng, T. Dolch, J. A. Ellis, R. D. Ferdman, E. Fonseca, N. Garver-Daniels, M. E. Gonzalez, F. Jenet, G. Jones, M. L.Jones, V. M. Kaspi, M. Koop, M. T. Lam, T. J. W. Lazio, L. Levin, A. N. Lommen, D. R. Lorimer, J. Luo, R. S. Lynch, D. R. Madison,M. A. McLaughlin, S. T. McWilliams, C. M. F. Mingarelli, D. J. Nice, N. Palliyaguru, T. T. Pennucci, S. M. Ransom, L. Sampson, S. A.Sanidas, A. Sesana, X. Siemens, J. Simon, I. H. Stairs, D. R. Stinebring, K. Stovall, J. Swiggum, S. R. Taylor, M. Vallisneri, R. vanHaasteren, Y. Wang, W. W. Zhu, and The NANOGrav Collaboration. The nanograv nine-year data set: Limits on the isotropic stochasticgravitational wave background.
The Astrophysical Journal , 821(1):13, 2016.[51] Justin Ellis and Rutger van Haasteren. jellis18/ptmcmcsampler: Official release, October 2017.[52] Justin Ellis and Rutger van Haasteren. jellis18/pal2: Pal2, January 2017.[53] Heikki Haario, Eero Saksman, and Johanna Tamminen. An adaptive metropolis algorithm.
Bernoulli , 7(2):223–242, 2001.[54] Heikki Haario, Eero Saksman, and Johanna Tamminen. Componentwise adaptation for high dimensional mcmc.
Computational Statistics ,20(2):265–273, 2005.[55] Cajo JF Ter Braak. A markov chain monte carlo version of the genetic algorithm differential evolution: easy bayesian computing for realparameter spaces.
Statistics and Computing , 16(3):239–249, 2006.[56] N. R. Lomb. Least-squares frequency analysis of unequally spaced data.
Astrophysics & Space Science , 39:447–462, February 1976.[57] J. D. Scargle. Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data.
Astrophys.J. , 263:835–853, December 1982.[58] Astropy Collaboration, T. P. Robitaille, E. J. Tollerud, P. Greenfield, M. Droettboom, E. Bray, T. Aldcroft, M. Davis, A. Ginsburg, A. M.Price-Whelan, W. E. Kerzendorf, A. Conley, N. Crighton, K. Barbary, D. Muna, H. Ferguson, F. Grollier, M. M. Parikh, P. H. Nair, H. M.Unther, C. Deil, J. Woillez, S. Conseil, R. Kramer, J. E. H. Turner, L. Singer, R. Fox, B. A. Weaver, V. Zabalza, Z. I. Edwards, K. AzaleeBostroem, D. J. Burke, A. R. Casey, S. M. Crawford, N. Dencheva, J. Ely, T. Jenness, K. Labrie, P. L. Lim, F. Pierfederici, A. Pontzen,A. Ptak, B. Refsdal, M. Servillat, and O. Streicher. Astropy: A community Python package for astronomy.
A&A , 558:A33, October2013.[59] Astropy Collaboration, A. M. Price-Whelan, B. M. Sip˝ocz, H. M. G¨unther, P. L. Lim, S. M. Crawford, S. Conseil, D. L. Shupe, M. W.Craig, N. Dencheva, A. Ginsburg, J. T. VanderPlas, L. D. Bradley, D. P´erez-Su´arez, M. de Val-Borro, T. L. Aldcroft, K. L. Cruz, T. P.Robitaille, E. J. Tollerud, C. Ardelean, T. Babej, Y. P. Bach, M. Bachetti, A. V. Bakanov, S. P. Bamford, G. Barentsen, P. Barmby,A. Baumbach, K. L. Berry, F. Biscani, M. Boquien, K. A. Bostroem, L. G. Bouma, G. B. Brammer, E. M. Bray, H. Breytenbach,H. Buddelmeijer, D. J. Burke, G. Calderone, J. L. Cano Rodr´ıguez, M. Cara, J. V. M. Cardoso, S. Cheedella, Y. Copin, L. Corrales,D. Crichton, D. D’Avella, C. Deil, ´E. Depagne, J. P. Dietrich, A. Donath, M. Droettboom, N. Earl, T. Erben, S. Fabbro, L. A. Ferreira,T. Finethy, R. T. Fox, L. H. Garrison, S. L. J. Gibbons, D. A. Goldstein, R. Gommers, J. P. Greco, P. Greenfield, A. M. Groener,F. Grollier, A. Hagen, P. Hirst, D. Homeier, A. J. Horton, G. Hosseinzadeh, L. Hu, J. S. Hunkeler, ˇZ. Ivezi´c, A. Jain, T. Jenness,G. Kanarek, S. Kendrew, N. S. Kern, W. E. Kerzendorf, A. Khvalko, J. King, D. Kirkby, A. M. Kulkarni, A. Kumar, A. Lee, D. Lenz,S. P. Littlefair, Z. Ma, D. M. Macleod, M. Mastropietro, C. McCully, S. Montagnac, B. M. Morris, M. Mueller, S. J. Mumford, D. Muna,N. A. Murphy, S. Nelson, G. H. Nguyen, J. P. Ninan, M. N¨othe, S. Ogaz, S. Oh, J. K. Parejko, N. Parley, S. Pascual, R. Patil, A. A. Patil,A. L. Plunkett, J. X. Prochaska, T. Rastogi, V. Reddy Janga, J. Sabater, P. Sakurikar, M. Seifert, L. E. Sherbert, H. Sherwood-Taylor,A. Y. Shih, J. Sick, M. T. Silbiger, S. Singanamalla, L. P. Singer, P. H. Sladen, K. A. Sooley, S. Sornarajah, O. Streicher, P. Teuben, S. W.Thomas, G. R. Tremblay, J. E. H. Turner, V. Terr´on, M. H. van Kerkwijk, A. de la Vega, L. L. Watkins, B. A. Weaver, J. B. Whitmore,J. Woillez, V. Zabalza, and Astropy Contributors. The Astropy Project: Building an Open-science Project and Status of the v2.0 CorePackage.