Search for ultralight scalar dark matter with NANOGrav pulsar timing arrays

TL;DR

This study searches for ultralight scalar dark matter using pulsar timing data, setting new upper limits on its amplitude and energy density, and finds that observed signals may be explained by solar system effects.

Contribution

It provides the first 11-year NANOGrav data analysis setting upper limits on ultralight scalar dark matter, improving constraints over previous studies and analyzing potential solar system effects.

Findings

Set 95% confidence upper limits on scalar dark matter amplitude.

Improved constraints compared to previous Bayesian limits.

Identified solar system ephemeris effects as a possible signal source.

Abstract

An ultralight scalar field is a candidate for the dark matter. The ultralight scalar dark matter with mass around $10^{-23}\,{\rm eV}$ induces oscillations of the pulse arrival time in the sensitive frequency range of the pulsar timing arrays. We search for the ultralight scalar dark matter using the North American Nanohertz Observatory for Gravitational Waves 11-year Data Set. We give the 95% confidence upper limit for the signal induced by the ultralight scalar dark matter. In comparison with the published Bayesian upper limits on the amplitude of the ultralight 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^{-8.34}$ to $10^{-8.19}\,{ \rm Hz}$ which corresponds to the mass range from $9.45\times10^{-24}$ to $1.34\times10^{-23}\,{\rm eV}$. In terms of the energy density of the dark matter, we find that the energy density near the Earth is less than $7\,{\rm GeV/cm^3}$ in the range from $10^{-8.55}$ to $10^{-8.01}\,{ \rm Hz}$ (from $5.83\times10^{-24}$ to $2.02\times10^{-23}\,{\rm eV}$). The strongest upper limit on the the energy density is given by $2\,{\rm GeV/cm^3}$ at a frequency $10^{-8.28}\,{ \rm Hz}$ (corresponding to a mass $1.09\times10^{-23}\,{\rm eV}$). We find that the signal of the ultralight scalar dark matter can be explained by the solar system ephemeris effect. Also, we reveal that the model of the solar system ephemeris effect prefers parameters which are contrary to the expectation that noise will be reduced on all pulsars.

Figures

Tables

Table 1. Table 1: Interpretation of the Bayes factor

$B_{12}$	Evidence in favor of ${ℳ}_1$ against ${ℳ}_2$
1-3	Not worth more than a bare mention
3-20	Positive
20-150	Strong
$>$150	Very strong

Table 2. Table 2: Prior Probability Distribution

parameter	description	prior probability distribution
FDM signal666$f$ is the center frequency of the bin: $f=\{{10}^{-9},{10}^{-8.97},\cdots ,{10}^{-8.01}\}$
$\Psi$	amplitude	Uniform[${10}^{-18}$,  ${10}^{-11}$] (for upper limit)
		logUniform[$-18$,  $-11$] (for model comparison)
$f\left[\mathrm{Hz}\right]$	frequency	logUniform[$\log f-0.015$,  $\log f+0.015$]
${\alpha }_e\left[\mathrm{rad}\right]$	phase at Earth	Uniform[$0$,  $2\pi$]
${\alpha }_p\left[\mathrm{rad}\right]$	phase at pulsar	Uniform[$0$,  $2\pi$]
red noise
$A$	amplitude	logUniform[$-20$,  $-11$]
$\gamma$	spectral index	Uniform[$0.02$,  $6.98$]
SSE noise777$\mu$ is the subscript for the set-III parameters : $\mu =\{1,2,\cdots ,6\}$
$\delta M^J\left[M_{\odot }\right]$	mass error of Jupiter	$𝒩$($0$,  $1.55\times {10}^{-11}$)
$\delta M^S\left[M_{\odot }\right]$	mass error of Saturn	$𝒩$($0$,  $8.17\times {10}^{-12}$)
$\delta M^U\left[M_{\odot }\right]$	mass error of Uranus	$𝒩$($0$,  $5.72\times {10}^{-11}$)
$\delta M^N\left[M_{\odot }\right]$	mass error of Neptune	$𝒩$($0$,  $7.96\times {10}^{-11}$)
$\delta a_{\mu }^J$	small offsets of set-III parameters	Uniform[$-0.05$,  $0.05$]
$\delta z\left[\mathrm{rad}/\mathrm{year}\right]$	rotation rate around ecliptic pole	Uniform[$-{10}^{-9}$,  ${10}^{-9}$]
white noise
EFAC	EFAC parameter	Uniform[$0.001$,  $10$] (for pre-analysis)
EQUAD$\left[s\right]$	EQUAD parameter	logUniform[$-10$,  $-4$]  (for pre-analysis)
ECORR$\left[s\right]$	ECORR parameter	logUniform[$-8.5$,  $-4$]  (for pre-analysis)

Table 3. Table 3: Result of pre-analysis

