# Parkes Pulsar Timing Array constraints on ultralight scalar-field dark matter

###### Abstract

It is widely accepted that dark matter contributes about a quarter of the critical mass-energy density in our Universe. The nature of dark matter is currently unknown, with the mass of possible constituents spanning nearly one hundred orders of magnitude. The ultralight scalar field dark matter, consisting of extremely light bosons with eV and often called “fuzzy” dark matter, provides intriguing solutions to some challenges at sub-Galactic scales for the standard cold dark matter model. As shown by Khmelnitsky and Rubakov, such a scalar field in the Galaxy would produce an oscillating gravitational potential with nanohertz frequencies, resulting in periodic variations in the times of arrival of radio pulses from pulsars. The Parkes Pulsar Timing Array (PPTA) has been monitoring 20 millisecond pulsars at two to three weeks intervals for more than a decade. In addition to the detection of nanohertz gravitational waves, PPTA offers the opportunity for direct searches for fuzzy dark matter in an astrophysically feasible range of masses. We analyze the latest PPTA data set which includes timing observations for 26 pulsars made between 2004 and 2016. We perform a search in this data set for evidence of ultralight dark matter in the Galaxy using Bayesian and Frequentist methods. No statistically significant detection has been made. We therefore place upper limits on the local dark matter density. Our limits, improving on previous searches by a factor of two to five, constrain the dark matter density of ultralight bosons with eV to be below with 95% confidence in the Earth neighborhood. Finally, we discuss the prospect of probing the astrophysically favored mass range eV with next-generation pulsar timing facilities.

^{†}

^{†}preprint: APS/123-QED

## I Introduction

Dark matter, a concept established in the early 1930s for the purpose of explaining the observed enigmatic dynamics of disk galaxies and motion of galaxies in clusters Zwicky (1933, 1937); Smith (1936), is nowadays considered to be an essential ingredient of the Universe. It is instrumental in explaining a wide range of astrophysical phenomena, such as strong gravitational lensing of elliptical galaxies Koopmans and Treu (2003), the dynamics of interacting clusters Clowe et al. (2004) and the large-scale structure of the Universe Tegmark et al. (2004). The latest analysis of temperature and polarization anisotropies of the cosmic microwave background Planck Collaboration et al. (2016) suggested that the Universe contains 26% dark matter, which is five times more than ordinary baryonic matter such as stars and galaxies.

The most popular dark matter candidates are weakly interacting massive particles (WIMPs) and QCD (quantum chromodynamics) axions. We refer to both as standard cold dark matter, or simply CDM. The CDM paradigm has met with impressive success in matching observational data on large cosmological scales (see Bertone et al. (2005); Primack (2012), for reviews). Recently, there has been an increased number of ideas about dark matter that go beyond the standard paradigm, building on old ideas in some cases (see e.g. Battaglieri et al. (2017) for an overview).

One such idea–an ultralight axion or axion-like particle–can be thought of as a generalization of the QCD axion. An axion is an angular field, i.e., the field range is finite and periodic with a periodicity with often referred to as the axion decay constant. A simple axion Lagrangian has a standard kinetic term, and a self-interaction potential generated by non-perturbative effects (that can be approximated by instanton potential):

(1) |

where is the mass of the axion .
The non-perturbative effects are typically highly suppressed
(e.g. exponentially suppressed by an instanton action), leading to
a fairly low energy scale .
In the early Universe, the scalar field is frozen at its primordial
value, generically expected to be order of . When the
Hubble expansion rate drops below the mass scale , the scalar field
oscillates with an amplitude that redshifts with the expansion of the
Universe. Averaging over oscillation cycles, behaves like CDM
with a relic density of (see e.g. Arvanitaki et al. (2010); Hui et al. (2017))^{1}^{1}1The relic density computation follows the classic arguments of
Preskill et al. (1983); Abbott and Sikivie (1983); Dine and Fischler (1983), which were developed
for the QCD axion.

(2) |

String theory contains many axion candidates with somewhere in the range
GeV Svrcek and Witten (2006).
Eq. (2) tells us that a very low is preferred
if the axion were to account for dark matter.
It should be emphasized though that there is a fairly large possible
range for ; in fact, the relic abundance is more sensitive
to than to .
A lighter mass, e.g., eV, can be easily accommodated by a slightly higher , though it is disfavored by astrophysical observations such as the existence and structure of dwarf galaxies^{2}^{2}2Note that the requisite is much less than the QCD scale; hence this is not the QCD axion..

Such an ultralight axion has a macroscopic de Broglie wavelength , given by

(3) |

where is the velocity, implying wave-like phenomena on astronomically accessible scales, unlike standard CDM. In linear perturbation theory, the wave-like property leads to a suppression of power on small scales (small compared to the Jeans scale, which is a geometric mean of the Compton and Hubble scale). It is this property that motivated Hu, Barkana and Gruzinov Hu et al. (2000) to propose an ultralight boson as an alternative to standard CDM, and to coin the term “fuzzy dark matter” (FDM). The term FDM refers generally to a scalar dark matter particle with a very small mass, such that its de Broglie wavelength is macroscopic. An ultralight axion is a particularly compelling realization. Our constraints derived in this paper apply to the ultralight axion, as well as the broader class of FDM.

The thinking was that the suppression of power on small scales would help resolve certain small-scale problems of CDM, which generally have to do with CDM predicting too much small-scale structure compared to that observed. There is a vast literature on this subject, but it remains a matter of debate as to whether the perceived small-scale structure problems of CDM are in fact amenable to astrophysical solutions (such as feedback processes modifying the mass distribution within Galactic halos); see Bullock and Boylan-Kolchin (2017) for a review.

There exist several different bounds on the FDM model. One class of bounds comes from measurements of the linear power spectrum at high redshifts, such as from the microwave background (e.g. Hložek et al. (2018)), and from the Lyman-alpha forest Iršič et al. (2017); Kobayashi et al. (2017). In particular, the Lyman-alpha forest data appear to disfavor a FDM mass lighter than about eV. Another example of a bound of this kind come from 21-cm observations – the recent detection of a global 21-cm absorption signal at redshift around Bowman et al. (2018) puts a lower limit on the FDM mass similar to the Lyman-alpha forest bound Schneider (2018); Lidz and Hui (2018); Sullivan et al. (2018). Yet another class of bounds comes from dynamical data on the density profiles of galaxies e.g. Calabrese and Spergel (2016); Deng et al. (2018); Bar et al. (2018). Many of these bounds are subject to their own astrophysical uncertainties. For instance, the Lyman-alpha forest bound is predicated upon the correct modeling of fluctuations from such as the ionizing background, the temperature and feedback processes. The 21-cm bound relies on assumptions about star formation (that it tracks the halo formation and that the fraction of baryons that form stars is less than about ), and of course assumes the validity of the detection. Constraints from rotation curve measurements generally make assumptions about how feedback processes, such as from stellar explosions, affect (or do not affect) density profiles.

Recently, a number of authors, based on numerical simulations and analytical arguments, pointed out additional testable astrophysical implications of FDM, especially in the nonlinear regime Schive et al. (2014); Mocz and Succi (2015); Veltmaat and Niemeyer (2016); Hui et al. (2017); Nori and Baldi (2018); Veltmaat et al. (2018). A particularly interesting probe of ultralight dark matter using pulsar timing arrays (PTAs) was pointed out by Khmelnitsky and Rubakov Khmelnitsky and Rubakov (2014). Through purely gravitational coupling, scalar field dark matter induces periodic oscillations in gravitational potentials with frequency twice the field mass . The oscillating gravitational potentials along the line of sight of pulsars cause sinusoidal variations in the times of arrival (ToAs) of radio pulses. The frequency of such variations lies right in the sensitivity band of PTAs. This way of detecting or constraining FDM is completely independent of other methods (and their assumptions), and provides a useful check. As shown in Khmelnitsky and Rubakov (2014); Porayko and Postnov (2014); De Martino et al. (2017); Blas et al. (2017) and later in this paper, the current PTA data can only be sensitive to very low-mass FDM ( eV). We will discuss what would be required to probe the higher and cosmologically more favorable masses.

The concept of a PTA is to regularly monitor ToAs of pulses from an array of the most rotationally stable millisecond pulsars Sazhin (1978); Detweiler (1979); Hellings and Downs (1983); Foster and Backer (1990). Measured ToAs are fitted with a deterministic timing model that accounts for the pulsar spin behavior and for the geometrical effects due to the motion of the pulsar and the Earth. The difference between the observed ToAs and those predicted by the best-fit timing model are called “timing residuals”. By analyzing the pulsar timing residuals, we can obtain the information about other physical processes that affect the propagation of radio pulses through the Galaxy, for instance, the presence of ultralight scalar field dark matter in the Galaxy.

The Parkes Pulsar Timing Array (PPTA) Manchester et al. (2013) uses the 64-m Parkes radio telescope in Australia. Building on earlier pulsar timing observations at Parkes, it started in 2005 to time 20 millisecond pulsars at a regular interval of two to three weeks. PPTA and its counterparts in North America (NANOGrav) McLaughlin (2013) and Europe (EPTA) Kramer and Champion (2013) have joined together to form the International Pulsar Timing Array (IPTA) Hobbs et al. (2010a); Verbiest et al. (2016), aiming for a more sensitive dataset. The IPTA currently observes around 70 pulsars using the world’s most powerful radio telescopes.

The first PPTA data release was published in 2013 Manchester et al. (2013). It included six years of observations for 20 pulsars. This data set was used to search for a stochastic gravitational wave (GW) background Shannon et al. (2013), continuous GWs Zhu et al. (2014) and GW bursts with memory Wang et al. (2015a). The second data release is still being actively developed, but for this paper, we have made use of a data set that contains observations made between 2004 and 2016 with five new pulsars added since 2010. An early subset of this data was used to place the most constraining limit to date on the amplitude of a stochastic GW background in the nHz regime Shannon et al. (2015).

In this work we search for evidence of ultralight scalar field dark matter in the Galaxy using the PPTA data. A similar study was carried out, through Bayesian analysis, by Porayko & Postnov Porayko and Postnov (2014) using the NANOGrav 5-yr 17-pulsar data set published in Demorest et al. (2013). Our work improves on that of Porayko and Postnov (2014) in several ways. First, we make use of an independent data set with much longer data span and smaller errors in the timing residuals. Second, we use an up-to-date Bayesian inference packages for PTA data analysis–PAL2 Ellis and van Haasteren (2017) and NX01 Taylor and Baker (2017)–and include proper treatment of the noise processes. Re-analyzing the NANOGrav data with the improved analysis, we find that the sensitivity was overestimated by a factor of ten in Porayko and Postnov (2014). Third, we also adopt a standard Frequentist searching method and obtain consistent results with Bayesian analysis.

Our paper is organized as follows. In Sec. II we describe pulsar timing residuals expected in the presence of ultralight scalar field dark matter in the Galaxy. In Sec. III we introduce our data set, the likelihood function and our Bayesian and Frequentist methods to model the noise properties of PPTA data. We also present results of our noise analysis. In Sec. IV we describe our search techniques and apply them to the PPTA data set. As we find no significant signals, we set upper limits on the local density of FDM in the Galaxy. In Sec. V we discuss how the sensitivity will be improved in the future. Finally, we provide concluding remarks in Sec. VI.

## Ii The pulsar timing residuals from fuzzy dark matter

In this section we briefly describe the magnitude and time dependence of timing residuals induced by the scalar field dark matter in the Galaxy. A detailed derivation can be found in Khmelnitsky and Rubakov (2014).

Because of the huge occupation number, the collection of ultralight dark matter particles behaves like a classical scalar field .
To a very good approximation, here we ignore quartic self-interaction and coupling of ultralight dark matter particles to other fields^{3}^{3}3In the axion context, the oscillation amplitude of gradually diminishes due to the expansion of the universe, making the quadratic an excellent approximation to the potential in Eq. (1). (Brdar et al., 2018).
The scalar action in this case can be written as:

(4) |

to which the standard Einstein-Hilbert action for the metric should be added. The equation of motion is the Klein-Grodon-Fock Equation: . We are interested in a computation of and the metric inside the Galaxy. The metric is approximately Minkowski plus corrections at the level of . To good approximation, everywhere in the Galaxy oscillates at an angular frequency (corrections due to the momentum of the particles and the gravitational potential are small). The energy-momentum tensor to the leading order diagonal and its spatial components (pressure) oscillate at twice the field particle mass. This produces time-dependent gravitational potentials and in the metric tensor (in the Newtonian covariant form) with leading oscillating contributions at a frequency

(5) |

The amplitude of oscillating parts of the potentials and are a factor of smaller than the time-independent parts , where is the local scalar field dark matter density. For cosmologically favored boson masses eV, the frequency is fortuitously located in the sensitivity range of PTAs.

As in the case of GWs Detweiler (1979), pulsar photons propagating in a time-dependent metric undergo a frequency shift , which is related to timing residuals Khmelnitsky and Rubakov (2014)

(6) |

where is the distance to the pulsar and is the amplitude of cosine component of the oscillating part of the energy-momentum tensor. The subsequent terms in Eq. (6) are suppressed with respect to by a factor , and to the leading order the signal does not depend on the oscillating part of the potential .

As one can see in Eq. (6), the dark matter signal also has “Earth” and “pulsar” terms. Oscillation frequencies at the Earth and at the pulsar are identical, which makes it analogous to the case of non-evolving continuous GWs Zhu et al. (2016). The scalar-field oscillation phases on the Earth and pulsar generally take different values; but they become correlated when the Earth and a pulsar are located within the coherence de Broglie wavelength .

The amplitude , which can be effectively probed in pulsar timing experiments, depends on the local density of dark matter :

(7) |

where is the measured local dark matter density Bovy and Tremaine (2012); Read (2014); Sivertsson et al. (2018a, b). The root-mean-square (rms) amplitude of induced pulsar-timing residuals is

(8) |

The expected signal amplitude scales strongly with the boson mass. At eV and above, the signal is negligibly small. For mass below eV, the induced rms residuals ( ns) are comparable to current timing precision for the best pulsars, as we discuss in Sec. III.1.

In this work, we assume the Earth term and pulsar terms have the same amplitude . This is a reasonable approximation since most PPTA pulsars are relatively close ( kpc) to the Earth (see Table 1). We discuss effects of the dark matter density variability in Sec. V. Under this assumption, Eq. (6) can be written into a more compact form:

(9) |

where we have defined and with . Defining this way allows us searching for a single phase parameter per pulsar. One should note, however, that the parameter pair is indistinguishable from .

## Iii PPTA data and noise properties

### iii.1 Observations and timing analysis

Here we provide a brief overview of the data set used in this work. The data set is available from the CSIRO pulsar data archive^{4}^{4}4https://doi.org/10.4225/08/5afff8174e9b3.
The observing systems and data processing techniques are similar to the first data release (DR1) as described in Ref. Manchester et al. (2013). Table 1 summarizes key characteristics of the PPTA data set, including the median ToA uncertainties, weighted rms values of timing residuals, data spans and the number of observations.