parameter	J0613-0200	J1012+5307	J1600-3053	J1713+0747	J1744-1134	J1909-3744
EFAC noise
Rcvr_800_GASP	${1.096}_{-0.059}^{+0.059}$	${1.138}_{-0.042}^{+0.049}$	${1.19}_{-0.11}^{+0.11}$	${1.134}_{-0.054}^{+0.050}$	${1.168}_{-0.064}^{+0.073}$	${0.984}_{-0.060}^{+0.081}$
Rcvr_800_GUPPI	${1.164}_{-0.032}^{+0.041}$	${1.182}_{-0.023}^{+0.023}$	${1.117}_{-0.040}^{+0.025}$	${1.066}_{-0.020}^{+0.023}$	${1.070}_{-0.024}^{+0.027}$	${1.041}_{-0.019}^{+0.023}$
Rcvr1_2_GASP	${1.061}_{-0.053}^{+0.056}$	${1.054}_{-0.041}^{+0.044}$	${1.155}_{-0.148}^{+0.096}$	${1.087}_{-0.060}^{+0.056}$	${0.991}_{-0.065}^{+0.080}$	${0.977}_{-0.060}^{+0.054}$
Rcvr1_2_GUPPI	${1.080}_{-0.023}^{+0.023}$	${1.081}_{-0.018}^{+0.021}$	${1.062}_{-0.019}^{+0.018}$	${1.040}_{-0.014}^{+0.015}$	${1.079}_{-0.026}^{+0.028}$	${1.049}_{-0.013}^{+0.017}$
L-wide_ASP				${1.013}_{-0.048}^{+0.062}$
L-wide_PUPPI				${1.091}_{-0.021}^{+0.026}$
S-wide_ASP				${1.114}_{-0.058}^{+0.066}$
S-wide_PUPPI				${1.110}_{-0.034}^{+0.038}$
EQUAD noise101010The values are expressed in logarithm with base 10.
Rcvr_800_GASP	$-{8.5}_{-1.5}^{+1.7}$	$-{6.68}_{-3.22}^{+0.18}$	$-{7.7}_{-2.2}^{+1.5}$	$-{6.97}_{-2.93}^{+0.18}$	$-{6.469}_{-0.085}^{+0.079}$	$-{6.620}_{-0.056}^{+0.046}$
Rcvr_800_GUPPI	$-{6.677}_{-0.188}^{+0.099}$	$-{6.348}_{-0.059}^{+0.062}$	$-{6.45}_{-3.45}^{+0.14}$	$-{7.14}_{-2.75}^{+0.14}$	$-{6.598}_{-0.045}^{+0.036}$	$-{7.344}_{-0.119}^{+0.090}$
Rcvr1_2_GASP	$-{9.74}_{-0.17}^{+3.16}$	$-{6.51}_{-3.38}^{+0.19}$	$-{6.34}_{-3.54}^{+0.17}$	$-{7.29}_{-2.59}^{+0.16}$	$-{6.382}_{-0.097}^{+0.072}$	$-{7.40}_{-2.41}^{+0.18}$
Rcvr1_2_GUPPI	$-{9.65}_{-0.27}^{+2.75}$	$-{6.50}_{-0.25}^{+0.11}$	$-{8.3}_{-1.7}^{+1.2}$	$-{8.01}_{-1.93}^{+0.24}$	$-{6.681}_{-0.051}^{+0.041}$	$-{8.04}_{-1.88}^{+0.19}$
L-wide_ASP				$-{7.54}_{-0.23}^{+0.11}$
L-wide_PUPPI				$-{7.89}_{-1.78}^{+0.13}$
S-wide_ASP				$-{8.46}_{-1.48}^{+0.93}$
S-wide_PUPPI				$-{7.64}_{-2.24}^{+0.16}$
ECORR noisea
Rcvr_800_GASP	$-{7.50}_{-0.96}^{+0.69}$	$-{8.08}_{-0.38}^{+1.37}$	$-{8.472}_{+0.019}^{+2.269}$	$-{7.45}_{-1.02}^{+0.49}$	$-{6.63}_{-1.38}^{+0.25}$	$-{7.95}_{-0.52}^{+0.76}$
Rcvr_800_GUPPI	$-{6.74}_{-0.32}^{+0.26}$	$-{8.480}_{0.010}^{+1.517}$	$-{6.20}_{-0.24}^{+0.13}$	$-{6.609}_{-0.142}^{+0.092}$	$-{6.40}_{-0.17}^{+0.12}$	$-{7.19}_{-0.24}^{+0.16}$
Rcvr1_2_GASP	$-{7.94}_{-0.51}^{+1.43}$	$-{6.59}_{-1.83}^{+0.28}$	$-{6.74}_{-1.70}^{+0.31}$	$-{7.11}_{-0.54}^{+0.25}$	$-{6.16}_{-0.20}^{+0.14}$	$-{8.31}_{-0.17}^{+0.89}$
Rcvr1_2_GUPPI	$-{6.77}_{-1.64}^{+0.15}$	$-{6.55}_{-1.45}^{+0.15}$	$-{6.81}_{-0.34}^{+0.16}$	$-{7.115}_{-0.091}^{+0.067}$	$-{6.372}_{-0.087}^{+0.090}$	$-{7.113}_{-0.079}^{+0.066}$
L-wide_ASP				$-{6.97}_{-0.11}^{+0.14}$
L-wide_PUPPI				$-{7.061}_{-0.069}^{+0.087}$
S-wide_ASP				$-{6.95}_{-0.18}^{+0.14}$
S-wide_PUPPI				$-{7.03}_{-0.10}^{+0.11}$
red noise
$A^a$	$-{13.15}_{-0.68}^{+0.17}$	$-{12.67}_{-0.16}^{+0.14}$	$-{13.45}_{-6.36}^{+0.16}$	$-{14.46}_{-3.74}^{+0.39}$	$-{13.44}_{-3.75}^{+0.33}$	$-{14.13}_{-1.96}^{+0.35}$
$\gamma$	${1.26}_{-0.99}^{+2.08}$	${0.96}_{-0.61}^{+0.79}$	${0.17}_{-0.10}^{+6.49}$	${2.3}_{-1.6}^{+3.9}$	${2.9}_{-1.6}^{+3.9}$	${6.899}_{-5.295}^{-0.021}$