Our data set consists of observations for 26 pulsars collected between 2004, February 5 and January 31, 2016 using the Parkes telescope. It includes DR1 data that were acquired between March 2005 and March 2011 for 20 pulsars, along with some earlier data for some pulsars that were observed in the 20-cm observing band prior to the official start of the PPTA project. Currently, the PPTA observes 25 pulsars, with PSR J17325049 having been removed from the pulsar sample in 2011 because we were unable to obtain high quality data sets for this pulsar.
The observing cadence is normally once every two to three weeks. In each session, every pulsar was observed in three radio bands (10, 20 and 50 cm) with a typical integration time of one hour. Twenty of these pulsars were monitored for more than ten years up to twelve years; only five pulsars have data spans less than five years. For this data set, the median ToA uncertainties vary from 149 ns (PSR J04374715) to 4.67 s (PSR J21243358); the weighted rms residuals in this data set vary from 152 ns (PSR J04374715) to 16.53 s (PSR J18242452A). PSRs J1939+2134 and J1824-2452A were excluded from the search analysis, as they show strong evidence for a large unmodeled red-noise component^{5}^{5}5This is evident as their rms residuals are much larger than the median ToA uncertainties given in Table 1. This may be accounted for using system- and band-specific noise terms (Lentati et al., 2016).. For our purpose, we find these two pulsars make little contribution to the sensitivity.

Pulsar Name | rms | Range | ||||
---|---|---|---|---|---|---|

(s) | (s) | (yr) | (kpc) | |||

J04374715 | 0.15 | 0.15 | 11.98 | 2004.022016.01 | 3820 | 0.16 |

J06130200 | 1.20 | 1.38 | 11.98 | 2004.022016.01 | 969 | 0.78 |

J07116830 | 3.29 | 1.58 | 11.98 | 2004.022016.01 | 1017 | 0.11 |

J10177156 | 0.97 | 0.76 | 5.54 | 2010.072016.01 | 524 | 0.26 |

J1022+1001 | 2.23 | 2.11 | 11.98 | 2004.022016.01 | 1008 | 1.13 |

J10240719 | 3.39 | 3.61 | 11.87 | 2004.022015.12 | 679 | 1.22 |

J10454509 | 3.82 | 3.35 | 11.98 | 2004.022016.01 | 854 | 0.34 |

J11256014 | 1.59 | 1.29 | 10.12 | 2005.122016.01 | 203 | 0.99 |

J14464701 | 1.81 | 1.47 | 5.19 | 2010.112016.01 | 161 | 1.57 |

J15454550 | 1.08 | 1.01 | 4.74 | 2011.052016.01 | 215 | 2.25 |

J16003053 | 0.91 | 0.71 | 11.98 | 2004.022016.01 | 969 | 1.80 |

J16037202 | 2.13 | 1.43 | 11.98 | 2004.022016.01 | 747 | 0.53 |

J16431224 | 1.75 | 2.96 | 11.98 | 2004.022016.01 | 713 | 0.74 |

J1713+0747 | 0.38 | 0.24 | 11.98 | 2004.022016.01 | 880 | 1.18 |

J17302304 | 2.01 | 1.48 | 11.98 | 2004.022016.02 | 655 | 0.62 |

J17325049 | 2.55 | 2.75 | 7.23 | 2004.032011.12 | 144 | 1.87 |

J17441134 | 0.68 | 0.61 | 11.98 | 2004.022016.01 | 855 | 0.40 |

J18242452A | 2.67 | 16.5 | 10.36 | 2005.052015.10 | 339 | 5.50 |

J18320836 | 0.53 | 0.25 | 2.86 | 2012.112015.10 | 68 | 0.81 |

J1857+0943 | 2.00 | 1.93 | 11.98 | 2004.022016.01 | 580 | 1.20 |

J19093744 | 0.25 | 0.16 | 11.98 | 2004.022016.01 | 1670 | 1.14 |

J1939+2134 | 0.36 | 1.43 | 11.87 | 2004.032016.01 | 591 | 3.50 |

J21243358 | 4.67 | 2.52 | 11.98 | 2004.022016.01 | 889 | 0.41 |

J21295721 | 1.82 | 1.19 | 11.65 | 2004.062016.01 | 540 | 3.20 |

J21450750 | 1.71 | 1.16 | 11.86 | 2004.032016.01 | 881 | 0.53 |

J22415236 | 0.44 | 0.28 | 5.98 | 2010.022016.01 | 615 | 0.96 |

During pulsar timing observations, ToAs are first referred to a local hydrogen maser frequency standard at the observatory. These ToAs are further transformed to Coordinated Universal Time (UTC) and then to a Terrestrial Time (TT) as published by the Bureau International des Poids et Mesures. For the current data set, we used TT(BIPM2015) and adopted the JPL DE418 Folkner et al. (2007) solar system ephemeris (SSE) model to project ToAs to the solar-system barycenter. Potential errors in SSE models are accounted for in our Bayesian analysis (Sec. IV.1).

Before performing the search for dark matter signals, we fit pulsar ToAs with a timing model using the standard TEMPO2 software package Hobbs et al. (2006); Edwards et al. (2006). Typical parameters included in this fit are the pulsar sky location (RAJ and DecJ), spin frequency and spin-down rate, dispersion measure, proper motion, parallax and (when applicable) binary orbital parameters. Additionally, constant offsets or jumps were fitted among ToAs collected with different receiver/backend systems. Below we describe our methods to estimate the noise properties of the PPTA data.

### iii.2 The likelihood function

The likelihood function for pulsar timing residuals, marginalized over the timing model parameters, can be written as van Haasteren et al. (2009); van Haasteren and Levin (2013):

(10) |

where is a vector of timing residuals with length , is the deterministic signal vector, including the dark matter signal as described in Sec. II and deterministic systematics, is the design matrix or regression matrix of the linear model Press et al. (1996) that describes how ToAs depend on timing model parameters^{6}^{6}6It can be obtained with the TEMPO2 designmatrix plugin..
The noise covariance matrix includes contributions from uncorrelated white noise (), time-correlated spin noise () and dispersion measure variations ().
In Eq. (10), we have defined . The covariance matrix depends on the set of noise parameters , and denotes deterministic signal parameters so that . As a result, this form of the likelihood, which was firstly implemented in van Haasteren et al. (2009), depends both on and , and provides the possibility of proper treatment of the signal in the presence of correlated noise and systematics.
The likelihood in the Eq. (10) can be further reduced to a more compact form (see Ref. van Haasteren and Levin (2013) for details):

(11) |

where the dimension matrix is obtained through the singular-value decomposition of the design matrix . Specifically, where and are unitary matrices with and dimension respectively, and is an diagonal matrix containing singular values of . The matrix is obtained such that with and consisting of the first and the remaining columns of respectively.

In this work, we assume that only the dark matter signal, noise processes (that will be described in the next subsection) and deterministic systematics, associated with SSE errors, contribute to the data. We neglect errors in terrestrial time standards and other common noise processes (such as a stochastic GW background). Therefore, the likelihood function for the full PTA can be expressed as a product:

(12) |

where is the number of pulsars in the timing array.

### iii.3 Noise modeling

For each pulsar in the PPTA data set, we estimate its noise properties using both Bayesian and Frequentist approaches. We present a general description of possible noise sources here.

Stochastic noise processes can be divided into the time-correlated and uncorrelated components. The uncorrelated (white) noise is represented by the uncertainties of the measured ToAs derived through cross-correlation of the pulsar template and the integrated profile. However, it is common that ToA uncertainties underestimate the white noise present in pulsar timing data. This might be caused by, e.g. radio frequency interference, pulse profile changes or instrumental artifacts. Two parameters, namely, EFAC (Error FACtor) and EQUAD (Error added in QUADrature), are included to account for excess white noise. They are introduced for each observing system used in the dataset. Following standard conventions, different parameterizations are used for EFAC and EQUAD. In TEMPO2 and for our Frequentist analysis, the re-scaled ToA uncertainties () are related to their original values () by

(13) |

In Bayesian analysis, we use the following relation

(14) |

Numerous studies Boynton et al. (1972); Blandford et al. (1984); Hobbs et al. (2010b) have found evidence for additional low-frequency noise in pulsar timing data. This time-correlated stochastic process is dominated by two components: achromatic (i.e, independent of radio frequency) spin noise and chromatic (i.e, dependent on radio frequency) such as dispersion measure (DM) variations. The former is intrinsic to the pulsar and might be related to pulsar rotational instabilities. The latter is associated with the interstellar medium which introduces time delays in pulsar ToAs. As pulsar travels in the tangent plane, the line of sight intersects spatially variable interstellar medium characterized by different column electron densities. For current receivers, the bandpass is generally not broad enough to resolve these kind of variations in each individual observation. Therefore, a typical strategy is to observe pulsars at widely separated radio bands, allowing the correction of DM variations.

Below we discuss details of noise modeling in the Bayesian and Frequentist frameworks.

#### iii.3.1 Bayesian framework

The Bayesian framework provides a consistent approach to the estimation of a set of parameters by updating the initial distribution of those parameters as more information becomes available. Bayes’ theorem states:

(15) |

where stands for the posterior (or updated) distribution of the parameters , given the data (or external information) , is the likelihood function, and is known as Bayesian evidence and defined as:

(16) |

The Bayesian evidence is a normalizing factor for parameter estimation problem and is a key criterion for model selection and decision making. Here does not depend on and it holds that . When applied for the case of PTAs, data includes an array of pulsar timing ToAs , includes and the likelihood is given by Eq. (10). The set of parameters, used for the Bayesian analysis, and the corresponding priors are described in Table 2.

Parameter | Description | Prior | Comments |
---|---|---|---|

Noise parameters () | |||

EFAC | White-noise modifier per backend | U[0, 10] | fixed for setting limits |

EQUAD | Additive white noise per backend | log-U[10, 4] | fixed for setting limits |

Spin-noise amplitude | log-U[20, 11] (search) | one parameter per pulsar | |

U[, ] (limit) | |||

Spin-noise spectral index | U[0, 7] | one parameter per pulsar | |

DM-noise amplitude | log-U[20, 6.5] (search) | one parameter per pulsar | |

U[, ] (limit) | |||

DM-noise spectral index | U[0, 7] | one parameter per pulsar | |

Signal parameters () | |||

Oscillation amplitude | log-U[20, 12] (search) | one parameter per PTA | |

U[, ] (limit) | |||

Oscillation phase on Earth | U[0, 2] | one parameter per PTA | |

U[0, 2] | one parameter per pulsar | ||

(Hz) | Oscillation frequency | log-U[9, 7] | delta function for setting limits |

BayesEphem parameters () | |||

Drift-rate of Earth’s orbit about ecliptic z-axis | U[, ] rad yr | one parameter per PTA | |

Perturbation of Jupiter’s mass | one parameter per PTA | ||

Perturbation of Saturn’s mass | one parameter per PTA | ||

Perturbation of Uranus’ mass | one parameter per PTA | ||

Perturbation of Neptune’s mass | one parameter per PTA | ||

Principal components of Jupiter’s orbit | U[0.05, 0.05] | six parameters per PTA |

For computational purposes, the noise covariance matrix from Eq. (10) can be split as a sum of a diagonal matrix and a large dense matrix , where is the diagonal matrix (), , where is the number of terms in the approximation sum. By using the Woodbury matrix lemma^{7}^{7}7 William (1989), the computationally heavy inversion of covariance matrix , involving operations, is reduced to lower rank diagonal matrix inversion . More details on this technique can be found in van Haasteren and Vallisneri (2015), Arzoumanian et al. (2014).

In this work we have used the so-called “Fourier-sum” prescription (or “time-frequency” method), introduced in Lentati et al. (2013). In this case, the Fourier transform matrix for achromatic processes can be written as:

(17) |

where , where is the whole timespan of the PPTA dataset, 11.98 years. The dimensionality of the Fourier matrix is (), where is number of frequency components, which in our case is 30. The noise vector for a specific noise process can be expressed as , where is the vector of Fourier coefficients.

The covariance matrix of Fourier coefficients can be derived from the covariance matrix of the theoretical power spectrum of a specific type of noise. Within Bayesian framework, we use the following parametrization for power-law noise:

(18) |

Therefore, the elements of the matrix , which are identical for both spin and DM noises, are expressed as:

(19) |

where iterates over different Fourier frequencies and is a Kronecker delta. If multiband observations are available, the degeneracy between the spin noise and DM contributions can be broken, because of the dependency of the amplitude of the DM variations on the observational frequency . This dependency enters in the Fourier transform matrix as:

(20) |

where Hzcmpc s and is the radio observing frequency at time . Using this terminology, the time delay between signal received at radio frequency and one received at is given by . Note that the linear and quadratic trends in DM variations get absorbed by timing model parameters DM1 and DM2, which are included in the Bayesian timing model. The inclusion of the DM derivatives in our analysis absolves us from the spectral leakage problem Lentati et al. (2014).

The formalism, described in this subsection, was implemented in a range of publicly available codes. For the single pulsar analysis we have used PAL2 Software – a package for the Bayesian processing of the pulsar timing data. Efficient sampling from the posteriors is performed by the Bayesian inference tool MULTINEST Feroz et al. (2009), running in constant efficiency mode – a computational technique that allows one to maintain the user-defined sampling efficiency for high-dimensional problems (see Ref. Feroz et al. (2013) for more details). For each PPTA pulsar we perform separately a full noise modeling analysis, simultaneously including all stochastic components discussed above. The noise parameters , estimated within single pulsar analysis, are given in Table 3. The marginalized posterior probabilities for the six most sensitive pulsars in PPTA (see Sec. IV.1) are presented in Appendix B.

As was shown in Keith et al. (2013); Coles et al. (2015), and later confirmed in Lentati et al. (2016), data for PSR J16037202 and PSR J1713+0747 show significant evidence for non-stationary extreme scattering events (ESEs), which are usually associated with the passage of high density plasma ‘blobs’ along the line of sight of a pulsar. ESEs are modeled as deterministic signals Lentati et al. (2016):

(21) |

by making use of shapelet basis function expansion:

(22) |

where is the epoch of ESE, stands for the characteristic length scale of ESE, is the -th Hermitian polynomial, is the number of terms used in the expansion, which is 3 in our case, is a vector of shapelet amplitudes. The inclusion of non-stationary ESEs in the noise model (see Table 3) leads to smaller DM spectral amplitudes and slightly steeper slopes, characterised by , which is consistent with results presented in Lentati et al. (2016).

#### iii.3.2 Frequentist methods

In the Frequentist framework, we use the method that was originally introduced in You et al. (2007) and further improved in Keith et al. (2013) for correcting DM variations. The basic idea works as follows. Timing residuals are separated into two components, one dependent on the radio wavelength, namely, dispersion measure variations – DM(t), and the other independent of the radio wavelength. The latter could contain red noise, GWs or dark matter signals. Since pulsar timing data are irregularly sampled, we use a linear interpolation scheme to estimate DM(t) at regular intervals. For the PPTA data, we estimate one DM(t) every 60-180 days using observations taken at three bands (10, 20, 50 cm). The time epochs and the estimated DM offsets are stored as DMOFF parameters in the TEMPO2 .par files. We model the red spin noise on data that have been corrected for DM variations, in which case, the noise covariance matrix contains only the white noise and spin noise terms.

Following the TEMPO2 convention, for our Frequentist analysis the intrinsic spin noise is parameterized using the following power-law spectrum

(23) |