Table 4. Table 4: Red noise and SSE noise used for fixed noise analysis

parameter	J0613-0200	J1012+5307	J1600-3053	J1713+0747	J1744-1134	J1909-3744
red noise
$A_{\mathrm{red}}$111111The values are expressed in logarithm with base 10.	$-{13.20}_{-0.73}^{+0.20}$	$-{12.68}_{-0.16}^{+0.14}$	$-{13.51}_{-6.34}^{+0.10}$	$-{19.939}_{0.055}^{+5.436}$	$-{13.51}_{-2.31}^{+0.33}$	$-{14.05}_{-5.77}^{+0.17}$
$\gamma$	${1.36}_{-0.98}^{+2.33}$	${1.07}_{-0.63}^{+0.74}$	${0.0897}_{-0.0069}^{+6.6316}$	${1.8}_{-1.6}^{+4.9}$	${2.8}_{-1.2}^{+4.0}$	${0.74}_{-0.62}^{+5.97}$
	common to all pulsars
SSE noise
$\delta M^J$	$-2\times 10_{}^{-12}{}_{-3.0\times {10}^{-11}}{}^{+3.3\times {10}^{-11}}$
$\delta M^S$	$1\times 10_{}^{-12}{}_{-1.6\times {10}^{-11}}{}^{+1.6\times {10}^{-11}}$
$\delta M^U$	$1\times 10_{}^{-11}{}_{-1.2\times {10}^{-10}}{}^{+1.0\times {10}^{-10}}$
$\delta M^N$	$1\times 10_{}^{-11}{}_{-1.7\times {10}^{-10}}{}^{+1.5\times {10}^{-10}}$
$\delta a_1^J$	$-{0.0046}_{-0.0099}^{+0.0109}$
$\delta a_2^J$	${0.004}_{-0.018}^{+0.019}$
$\delta a_3^J$	$-{0.0128}_{-0.0083}^{+0.0071}$
$\delta a_4^J$	$-{0.0123}_{-0.0077}^{+0.0117}$
$\delta a_5^J$	${0.0044}_{-0.0087}^{+0.0080}$
$\delta a_6^J$	${0.008}_{-0.019}^{+0.012}$
$\delta z$	$-2.0\times 10_{}^{-10}{}_{-7.4\times {10}^{-10}}{}^{+1.15\times {10}^{-09}}$

Equations

$$s(t)=-\frac{1}{2\pi f}[\Psi ({\mathbf{x}}_{\mathbf{e}})\sin (2\pi ft+\alpha ({\mathbf{x}}_{\mathbf{e}}))-\Psi ({\mathbf{x}}_{\mathbf{p}})\sin (2\pi f(t-D)+\alpha ({\mathbf{x}}_{\mathbf{p}}))],$$

$$\Psi (\mathbf{x})\equiv \frac{\pi G\rho (\mathbf{x})}{m^2},$$

$$s(t)=-\frac{\Psi }{2\pi f}[\sin (2\pi ft+{\alpha }_e)-\sin (2\pi ft+{\alpha }_p)],$$

$${\alpha }_e\equiv \alpha ({\mathbf{x}}_{\mathbf{e}}),{\alpha }_p\equiv \alpha ({\mathbf{x}}_{\mathbf{p}})-2\pi fL.$$

$$\Psi$$

$$f\simeq 4.84\times 10^{-9}\mathrm{Hz}\left(\frac{m}{10^{-23}\mathrm{eV}}\right),$$

$$p(\theta |D)=\frac{p(D|\theta )p(\theta )}{p(D)},$$

$${\int }_0^{R}p(\Psi |D)d\Psi =\mathrm{0.95.}$$

$$B_{12}=\frac{p(D|{\mathcal{M}}_1)}{p(D|{\mathcal{M}}_2)},$$

$$B_{12}=\frac{p(\Psi =10^{-18}|{\mathcal{M}}_1)}{p(\Psi =10^{-18}|D,{\mathcal{M}}_1)}.$$