where is an overall amplitude (normally expressed in yr), is the so-called corner frequency, is the power-law exponent. The covariance matrix for such a red noise process is given by

where with and being the ToA at the -th and -th observation respectively, is the modified Bessel function of second kind and is the Gamma function.

We follow the method described in Coles et al. (2011) to estimate red noise properties iteratively. We fit a power-law model of the form given by Eq. (23) to the power spectrum of timing residuals, leading to an initial estimate of the noise covariance matrix. We then use the Cholesky decomposition of this matrix to transform the data. The power spectrum of the transformed residuals should be white. We repeat the above procedure to obtain improved estimates of the spectrum. The iteration is considered converged if the whitened data show a sufficiently flat spectrum for which the spectral leakage is not dominant. The results are usually validated with simulations. We list our best estimates of red noise parameters in Table 3.

Pulsar Name | Bayesian | Frequentist | |||||||

note | note | yr | yr | ||||||

J04374715 | C | C | 3.5 | 0.08 | |||||

J06130200 | SC | C | 2.5 | 0.08 | |||||

J07116830 | C | SC | 4.0 | 0.08 | |||||

J10177156 | C | C | 6.0 | 1.0 | |||||

J10221001 | SC | C | 2.0 | 0.08 | |||||

J10240719 | SC | C | 3.0 | 0.08 | |||||

J10454509 | C | C | 3.0 | 0.3 | |||||

J11256014 | C | C | 3.0 | 0.2 | |||||

J14464701 | |||||||||

J15454550 | C | 3.0 | 0.1 | ||||||

J16003053 | SC | C | 2.0 | 0.08 | |||||

J16037202 | C | C | 3.0 | 0.08 | |||||

J16431224 | C | C | 1.5 | 0.08 | |||||

J17130747 | C | C | |||||||

J17302304 | C | C | 2.0 | 0.08 | |||||

J17325049 | SC | SC | |||||||

J17441134 | SC | SC | 6.0 | 1.0 | |||||

J18242452A | SC | C | 4.0 | 0.1 | |||||

J18320836 | |||||||||

J18570943 | SC | C | |||||||

J19093744 | C | C | 2.5 | 0.07 | |||||

J19392134 | C | C | 4.0 | 0.08 | |||||

J21243358 | SC | 5.0 | 1.0 | ||||||

J21295721 | SC | C | 2 | 0.08 | |||||

J21450750 | C | C | 1.0 | 0.08 | |||||

J22415236 | C | SC | |||||||

Including extreme scattering events | |||||||||

J16037202 | C | C | |||||||

J17130747 | C | C |

## Iv Search techniques and Results

### iv.1 Bayesian analysis

Within a Bayesian framework, the signal detection problem is addressed through model selection. Given the observational data, we wish to choose between two mutually exclusive hypotheses: the null hypothesis that the signal is absent and the alternative hypothesis that the signal is present. We compute the evidences , defined in Eq. (16), of the two hypotheses, and . Assuming a priori equal probability for both hypotheses, the following evidence ratio (commonly called Bayes factor) quantifies the support of against

(25) |

where are the parameters of the deterministic systematics, SSE errors in our case, which should be distinguished from dark matter signal parameters . In order to obtain accurate evidence estimates, we carry out numerical integration via MULTINEST with enabled importance nested sampling in constant efficiency mode. With the current PPTA data, we find a log Bayes factor of 2.1 in the frequency range [] Hz, implying that our data are consistent with containing only noise. When we extend the search frequency to Hz, the signal hypothesis is favored against the null hypothesis with . We suspect this is caused by the unmodeled perturbations of the mass and orbital elements of Mercury, for which the synodic period is days, corresponding to a frequency of Hz. We defer the investigation of this feature to a future work.

In order to set an upper limit on the signal amplitude within Bayesian framework, we perform the parameter estimation routine. By sampling from the posterior probabilities of model parameters, we can numerically marginalize over nuisance parameters, and get the posterior distribution for the amplitude . We define the Bayesian upper limit , such that of the samples from the posterior probability lie within the range :

(26) | |||||

We split the frequency range between and Hz into a number of small bins and find for each bin separately.

To reduce the computational costs of numerical marginalization, a common practice is to fix the noise model parameters to their maximum likelihood values Babak et al. (2016); Arzoumanian et al. (2014), determined from single pulsar analysis. However such a procedure can possibly lead to upper limits biased by a factor of Arzoumanian et al. (2014). In this work we allow both signal and correlated noise parameters to vary simultaneously. The white noise EFACs and EQUADs, which should have little or no correlation with dark matter parameters, are fixed to the maximum-likelihood values obtained from single pulsar analysis.

Recently, it was shown that the search for a stochastic GW background can be seriously affected by the uncertainties in the SSE Arzoumanian et al. (2018); Tiburzi et al. (2016). We employ a physical model BayesEphem to account for the SSE uncertainties that are most relevant for pulsar timing. The BayesEphem model has 11 parameters, including 4 parameters which describe the perturbations in the masses of outer planets, 1 parameter which is associated with the uncertainty in the semi-major axis of Earth-Moon barycenter orbit, and 6 parameters that characterize the perturbation of the Earth’s orbit due to errors in the Jovian average orbital elements. The BayesEphem modeling technique is described in Arzoumanian et al. (2018) in detail, and implemented in publicly available software packages, such as enterprise and NX01. The latter was used to put robust constrains on the amplitude of the FDM in this work.

The number of free parameters for the PPTA dataset is (see Table 2), where is the number of pulsars in PTA. In order to further reduce the computational costs, we have formed the “restricted dataset” by choosing the five best pulsars. As shown in Fig. 1, they contribute to more than 95% sensitivity of the full PPTA. Here pulsars are ranked according to their contribution to the squared signal-to-noise ratio ; see Eq. (29) in the next section. We carry out the calculations by adding detectable signals to 1000 noise realizations, sampled from individual pulsar noise posterior distribution obtained in Sec. III.3.1.

#### iv.1.1 Validation of the search results

In order to validate our upper limits and test the robustness of our algorithms, we have injected a signal with Hz and amplitude into our restricted data set. At this frequency, the amplitude of the injected signal is comparable to the Bayesian upper limit. In order to recover this signal we run the full Bayesian analysis, simultaneously accounting for both dark matter signal and noise. The posterior probabilities are demonstrated in Fig. 2, indicating successful recovery of the injected signal.

### iv.2 Frequentist analysis

In a Frequentist framework, signal detection is essentially a statistical hypothesis testing problem; we wish to choose between the null hypothesis and the signal hypothesis based on the observations. The task is to find an optimal statistic that maximizes the signal detection probability at a fixed false alarm probability. Following the Neyman-Pearson criterion, the log-likelihood ratio is an optimal statistic

(27) |

where we have used Eq. (11-12) to derive the second equality above, and the inner product between two time series x and y is defined as

(28) |

It is useful to define the signal-to-noise ratio in the following form

(29) |

where stands for the expectation value over a large number of noise realizations. In this work, we adopt as our detection statistic. For our Frequentist analysis, noise model parameters are fixed at their maximum likelihood values. The signal parameters in question are: the amplitude of dark matter induced gravitational-potential oscillations , oscillation frequency , phase parameters and ; see Eq. (9). It turns out that the statistic can be analytically maximized over and thus the parameter space that needs to be numerically searched over is dimensional. For our data this corresponds to 28 dimensions, making a grid-based search unfeasible. We employ the Particle Swarm Optimization technique Eberhart and Kennedy (1995), which has been demonstrated to be very effective for searches for continuous GWs with PTAs Wang et al. (2015b); Zhu et al. (2016). The detection statistic follows a distribution with one degree of freedom for noise-only data.

Since we find no evidence for statistically significant signals in the data, which is consistent with results from the Bayesian analysis as described in the previous subsection, we set upper limits on the dimensionless amplitude . We compute the 95% confidence upper limits for a number of frequency bins between and Hz. At each frequency, we compute the for simulated signals with random phase parameters and a fixed . The 95% confidence upper limit on corresponds to the amplitude at which 95% of signals result in . Here the threshold is chosen such that the expectation value for our detection statistic in the presence of signals, given by , is greater than the detection threshold that corresponds to 1% false alarm probability. It implies that: if there was a signal with an amplitude higher than present in the data, it would have been detectable with more than 95% probability.

### iv.3 Upper limits

Fig. 3 shows the 95% upper limits on the amplitude , calculated within Bayesian (black solid line) and frequentist frameworks (purple solid line). As one can see, Bayesian upper limits are a factor of 2-3 worse than frequentist upper limits in the low-frequency regime, while in the mid-to-high frequency range both methods produce comparable results. The difference might be predominantly attributed to the covariance between signal and noise (especially the red spin noise). Frequentist upper limits were calculated by fixing noise parameters at their maximum likelihood values, whereas we search simultaneously over signal and noise parameters in the Bayesian analysis.

The Bayesian upper limits, obtained with 5-year NANOGrav dataset Demorest et al. (2013), are also plotted as the thin dash-dotted (taken from Porayko and Postnov (2014)) and dashed (recalculated in this paper) lines. We note that upper limits presented in Ref. Porayko and Postnov (2014) were underestimated by a factor of ten due to the less conservative^{8}^{8}8We note that uniform priors result in upper limits that are a factor of higher than log-uniform priors. choice of prior (log-uniform) probability of the amplitude , as well as the non-inclusion of DM variations and additional white noise terms (EFAC and EQUAD).
From Fig. 3, one can see that our data set is a factor of five more sensitive to the dark matter signal than NANOGrav 5-year data at low frequencies, corresponding to boson masses eV. In the intermediate regime, the improvement is about a factor of two. This is expected because of our much longer data span and higher observing cadence. It is interesting to note that the upper limit curves in Fig. 3 exhibit similar frequency dependencies to the sky-averaged upper limits for continuous GWs (see, e.g., Zhu et al. (2014)). In Appendix A we present Frequentist upper limits obtained by including in the analysis only Earth terms. We also show how Bayesian upper limits are modified if different fixed SSE models are used.

## V Future prospects

In this section we discuss the future improvement in sensitivity of PTAs to the dark matter signal. In particular, the Five-hundred-meter Aperture Spherical Telescope (FAST Nan et al. (2011)) in China, MeerKAT Bailes et al. (2018) – a precursor for the planned Square Kilometre Array (SKA Lazio (2013)) – and ultimately the SKA, are expected to significantly increase the sensitivities of PTAs. With broad frequency bands and massive collecting areas, the radiometer noise for some of the brightest pulsars can be reduced from current 100 ns level down to below 10 ns Hobbs et al. (2014). However, it might be too optimistic to assume a white noise level of 10 ns because of the so-called jitter noise, which is thought to be associated with the intrinsic and stochastic variability in the shape of individual pulses Osłowski et al. (2011). Such a limitation implies that the timing precision stops improving for the brightest pulsars even when better instruments are used. The level of jitter noise can be approximately estimated with the following relation Shannon and Cordes (2012)

(30) |

where is the time of integration, and are the pulse width and pulse period, respectively. Note that the only way to reduce jitter noise is to increase . In comparison, the radiometer noise is given by Hobbs et al. (2014)

(31) |

where is the pulse profile signal-to-noise ratio, is the system-equivalent flux density (SEFD),
is the pulsar mean flux density and is the observing bandwidth. We adopt nominal SKA parameters^{9}^{9}9SKA1 system baseline V2 description https://www.skatelescope.org/, Jy,
MHz and set a fiducial minutes.

Pulsar Name | (ns) | (ns) | (ns) |
---|---|---|---|

J04374715 | 0.06 | 50.4 | 50.4 |

J10177156 | 4.6 | 13.7 | 14.5 |

J14464701 | 26.0 | 22.1 | 34.1 |

J15454550 | 15.6 | 36.1 | 39.3 |

J16003053 | 2.9 | 26.6 | 26.8 |

J1713+0747 | 0.8 | 35.1 | 35.1 |

J17441134 | 3.9 | 41.2 | 41.4 |

J18320836 | 3.7 | 14.2 | 14.8 |

J19093744 | 1.2 | 11.2 | 11.3 |

J22415236 | 1.5 | 15.4 | 15.5 |

Table 4 lists white noise budgets (, and the total white noise ) expected in the FAST/SKA era for ten PPTA pulsars that have the lowest value of . As one can see, for the SKA, jitter noise will dominate over the radiometer noise for the majority of bright pulsars. In order to realistically estimate the PTA sensitivity in the FAST/SKA era, we use the total white noise given in Table 4 plus the intrinsic spin noise (where appropriate) with parameters determined from the Bayesian analysis.

Fig. 4 shows forecasted upper limits on the density of FDM in the Galaxy for three cases, all assuming a data span of ten years.
Case a) is a conservative PTA that includes only ten pulsars as listed in Table 4 and an observing cadence of once every 14 days. Upper limits in this case are obtained by running full Bayesian analysis of simulated data.
We analytically scale this limit curve to two more ambitious cases^{10}^{10}10Note that the scaling factor should be a good approximation at high frequencies where red noise plays a less important role..
We increase the number of pulsars to 100 in case b), leading to a factor of improvement.
For case c), we further increase the cadence to once every day and adopt an integration time of two hours, providing another factor of improvement. Case c) might be an interesting option in the SKA era since small radio telescopes (compared to SKA/FAST), such as Parkes, can be dedicated for high-cadence and long integration observations of the brighter pulsars.

As one can see from Fig. 4, we will be able to constrain the contribution of FDM to the local dark matter density below 10% for eV in ten years under the conservative assumption for SKA sensitivity. However, it is more challenging for boson masses above eV; we estimate that decade-long observations of hundreds of pulsars timed at nearly daily cadence with precision ns are necessary to place interesting limits.

There are a couple of ways to improve our analysis. First, the coherence between pulsar terms and Earth terms can be used to enhance the sensitivity. When a pulsar and the Earth are located within a de Broglie wavelength , the oscillation phases, which have been assumed to be independent in the current analysis, are correlated. However, for eV, this effect will have no impact on the current results, since and no pulsars have been found within 60 pc to the Earth. Another interesting point is that pulsars that are close to each other within also experience phase-coherent oscillations De Martino et al. (2017). We plan to explore these features in a future work.

Second, the oscillation amplitude is proportional to the local dark matter density. Thus, in contrast to the amplitude of the Earth term, the amplitude of the pulsar term varies from pulsar to pulsar; see Eq. (7). In -FDM cosmological simulations De Martino et al. (2017); Schive et al. (2014), it was shown that due to wave interference the dark matter forms gritty pattern with typical granule size of around . When averaged over scales, the periphery ( kpc) density profile is similar to the classical Navarro-Frenk-White (NFW) profile, whereas a distinct density peak is seen in the central regions (usually called solitonic core, see Schive et al. (2014) for details).

Fig. 5 shows the expected signal amplitude for PPTA pulsars assuming the NFW dark matter density profile Navarro et al. (1996) with parameters from Sofue (2012). As one can see, pulsars closer to the Galactic Center provide better sensitivity to the dark matter signal. The amplitude of the dark matter signal becomes even larger than NFW prediction within the central solitonic core (kpc) De Martino et al. (2017).
For the current PPTA sample, PSR J18242452A is expected to have the largest signal amplitude, a factor of larger than other pulsars^{11}^{11}11The density of the scalar field dark matter in globular clusters is not expected to deviate significantly from the general trend as is larger than typical sizes of globular clusters. Thus, the amplitude of the oscillation at J18242452A, located in a globular cluster, is expected to follow the NFW prediction.. However, this pulsar is nearly the worst timer in PPTA (see Table 1 and Fig. 1). Existing and future pulsar surveys might help find high quality millisecond pulsars close to the Galactic Center and thus provide better sensitivity to the dark matter searches Kramer et al. (2004).