$$P(f)=\frac{A^2}{12{\pi }^2}{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{3-\gamma }{f}^{-3},$$

$$p(\delta \mathbf{t}|\mathbf{\theta })$$

$$\delta \mathbf{r}$$

$$\mathbf{T}$$

$$\mathbf{B}$$

Full text

††preprint: KOBE-COSMO-19-07

I Introduction

The dark matter problem is clearly one of the most important issues in modern cosmology. Recently, motivated by string theory, an ultralight scalar dark matter has been intensively studied [1, 2]. In particular, an ultralight scalar field with mass $10^{-23}\,{\rm 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 occupation number of the FDM accounting for the energy density of the dark matter is very large. The main difference between FDM and CDM is that the pressure of the FDM is coherently oscillating, while that of CDM almost vanishes. Khmelnitsky and Rubakov have pointed out that the effect of oscillating pressure might be detected with the pulsar timing arrays (PTAs) [5]. Indeed, the oscillation of the pressure induces the oscillation of the gravitational potential, and as a result, it induces the oscillation of the arrival 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 of Khmelnitsky 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 many pulsars was proposed in the articles [8, 9, 10]. Nowadays, the PTAs are most sensitive to the gravitational waves with a few nanohertz frequency. There are three major pulsar timing projects aimed at observing the pulsars and searching for the gravitational waves: the European Pulsar Timing Array (EPTA) [11], the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [12], and the Parkes Pulsar Timing Array (PPTA) [13]. The collaboration of the three projects is called the International 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 domain using the NANOGrav 11-year Data Set. We quantitatively investigate whether the ultralight scalar dark matter is detectable or not using the Bayesian model selection approach. We clarify the prior dependence of constraints on the amplitude of the FDM and obtain three times stronger constraints on the amplitude of the FDM in the frequency range from $10^{-8.34}$ to $10^{-8.19}\,{\rm Hz}$. We also discuss how the results of Bayesian analysis depend on the the solar system ephemeris noise in the model describing the observation data.

This article is organized as follows. In Section II we describe a model of FDM signal. In Section III we briefly review the Bayesian statistics, and explain how to use it for our analysis. In Section IV we describe the model, the data, and the function used in the Bayesian analysis. In Section V we briefly review the MCMC simulation. In Section VI we describe the analysis of the white noise that is performed before the main analysis. In Section VII we summarize the results of Bayesian analysis using NANOGrav 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 gravitational potential and the oscillation of the arrival time of the pulse passing through the gravitational potential. The oscillation of the arrival time of the pulse induced by the FDM is given by [17]

[TABLE]

where $f=m/\pi$ is the frequency, $m$ is the mass, $D$ is the distance between the pulsar and the earth, each ${\bf x_{e}}$ and ${\bf x_{p}}$ are the position of the Earth and the pulsar, and $\alpha$ denotes the phase. Here, we used the gravitational potential

[TABLE]

where $\rho({\bf 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:

[TABLE]

where we assumed $\Psi\equiv\Psi({\bf x_{e}})=\Psi({\bf x_{p}})$ and defined

[TABLE]

Here, since we do not aim to estimate the distance, we put together the phase $\alpha({\bf x_{p}})$ and the distance $L$. In fact, the distance has an uncertainty of tens to hundreds of parsec currently that is too large to determine the phase [18]. Since, as is mentioned in the article [16], the distance between the Earth and the pulsar $L$ is not so large, it is reasonable to assume that the amplitudes at the earth $\Psi({\bf x_{e}})$ and the pulsar $\Psi({\bf x_{p}})$ are equal. The FDM signal is a simple sine wave superposition, and so some noise in the data may have a similar waveform. Therefore, it is necessary to consider whether the FDM signal obtained by data analysis is the actual FDM signal or some noise.

If all dark matter is occupied by the FDM, the predicted amplitude at the position of the Earth is

[TABLE]

and

[TABLE]

where $\rho=0.4{\rm Gev/cm^{3}}$ 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 Bayesian data analysis, see for example [20, 21, 22].

The purpose of the Bayesian parameter estimation is to estimate the posterior probability distribution $p(\theta|D)$ of the parameters $\theta$ given the data $D$. Having the data, we can update our belief about the parameters using Bayes’ rule, namely

[TABLE]

where $p(\theta)$ is the prior probability distribution, $p(D|\theta)$ is the likelihood function, and $p(D)$ is the evidence. In the parameter estimation, the evidence can be regarded as a normalization constant which can be omitted, because it does not involve the parameter. We will define the data $D$, the parameter $\theta$, the likelihood function $p(D|\theta)$, and the prior probability distribution $p(\theta)$ in Section IV.

Using the posterior probability distribution for the amplitude of the FDM signal $\Psi$, we define the 95% upper limit $R$ by

[TABLE]

We will calculate this 95% upper limit with the samples of the posterior probability distribution which is generated by the MCMC in Section VII.1.

For the Bayesian model comparison, it is often considered the Bayes factor:

[TABLE]

where $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are competing models which assign a meaning to the parameters. The Table I gives the interpretation of the Bayes factor in terms of the strength of the evidence [23, 24, 25].

If $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are nested models, we can use the Savege-Dickey density ratio to calculate the Bayes factor [26, 27, 28, 24]. We assume that $\mathcal{M}_{1}$ is a model described in Section IV.2, and $\mathcal{M}_{2}$ is a model in which the amplitude $\Psi$ of the model $\mathcal{M}_{1}$ is fixed to a very small value $10^{-18}$. In this case, the model $\mathcal{M}_{2}$ can be regarded as a model with no FDM signal. Then the Savege-Dickey density ratio is

[TABLE]

We will calculate this Bayes factor with the samples of the posterior probability distribution which is generated by the MCMC in Section IV. Specifically, we calculate the Bayes factor in multiple small bins around a fixed value $\Psi_{0}$, then derive the mean and unbiased standard deviation.

IV Bayesian analysis in the time domain

In this section, we explain the data $D$, the model $\mathcal{M}$, and the parameter $\theta$, and define the likelihood function $p(D|\theta)$, and the prior probability distribution $p(\theta)$.

IV.1 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) precision and long observation period, which would be suitable for detecting the FDM signal which becomes larger as the frequency becomes 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 parameters included in the timing model are spin parameters, astrometry parameters, binary parameters (if pulser is binary), dispersion measure parameters, frequency dependency parameters, jump parameters, see the articles [30, 18] for details. Note that, by fitting 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)111https://www.bipm.org/en/bipm-services/timescales/time-ftp/ttbipm.html is used for the Terrestrial Time, and JPL DE436 [33] is used for the planetary ephemeris. It is 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.2 was first introduced to the model by the article [34].

In order to obtain the timing residuals, we use the libstempo222http://vallis.github.io/libstempo which is the PYTHON interface to TEMPO2333https://bitbucket.org/psrsoft/tempo2.git [35] timing package.444We 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 the file ‘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, TEMPO2 was installed as expected. We used the identical data set except for the parameter file of PSR J1713+0747. In the parameter file of the PSR J1713+0747, we changed only a parameter EPHEM from DE430 [36] to DE436, where this parameter specifies which ephemerides to be used. Then we used libstempo to fit the timing parameters of the PSR J1713+0747 and created a new parameter file. We iterated the timing fitting five times, which would be sufficient for parameters to converge to certain values. We found that the LAMBDA and BETA parameters of the astrometry parameters shift by more than one sigma away from the older value. This result can be inferred from the deviation of the position parameter due to the difference in ephemeris as pointed out in the article [37].

IV.2 Model

Following the paper [30, 34], the model $\mathcal{M}_{1}$ of the timing residuals $\delta\bm{t}$ for each pulsar includes the FDM signal, the linearized timing 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 noise is represented by the product of a design matrix and small offsets. The design matrix describes the dependence of the timing residuals on respective timing model parameters. We obtain the design matrix using the TEMPO2 via libstempo. The timing model 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 components are thought to be caused by a random walk in one of the pulsar spin parameters [39, 40, 41, 42] and contributions to TOAs by an asteroid belt around the pulsar [43]. The chromatic components are thought to be caused by the pulse propagating through the ionized interstellar medium [44, 38]. Although the origins of red noise are various, simple power-law spectrum form is often used as the power spectral density. Specifically, the power spectral density of the red noise $P(f)$ is defined as

[TABLE]

where $f$ is a red noise frequency, $f_{\rm yr}$ is $1{\rm yr}^{-1}$, $A$ is a dimensionless amplitude of the red noise, and $\gamma$ is a spectral index of the red noise.

Third, the Solar System ephemeris noise characterizes the inaccuracy in the vector between the geocenter and the Solar System barycenter caused by the inaccuracies of the Solar System ephemeris. We will refer to this noise as the SSE noise. It is known that SSE noise affect upper limits and Bayes factors for amplitudes of the stochastic gravitational wave background [34]. The components of the SSE noise are considered to be the uncertainties of the planet mass and the planet orbit, and the uncertainty of the rotation rate around the ecliptic pole. The uncertainties of the planet mass are taken into account in Jupiter, Saturn, Uranus and Neptune. The uncertainties of the planet orbit are taken into account in Jupiter. We used the values and the data implemented in ENTERPRISE (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE) which is a pulsar timing analysis code555https://github.com/nanograv/enterprise. Thus, regarding the SSE noise caused by the uncertainties of the planet orbit, the Jupiter mass is the value of the IAU 2009 system of astronomical constants [45]. Also, the data describing the dependence of the Jupiter orbit on respective set-III parameters [46] is the same as ENTERPRISE. Regarding the SSE noise caused by the uncertainty of the rotation rate around the ecliptic pole, the reference date corresponds to MJD 55197.

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 additional white noise. The EFAC and the EQUAD are used independently for each combination of receivers and backend systems. The ECORR differs from the EQUAD in that there is a correlation between the TOAs obtained during a single observation. This ECORR characterizes pulse jitter caused by stochastic amplitude and phase variations in pulse, which correlates in a certain frequency band and doesn’t correlate in time [44].

IV.3 Likelihood Function

The likelihood for each pulsar can be written as [30, 47]:

[TABLE]

where

[TABLE]

$\delta\bm{t}$ is the $N_{\rm TOA}$ dimensional timing residuals, $N_{\rm TOA}$ denotes the number of TOAs of the pulsar, $\bm{s}$ is the $N_{\rm TOA}$ dimensional FDM signal, $\bm{n}_{\rm SSE}$ is the $N_{\rm TOA}$ dimensional linearized timing model noise, $\bm{M}$ is a $N_{\rm TOA}\times N_{\rm TM}$ design matrix, $N_{\rm TM}$ denotes the number of the timing model parameters, $\bm{F}$ is a $N_{\rm TOA}\times 2N_{\rm red}$ matrix which has columns of alternating cosine and sine functions for the red noise, $N_{\rm red}$ denotes the number of frequencies, $10^{80}$ are the variance of the TM noise parameters, $\Xi$ is the power spectrum of the corresponding frequency times the inverse of the observation period of the pulsar, and ${\bm{C}}_{\rm white}$ is the $N_{\rm TOA}\times N_{\rm TOA}$ covariance matrix which is composed by the EFAC, EQUAD, and ECORR. The likelihood function using multiple pulsars is represented by the product of the likelihood function of each pulsar.

IV.4 Prior Probability Distribution

We use specific knowledge only for the mass errors of each planet as the prior information. Using the propagation of uncertainty law, the variances of Jupiter mass, Saturn mass, Neptune mass are calculated from the IAU 2009 system of astronomical constants. Also, the variance of Uranus mass is calculated from the values in the article [48] which is newer than the IAU 2009 system of astronomical constants. Then we assume a normal distribution for the uncertainty of each planet mass and apply the obtained variances. For parameters without specific knowledge, we use a log-uniform distribution for parameters which are need to be searched over several orders of magnitude with only positive values, and we use a uniform distribution for the other parameters. The range of the log-uniform distribution and the uniform distribution is taken sufficiently wider than the value that the parameter would take. Regarding the amplitude of the FDM signal, we especially consider both cases of uniform distribution and log-uniform distribution as in the article [34]. The uniform distribution is used for the upper limit, and the log-uniform distribution is used for the model comparison. The parameters and their prior probability distribution are given in Table 2.

As is done in [49, 50, 34], we analyze the white noises in advance before the main analysis. The resulting MCMC chains are used to calculate the value that maximizes the one-dimensional posterior probability distribution corresponding to each white noise, where this value is called the maximum a posteriori (MAP) value. The main analysis is performed by fixing the possible values 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 we used is called parallel tempering. In the parallel tempering method, a concept of temperature is introduced, and the MCMC simulations of different temperatures are executed in parallel. The advantage of parallel tempering is that it is possible to 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=1.00,4.64,21.5,100$.

In order to perform the parallel tempering, we use the software package PTMCMCSampler888https://github.com/jellis18/PTMCMCSampler [51] via PAL2999https://github.com/jellis18/PAL2 [52] which is a Bayesian inference package for PTA and can include the PTMCMCSampler. Regarding models not implemented in the PAL2, the FDM signal is implemented like the continuous gravitational waves and the SSE noise is implemented like any other noises. Following the article [49], we use adaptive Metropolis [53], single component adaptive Metropolis [54], and differential evolution [55], as a proposal algorithm which is used to generate next samples using past samples. Furthermore, we also use a simple proposal algorithm to generate the next sample of each parameter by proposal distribution which is the same distribution as the probability distribution. All of these proposal algorithms are used in a single MCMC simulation and which one is used is chosen randomly for each proposal in the MCMC simulation. In this article, we use the value written in the PAL2 for each variable used in PTMCMCSampler, unless specifically mentioned.

VI Pre-analysis

As is usual [49, 50, 34], in order to obtain the MAP values of the parameters of the white noise, we analyze the white noise first before the main analysis. By doing this, the number of the free parameters can be reduced in the main analysis. In the pre-analysis, we performed independent analysis for each pulsar, and we used the model which contains the red noise in addition to white noise. We ran the MCMC simulation with $10^{6}$ iterations and removed the first 25% as a burn-in period, where the burn-in period is the period during which samples have not yet been obtained from the target distribution.

Table 3 shows the results of the white noise and the red noise obtained by the pre-analysis. The obtained posterior distribution is expressed by the MAP value and the 95% confidence interval. At the Green Bank Observatory, the receivers are Rcvr_800 and Rcvr1_2, and the backend system is GASP in early observations and GUPPI in later observations. Similarly, at the Arecibo Observatory, the receivers are L-wide and S-wide, and the backend system is ASP in early observations and PUPPI in later observations. 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 absorbed by 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.

VII.1 Upper limits

We calculated the 95% confidence upper limits on the amplitude of the FDM signal by the Bayesian analysis. We ran all the MCMC simulation with $10^{6}$ iterations and removed the first 25% as a burn-in period. As the prior probability distribution of the amplitude of the FDM signal, we considered two cases: the uniform distribution and the log-uniform distribution. The uniform prior was used to place the conservative upper limits, the log-uniform prior was used to calculate the Bayes factors, where the upper limits were calculated using the Eq. (8) and the Bayes factors were calculated using Eq. (10). In order to see the effect of including the SSE noise in the model on the results, we also calculated the upper limits and the Bayes factors when the SSE noise is not included in the model.

In Figure 1, we show the upper limits and the Bayes factors of the amplitude $\Psi$ as a function of the frequency $f$ and the mass $m$. The relation between the frequency $f$ and the mass $m$ is given by Eq. (6). First, in the above plot, the black solid and dashed lines denote the upper limit using uniform prior and log-uniform prior, respectively. Here, we plotted the results obtained using log-uniform prior, but as we mentioned in the section IV.4, we regard the results obtained by the uniform prior as conservative upper limits. The red solid and dashed lines denote the upper limit obtained when the SSE noise was not included in the model, and using uniform-prior and log-uniform prior, respectively. The bold black line denotes the upper limit of the Bayesian analysis obtained in [16] ( taken from Figure 3). The green line denotes the predicted amplitude of the FDM signal given by Eq. (5) with $\rho=0.4{\rm Gev/cm^{3}}$. The purple vertical lines denote the inverse of the observation periods of pulsars and corresponds 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 the SSE noise, respectively. The unbiased standard deviation is used for error bars. Only to make the plot easier to see, when the Bayes factor exceeds 20, it is represented by the upper triangle and the mean value and the unbiased standard deviations of the Bayes 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 limits for the log-uniform prior gives stronger limits than for the uniform prior, but the difference is small. The reason the difference is small is that the Bayes factor exceeds 3 when the frequency becomes $10^{-8.19}\,{\rm Hz}$ ($1.34\times 10^{-23}\,{\rm eV}$) or lower. According to the Table 1, the Bayes factor exceeds 3 means that there is a signal that is somewhat similar to the FDM signal. Therefore, whichever prior probability distribution is used, the value of the posterior probability distribution tends to be large at the parameter values of that signal, and as a result the upper limits does not change so much. From the above, it was found that the FDM signal is positively supported in the wide frequency range. However, from a physical point of view, it is hard to think that FDM signals has been found, because the upper limit obtained is about an order of magnitude greater than the expected amplitude and also FDM 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 of the FDM signal obtained using uniform distribution does not change much. However, the difference between the upper limits obtained from the different prior probability distributions is large. The reason for this is that all Bayes factors are smaller than 3 and in most cases they do not exceed 1. Therefore, in this case we conclude that no FDM signal has been detected. The red and black results show that the Bayes factor is smaller when the SSE noise is considered. Thus, we conclude that the SSE noise can mimic the FDM signal. This result would be inferred from the result that the SSE noise can mimic the stochastic gravitational wave background [34]. In comparison with the published Bayesian upper limits of the amplitude of the FDM signal using the PPTA 12-year data set [16], i.e. comparing the black and the bold black lines, we found that stronger upper limits were obtained when the frequency was in the range from $10^{-8.34}$ to $10^{-8.19}{\rm Hz}$ (from $9.45\times 10^{-24}$ to $1.34\times 10^{-23}\,{\rm eV}$). In this range, up to three times 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 amplitude of the FDM signal. Thus, we convert the amplitude of the FDM signal into the energy density using Eq. (5), and the result is 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 the range from $10^{-8.55}$ to $10^{-8.01}\,{\rm Hz}$ (from $5.83\times 10^{-24}$ to $2.02\times 10^{-23}\,{\rm eV}$) where we analyzed. The strongest upper limit on the the energy density is $2\,{\rm GeV/cm^{3}}$ at the frequency $10^{-8.28}\,{\rm Hz}$ ($1.09\times 10^{-23}\,{\rm eV}$).

VII.2 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 using the obtained MAP values of the parameters. We ran the MCMC simulation with $10^{6}$ iterations for analysis of red noise and SSE noise and with $10^{5}$ iterations for analysis of the FDM signal, and in both cases we removed the first 25% as a burn-in period. The result of the red noise and the SSE noise used for the fixed noise analysis are summarized in Table 4. The obtained posterior distribution 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 distribution of the amplitude of the FDM signal. The results are plotted in Figure 3 which is the similar plot as Figure 1. The solid and dashed 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 FDM signal. 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 analyzed in advance. We consider that this analysis is not suitable for giving an upper limit to the amplitude of the FDM signal, and the obtained 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 previous section. In particular, when using the log-uniform prior, surprisingly, the upper limits are smaller than the predicted amplitude in some range. As for Bayesian factors, they are all smaller than 1, which is consistent with the fact that the upper limits are strongly influenced by the prior probability distribution. From this result, it is inferred that the FDM signal is well absorbed by the 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 noise Then we calculated the Lomb-Scargle periodogram of the timing residual by subtracting it, where the Lomb-Scargle periodogram can be used to search for periodic signals in non-uniformly spaced time series data [56, 57]. The reason for not subtracting the red noise from the original timing residual is that it is difficult to create the noise included in the actual data because the red noise that we can create is only one realization of the stochastic process. For comparison, we also calculated the Lomb-Scargle periodogram of the original timing residual and the simulated timing residual induced by the red noise only. In order to calculate the Lomb-Scargle periodogram, we used Astropy121212http://www.astropy.org [58, 59] which is a Python package for the astronomy. For the purpose of expressing red noise, the original residual has short observation period and lacks frequency resolution, so that we made red noise using simulated observation data that the observation period is $10^{5}$ 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 the Lomb-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 red noise is also displayed below it, because the red noise is large and the other lines are difficult to see. The purple vertical line represents 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 in the frequency smaller than $10^{-8}{\rm Hz}$ is larger than the power in the frequency larger than $10^{-8}{\rm Hz}$. Therefore, it can be understood that 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 to characterize low frequency noise in PSRs J0613-0200 and J1012+5307, but it can not be considered so in the other pulsars. In particular, the amplitude of red noise is large in PSR J1744-1134, so that in the analysis in which the MAP value of the red noise is fixed, it can be seen that this pulsar will contribute little to the analysis result of the amplitude of the FDM signal. The reason why the red noise can be detected properly in PSRs J0613-0200 and J1012+5307 is that, according to the result of analyzing only 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 that the low frequency noise are reduced in PSRs J0613-0200 and J1909-3744. On the other hand, this is not the case with the other pulsars, and it can be seen that the noise increases rather than decreases. We conclude that the reason why the upper limits are smaller in fixed noise analysis than in the main analysis is mainly due to PSR J1909-3744, because only this pulsar has small red 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 is a problem that the model of the SSE noise may favor parameters which are contrary to the expectation that noise will be reduced on all pulsars. More precise measurements of the SSE would be needed to detect or exclude the FDM signal. However, these issues are beyond the scope of this article.

VIII Conclusion

We searched for the FDM by performing the Bayesian analysis in the time domain using the NANOGrav 11-year Data Set. In Section VII.1, we gave the 95% confidence upper limit on the amplitude of the FDM signal. we found that probability that the FDM signal should be included in the model was less than 75% in all frequency region. Compared with the published Bayesian upper 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 than the previous study when the frequency was in the range from $10^{-8.34}$ to $10^{-8.19}{\rm Hz}$ ($9.45\times 10^{-24}$ to $1.34\times 10^{-23}\,{\rm eV}$ in terms of the FDM mass). In other region, we also obtained the similar upper limit on the amplitude of the FDM signal. Since the amplitude of FDM signal can be converted to the energy density of the dark matter near the Earth, it is easy to obtain the upper limit of the energy density. The upper limit on the energy density was lower than $7\,{\rm GeV/cm^{3}}$ in the range from $10^{-8.55}$ to $10^{-8.01}\,{\rm Hz}$ (from $5.83\times 10^{-24}$ to $2.02\times 10^{-23}\,{\rm eV}$) where we analyze. In particular, at a frequency of $10^{-8.28}\,{\rm Hz}$ (a mass of $1.09\times 10^{-23}\,{\rm eV}$), we obtained the strongest upper limit $2\,{\rm GeV/cm^{3}}$. In addition to the main analysis, we also investigated the case where the SSE noise was not included in the model. In this case, we showed that we can not exclude the existence of the FDM signal, because the probability that the FDM signal should be included in the model was more than 75% in the frequency region $10^{-8.19}\,{\rm Hz}$ or less. This results show that the Bayes factor is smaller when the SSE noise is considered. Thus, we conclude that the SSE noise can mimic the FDM signal, which would be inferred from the result that the SSE noise can mimic the stochastic gravitational wave background [34].

In Section VII.2, by analyzing the noise in advance, we examined how much the FDM signal was absorbed. In this case, we clarified that the probability that the FDM signal should be included in the model was much lower than 50% in all frequency region. 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 a simulated noise from the parameters obtained by the analysis earlier, and investigated whether removing the SSE noise from the observed data would reduce the power of the low frequency region of the data. We found that the power of the PSR J1909-3744 has become smaller, and we conclude that this pulsar contributed to the improvement of the sensitivity. With other pulsars, we also found that the power increased on the contrary. If the SSE noise is properly found, it is thought that the noise will be reduced for all pulsars. Hence, this result seems strange. In order to solve this problem, it is necessary to reduce the uncertainty of the SSE by more accurate observation of the SSE.

Acknowledgements.

R.K. was supported by Grant-in-Aid for JSPS Research Fellow and JSPS KAKENHI Grant Numbers 17J00496. J. S. was in part supported by JSPS KAKENHI Grant Numbers JP17H02894, JP17K18778, JP15H05895, JP17H06359, JP18H04589. J. S is also supported by JSPS Bilateral Joint Research Projects (JSPS-NRF collaboration) “String Axion Cosmology.”