## Vi Conclusions

Pulsar timing is a powerful tool to study a wide variety of astrophysical phenomena. By exploiting precision timing observations from an array of the most stable millisecond pulsars, PTAs allow us to measure minute correlations in the ToAs of different pulsars. Like continuous GWs from individual supermassive binary black holes, FDM in the Galaxy produces periodic variations in pulsar ToAs. We perform a search for evidence of ultralight dark matter in the latest PPTA data set. Finding no statistically significant signals, we place upper limits on the dark matter density: for boson mass eV, our analysis constrains the density below with 95% confidence; at eV, our upper limits remain three orders of magnitude above the local dark matter density inferred from kinematics measurements of stars in the Galaxy Sivertsson et al. (2018a).

We derived the noise properties of PPTA data and obtain dark matter constraints using both Bayesian and Frequentist methods. Our upper limits from the two methods are broadly consistent. We reanalyzed the NANOGrav 5-yr data set and found that the PPTA data result in a factor of two to five improvement in dark matter constraints. We studied potential systematics due to SSE errors in our analysis and found that the search for ultralight dark matter is insensitive to such errors. We have ignored effects from instabilities in terrestrial time standards; such clock errors produce a monopolar broad-band noise Hobbs et al. (2012). Whereas this effect should be distinguishable from the sinusoidal ToA variations due to ultralight dark matter, one needs to include it in a future study to quantitatively assess the impact.

We forecasted the PTA sensitivity in the FAST/SKA era while accounting for realistic noise levels. We found that observing the ten best PPTA pulsars for ten years would constrain the density of FDM below for eV, about 10% of measured total dark matter density. At eV, our projected limit is around ; for higher boson masses, the upper limits increase as . Above eV, the projected limits are more than one order of magnitude above the local dark matter density. To place interesting limits in this mass range, an ambitious timing program in which hundreds of pulsars timed with daily cadence and high precision ( ns) precision for more than a decade is required. Finally, we point out that high-quality pulsars in the vicinity of the Galactic Center will be ideal tools to test the fuzzy dark matter hypothesis.

## Acknowledgements

The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. N.K.P. acknowledges the support from IMPRS Bonn/Cologne and the Bonn-Cologne Graduate School (BCGS). X.Z., M.B., D.J.R., R.M.S. & L.W. are supported by ARC CE170100004. X.Z. & L.W. are additionally supported by ARC DP150102988. L.H. is supported in part by NASA NXX16AB27G and DOE DE-SC0011941. Work at NRL is supported by NASA. P.D.L. is supported through ARC FT160100112 and ARC DP180103155. M.B., S.O. & R.S. acknowledge support through the ARC Laureate Fellowship grant FL150100148. J.W. is supported by Qing Cu Hui of Chinese Academy of Sciences (CAS). We acknowledge the Institute for Theoretical and Experimental physics and, in particular, Sergey Blinnikov for providing computing facilities. The authors would like to thank Dr. Maxim Pshirkov, Michael Ivanov, Dr. Nicolas Caballero and Dr. David Champion for fruitful discussions. We would like to thank the anonymous referees for useful comments.

## References

- Zwicky (1933) F. Zwicky, Helvetica Physica Acta 6, 110 (1933).
- Zwicky (1937) F. Zwicky, Astrophys. J. 86, 217 (1937).
- Smith (1936) S. Smith, Astrophys. J. 83, 23 (1936).
- Koopmans and Treu (2003) L. V. E. Koopmans and T. Treu, Astrophys. J. 583, 606 (2003), astro-ph/0205281 .
- Clowe et al. (2004) D. Clowe, A. Gonzalez, and M. Markevitch, Astrophys. J. 604, 596 (2004), astro-ph/0312273 .
- Tegmark et al. (2004) M. Tegmark et al., Astrophys. J. 606, 702 (2004), astro-ph/0310725 .
- Planck Collaboration et al. (2016) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al., A&A 594, A13 (2016), arXiv:1502.01589 .
- Bertone et al. (2005) G. Bertone, D. Hooper, and J. Silk, Phys. Rep. 405, 279 (2005), hep-ph/0404175 .
- Primack (2012) J. R. Primack, Annalen der Physik 524, 535 (2012).
- Battaglieri et al. (2017) M. Battaglieri et al., (2017), arXiv:1707.04591 [hep-ph] .
- Arvanitaki et al. (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
- Hui et al. (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 .
- Preskill et al. (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B120, 127 (1983), [,URL(1982)].
- Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. B120, 133 (1983), [,URL(1982)].
- Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B120, 137 (1983).
- Svrcek and Witten (2006) P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hep-th/0605206 [hep-th] .
- Hu et al. (2000) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), astro-ph/0003365 .
- Bullock and Boylan-Kolchin (2017) J. S. Bullock and M. Boylan-Kolchin, Ann. Rev. Astron. Astrophys. 55, 343 (2017), arXiv:1707.04256 [astro-ph.CO] .
- Hložek et al. (2018) R. Hložek, D. J. E. Marsh, and D. Grin, Mon. Not. R. Astron. Soc. (2018), 10.1093/mnras/sty271, arXiv:1708.05681 .
- Iršič et al. (2017) V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, Phys. Rev. Lett. 119, 031302 (2017), arXiv:1703.04683 .
- Kobayashi et al. (2017) T. Kobayashi, R. Murgia, A. De Simone, V. Iršič, and M. Viel, Phys. Rev. D 96, 123514 (2017), arXiv:1708.00015 .
- Bowman et al. (2018) J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen, and N. Mahesh, Nature 555, 67 (2018).
- Schneider (2018) A. Schneider, (2018), arXiv:1805.00021 [astro-ph.CO] .
- Lidz and Hui (2018) A. Lidz and L. Hui, (2018), arXiv:1805.01253 [astro-ph.CO] .
- Sullivan et al. (2018) J. M. Sullivan, S. Hirano, and V. Bromm, ArXiv e-prints (2018), arXiv:1809.01679 .
- Calabrese and Spergel (2016) E. Calabrese and D. N. Spergel, Mon. Not. Roy. Astron. Soc. 460, 4397 (2016), arXiv:1603.07321 [astro-ph.CO] .
- Deng et al. (2018) H. Deng, M. P. Hertzberg, M. H. Namjoo, and A. Masoumi, (2018), arXiv:1804.05921 [astro-ph.CO] .
- Bar et al. (2018) N. Bar, D. Blas, K. Blum, and S. Sibiryakov, ArXiv e-prints (2018), arXiv:1805.00122 .
- Schive et al. (2014) H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics 10, 496 (2014), arXiv:1406.6586 .
- Mocz and Succi (2015) P. Mocz and S. Succi, Phys. Rev. E91, 053304 (2015), arXiv:1503.03869 [physics.comp-ph] .
- Veltmaat and Niemeyer (2016) J. Veltmaat and J. C. Niemeyer, Phys. Rev. D94, 123523 (2016), arXiv:1608.00802 [astro-ph.CO] .
- Nori and Baldi (2018) M. Nori and M. Baldi, (2018), arXiv:1801.08144 [astro-ph.CO] .
- Veltmaat et al. (2018) J. Veltmaat, J. C. Niemeyer, and B. Schwabe, (2018), arXiv:1804.09647 [astro-ph.CO] .
- Khmelnitsky and Rubakov (2014) A. Khmelnitsky and V. Rubakov, JCAP 2, 019 (2014), arXiv:1309.5888 .
- Porayko and Postnov (2014) N. K. Porayko and K. A. Postnov, Phys. Rev. D 90, 062008 (2014), arXiv:1408.4670 .
- De Martino et al. (2017) I. De Martino, T. Broadhurst, S.-H. H. Tye, T. Chiueh, H.-Y. Schive, and R. Lazkoz, Physical Review Letters 119, 221103 (2017), arXiv:1705.04367 .
- Blas et al. (2017) D. Blas, D. L. Nacir, and S. Sibiryakov, Physical Review Letters 118, 261102 (2017), arXiv:1612.06789 [hep-ph] .
- Sazhin (1978) M. V. Sazhin, Soviet Astronomy 22, 36 (1978).
- Detweiler (1979) S. Detweiler, Astrophys. J. 234, 1100 (1979).
- Hellings and Downs (1983) R. W. Hellings and G. S. Downs, Astrophys. J. 265, L39 (1983).
- Foster and Backer (1990) R. S. Foster and D. C. Backer, Astrophys. J. 361, 300 (1990).
- Manchester et al. (2013) R. N. Manchester et al., Publ. Astron. Soc. Aust. 30, e017 (2013), arXiv:1210.6130 [astro-ph.IM] .
- McLaughlin (2013) M. A. McLaughlin, Classical and Quantum Gravity 30, 224008 (2013), arXiv:1310.0758 [astro-ph.IM] .
- Kramer and Champion (2013) M. Kramer and D. J. Champion, Classical and Quantum Gravity 30, 224009 (2013).
- Hobbs et al. (2010a) G. Hobbs et al., Classical and Quantum Gravity 27, 084013 (2010a), arXiv:0911.5206 [astro-ph.SR] .
- Verbiest et al. (2016) J. P. W. Verbiest et al., Mon. Not. R. Astron. Soc. 458, 1267 (2016), arXiv:1602.03640 [astro-ph.IM] .
- Shannon et al. (2013) R. M. Shannon et al., Science 342, 334 (2013), arXiv:1310.4569 .
- Zhu et al. (2014) X.-J. Zhu et al., Mon. Not. R. Astron. Soc. 444, 3709 (2014), arXiv:1408.5129 .
- Wang et al. (2015a) J. B. Wang et al., Mon. Not. R. Astron. Soc. 446, 1657 (2015a), arXiv:1410.3323 .
- Shannon et al. (2015) R. M. Shannon et al., Science 349, 1522 (2015), arXiv:1509.07320 .
- Demorest et al. (2013) P. B. Demorest et al., Astrophys. J. 762, 94 (2013), arXiv:1201.6641 .
- Ellis and van Haasteren (2017) J. Ellis and R. van Haasteren, “jellis18/pal2: Pal2,” (2017).
- Taylor and Baker (2017) S. Taylor and P. T. Baker, “stevertaylor/nx01 1.2,” (2017).
- Brdar et al. (2018) V. Brdar, J. Kopp, J. Liu, P. Prass, and X.-P. Wang, Phys. Rev. D 97, 043001 (2018), arXiv:1705.09455 [hep-ph] .
- Zhu et al. (2016) X.-J. Zhu, L. Wen, J. Xiong, Y. Xu, Y. Wang, S. D. Mohanty, G. Hobbs, and R. N. Manchester, Mon. Not. R. Astron. Soc. 461, 1317 (2016), arXiv:1606.04539 [astro-ph.IM] .
- Bovy and Tremaine (2012) J. Bovy and S. Tremaine, Astrophys. J. 756, 89 (2012), arXiv:1205.4033 [astro-ph.GA] .
- Read (2014) J. I. Read, Journal of Physics G Nuclear Physics 41, 063101 (2014), arXiv:1404.1938 .
- Sivertsson et al. (2018a) S. Sivertsson, H. Silverwood, J. I. Read, G. Bertone, and P. Steger, Mon. Not. R. Astron. Soc. 478, 1677 (2018a), arXiv:1708.07836 .
- Sivertsson et al. (2018b) S. Sivertsson, H. Silverwood, J. I. Read, G. Bertone, and P. Steger, Mon. Not. R. Astron. Soc. 478, 1677 (2018b), arXiv:1708.07836 .
- Lentati et al. (2016) L. Lentati et al., Mon. Not. R. Astron. Soc. 458, 2161 (2016), arXiv:1602.05570 [astro-ph.IM] .
- Manchester et al. (2005) R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, AJ 129, 1993 (2005), arXiv:astro-ph/0412641 .
- Folkner et al. (2007) W. Folkner, E. Standish, J. Williams, and D. Boggs, Jet Propulsion Laboratory, Memorandum IOM 343R-08-003 (2007).
- Hobbs et al. (2006) G. B. Hobbs, R. T. Edwards, and R. N. Manchester, Mon. Not. R. Astron. Soc. 369, 655 (2006), astro-ph/0603381 .
- Edwards et al. (2006) R. T. Edwards, G. B. Hobbs, and R. N. Manchester, Mon. Not. R. Astron. Soc. 372, 1549 (2006), astro-ph/0607664 .
- van Haasteren et al. (2009) R. van Haasteren, Y. Levin, P. McDonald, and T. Lu, Mon. Not. R. Astron. Soc. 395, 1005 (2009), arXiv:0809.0791 .
- van Haasteren and Levin (2013) R. van Haasteren and Y. Levin, Mon. Not. R. Astron. Soc. 428, 1147 (2013), arXiv:1202.5932 [astro-ph.IM] .
- Press et al. (1996) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes in C, Vol. 2 (Cambridge university press Cambridge, 1996).
- Boynton et al. (1972) P. E. Boynton, E. J. Groth, D. P. Hutchinson, G. P. Nanos, Jr., R. B. Partridge, and D. T. Wilkinson, Astrophys. J. 175, 217 (1972).
- Blandford et al. (1984) R. Blandford, R. Narayan, and R. W. Romani, Journal of Astrophysics and Astronomy 5, 369 (1984).
- Hobbs et al. (2010b) G. Hobbs, A. G. Lyne, and M. Kramer, Mon. Not. R. Astron. Soc. 402, 1027 (2010b), arXiv:0912.4537 .
- William (1989) H. W. William, SIAM Revew 31, 221 (1989).
- van Haasteren and Vallisneri (2015) R. van Haasteren and M. Vallisneri, Mon. Not. R. Astron. Soc. 446, 1170 (2015), arXiv:1407.6710 [astro-ph.IM] .
- Arzoumanian et al. (2014) Z. Arzoumanian et al., Astrophys. J. 794, 141 (2014), arXiv:1404.1267 .
- Lentati et al. (2013) L. Lentati, P. Alexander, M. P. Hobson, S. Taylor, J. Gair, S. T. Balan, and R. van Haasteren, Phys. Rev. D 87, 104021 (2013), arXiv:1210.3578 [astro-ph.IM] .
- Lentati et al. (2014) L. Lentati, P. Alexander, M. P. Hobson, F. Feroz, R. van Haasteren, K. J. Lee, and R. M. Shannon, Mon. Not. R. Astron. Soc. 437, 3004 (2014), arXiv:1310.2120 [astro-ph.IM] .
- Feroz et al. (2009) F. Feroz, M. P. Hobson, and M. Bridges, Mon. Not. R. Astron. Soc. 398, 1601 (2009), arXiv:0809.3437 .
- Feroz et al. (2013) F. Feroz, M. P. Hobson, E. Cameron, and A. N. Pettitt, ArXiv e-prints (2013), arXiv:1306.2144 [astro-ph.IM] .
- Keith