The second data release from the European Pulsar Timing Array: VI. Challenging the ultralight dark matter paradigm

Pulsar Timing Array experiments probe the presence of possible scalar or pseudoscalar ultralight dark matter particles through decade-long timing of an ensemble of galactic millisecond radio pulsars. With the second data release of the European Pulsar Timing Array, we focus on the most robust scenario, in which dark matter interacts only gravitationally with ordinary baryonic matter. Our results show that ultralight particles with masses $10^{-24.0}~\text{eV} \lesssim m \lesssim 10^{-23.3}~\text{eV}$ cannot constitute $100\%$ of the measured local dark matter density, but can have at most local density $\rho\lesssim 0.3$ GeV/cm$^3$.

Introduction.-Thenature of Dark Matter (DM) is arguably one of the most fascinating and mysterious questions that we are struggling to answer.Galaxy rotation curves [1,2], the peculiar motion of clusters [3,4], the Bullet Cluster system [5], measurements of cosmological abundances from Cosmic Microwave Background (CMB) and Baryonic Acoustic Oscillation (BAO) observations [6,7] all point to the existence of a hithertounseen type of matter, constituting roughly 26% of the current energy density of the Universe and interacting mostly gravitationally with the Standard Model of particle physics.The standard Cold Dark Matter (CDM) paradigm describes successfully many aspects of the large-scale structure of the Universe, but struggles to predict what we observe at scales smaller than the ∼kpc.For instance, observations favor a constant density profile in the inner part of galaxies, while CDM predicts a steep power-law-like behavior (cusp-core problem) [8][9][10].Furthermore, well-known issues are associated with the discrepancy between the observed and expected number of Milky Way (MW) satellites (missing satellite problem) [11,12] and with ΛCDM simulations showing that the most massive subhaloes of the MW would be too dense to host any of its bright satellites (too-big-to-fail problem) [13].Moreover, recent anomalies in gravitationally lensed images [14] seem to disfavor the long-standing Weakly-Interacting-Massive-Particles (WIMPs) hypothesis for CDM.While some of these issues might be alleviated by invoking baryonic feedback mechanisms [15], e.g.Active Galactic Nuclei (AGN) [16] and/or supernova feedback [17][18][19][20][21][22], it is still unclear how the flat density profile of dwarf galaxies (e.g.Fornax [23]), with almost no baryonic activity in the center, can be explained without invoking a novel mechanism.An intriguing alternative is to consider the possibility that DM is fuzzy, i.e. an ultra-light scalar field (m ϕ ∼ 10 −22 eV) or axion-like particle, whose wavelike nature suppresses structure formation on scales smaller than the de Broglie wavelength, while maintaining all the achievements of the CDM paradigm on large scales.Moreover, the existence of ultralight scalars can also be motivated on a * csmarra@sissa.itmore theoretical ground: in particular, axion-like particles generically arise in string theory compactifications as Kaluza-Klein zero modes of antisymmetric tensor fields [24][25][26].
A wealth of studies have been carried out to probe the existence of ultra-light dark matter (ULDM), ranging from CMB observables to Lyman-α and stellar kinematics.Specifically, the integrated Sachs-Wolfe effect on CMB anisotropies rules out masses m ϕ ≲ 10 −24 eV [27], while Lyman-α gives a lower bound m ϕ ≳ 10 −21 eV for ultra-light candidates constituting more than ∼ 30% of DM1 [28][29][30][31][32]. Stellar orbit kinematics in ultra-faint dwarf (UFD) galaxies might even be able to bound the scalar field mass to be m ϕ ≳ 10 −19 eV, although this is still under debate [33,34].However, the sensitivity of non-CMB constraints to uncertainties in the modeling of small scale structure properties [35,36] makes it compelling to rely on complementary and independent probes.It was shown by Khmelnitsky and Rubakov [37] that the presence of ULDM induces an oscillating gravitational potential that affects the light travel time of radio pulses emitted by pulsars.Therefore, Pulsar Timing Arrays (PTAs) stand out as promising experiments to test the presence of ULDM particles in the MW.Previous PTA searches placed 95% upper limits on the local energy density of ULDM at 3 × 10 −24 eV to ≲ 1 GeV/cm 3 [38][39][40].
In this work, which is complementary to the European Pulsar Timing Array (EPTA) interpretation effort [41], we focus on a specific range of ULDM masses and constrain the local ULDM density to values below the observed local DM density.We do so by analyzing the effect of ULDM on the times of arrival (TOAs) of pulsar radio beams.Therefore, if ULDM particles exist in the mass range that we consider, they cannot constitute all of the observed DM.
Models.-As we only have gravitational evidence of DM, we focus on an ultra-light scalar field with negligible self-interactions and no couplings with the Standard Model.The action for this field can be written as Due to its high occupation number and non-relativistic nature, the ULDM scalar field can be thought as a classical wave [37]: where m ϕ is the mass of the scalar field, γ(⃗ x) is a spacedependent phase and φ(⃗ x) accounts for the pattern of interference in the proximity of ⃗ x caused by the wave-like nature of ULDM.The scalar field density ρ ϕ is conveniently normalized to the local DM density ρ DM , which can be determined e.g. by fitting the MW rotation curve or, in a more refined way, by studying the vertical oscillations of disc stars [42][43][44][45].In the following, we assume a fiducial value ρ DM ≈ 0.4 GeV/cm 3 .The oscillating nature of ULDM induces an oscillating gravitational potential leading to a periodic displacement δt DM in the TOAs of radio pulses emitted by pulsars, which can be written as [37,39]: )) are related to the phases of Eq. ( 2) evaluated at the pulsar (Earth) location, with d p standing for the pulsar-Earth distance.The amplitude in Eq. ( 4) is computed assuming a constant DM density background and possible deviations caused by the wave-like nature of the ultralight scalar field are parametrized in terms of the pulsar (Earth) dependent phase factors φ2 ( ⃗ x p ) ≡ φ2 P ( φ2 ( ⃗ x e ) ≡ φ2 E ).The approximation of constant DM density across pulsars is sufficient, as their distances from Earth are all ∼kpc and subject to large measurement uncertainties [39].Notice that accurate measurements of pulsar-Earth distances might help to reduce the number of free parameters in the limit in which γ( ⃗ x p ) = γ( ⃗ x e ).Moreover, precise determination of pulsar positions could provide us with more information about the dark matter density in its surroundings.On scales smaller than the de Broglie wavelength, the ULDM scalar field oscillates coherently, with the same amplitude φ (see Eq. ( 2)).Since the typical ULDM velocity is expected to be v ϕ ∼ 10 −3 , the coherence length is approximately Therefore, φ2 E and φ2 P are: • uncorrelated if the coherence length of ULDM is less than the average inter-pulsar and pulsar-Earth separation.In this case, φ2 E and φ2 P will thus be separate parameters; • correlated if the coherence length of ULDM is larger than the inter-pulsar and pulsar-Earth separations and encloses the typical Galacto-centric region tested by the most precise MW rotation curves measurements (roughly the inner ∼ 20 kpc [46]).In this case, φ2 E = φ2 P for all the pulsars.Moreover, rotation curves also sample from the same coherence patch, and thus measure the local abundance ρ DM of DM.Therefore, the stochastic parameter φ2 can be safely absorbed in a redefinition of Ψ c .
• pulsar-correlated if the coherence length of ULDM is larger than the inter-pulsar and pulsar-Earth separations, but smaller than the typical Galactocentric radius sampled by rotation curves.In this case, φ2 E = φ2 P for all the pulsars.However, DM density estimates from rotation curves average over different patches.We therefore keep φ2 as a free parameter and consistently marginalize over it.In this way, the limits on ρ DM obtained from pulsars will constrain the same quantity measured by rotation curves.
We perform the analysis in the three limits above, noting that the fully correlated limit has not been considered in the previous studies [47,48].From Eq. ( 5), recalling that the pulsar-Earth distance is O(kpc) for the observed systems, it follows that the correlated regime is an excellent approximation for masses lower than m ϕ ∼ 2 × 10 −24 eV; the pulsar-correlated regime holds for 2 × 10 −24 eV ≲ m ϕ ≲ 5 × 10 −23 eV and the uncorrelated regime is valid for m ϕ ≳ 5 × 10 −23 eV.Dataset and methodology.-TheEPTA monitors 42 millisecond radio pulsars with five telescopes located in France, Germany, Italy, the Netherlands, and the United Kingdom.The second data release (DR2) of the EPTA contains 24.7 years of observations of pulse arrival times of 25 pulsars, surveyed with an approximate cadence of once every 3 weeks [49], which translates into a Nyquist frequency of approximately 3 × 10 −7 Hz.TOAs are measured at the position of the Solar System Barycenter, and are primarily described by pulsar-specific deterministic timing models accounting for the position of each pulsar in the sky, spin down, proper motion, the presence of a binary companion, etc.The timing models are provided with the data set.Deviations between the TOAs predicted by these models and the measured TOAs are referred to as the timing residuals, δt.The residuals contain contributions from various sources of noise, from variations in the dispersion measure to irregularities in pulsar rotation, but they may also contain signals of astrophysical interest.The sources of noise are identified as part of the noise analysis of the EPTA DR2 [50].The DR2 data set also hints at growing evidence for the stochastic gravitational-wave background, which manifests itself as a temporally-correlated stochastic process with a hallmark inter-pulsar correlation signature of General Relativity [51].Our work is complementary to the EPTA-wide interpretation effort of the spatial and temporal correlations in DR2 [41].
We use Bayesian inference techniques to search for the ULDM signal while simultaneously fitting timing model parameters and all known sources of noise to the data, to correctly marginalize over the associated uncertainties.The likelihood of the timing residuals, L(δt|θ), given the parameters of the models, θ, is [52][53][54][55][56] ln L(δt|θ This is a time-domain Gaussian likelihood, multivariate with respect to a number of observations, i.e., δt has dimension equal to the number of observations.The contribution of ULDM from Eq. ( 3) is added to µ, which also contains contributions from the timing model [49] and noise processes, according to the noise analysis of Ref. [50].The diagonal part of the covariance matrix C contains TOA measurement uncertainties that include temporally-uncorrelated "white" noise.Contributions from temporally correlated "red" noise may be added as off-diagonal elements in C.However, for computational efficiency, red noise contributions are modeled in µ [53,54].The priors π(θ) are set based on Table I, see Ref. [47] for further details.To obtain a sufficient amount of posterior samples across the mass-frequency parameter space, the search is effectively performed across equallyspaced segments of π(m ϕ ), which we refer to as bins.
The measurements of parameters are obtained as posterior distributions, P(θ|d) ∝ L(δt|θ)π(θ).The posteriors are evaluated using the Parallel-Tempering-Markov-Chain Monte-Carlo sampler [57] implemented in enterprise [55] and enterprise extensions [56] .The conclusion about the presence or absence of the ULDM signal in the data is based on the Bayesian odds ratio.In our case, that is equal to the Bayes factor, because we assume the prior odds of both scenarios to be equal.We evaluate Bayes factors, B, using the Savage-Dickey density ratio [58].In particular, finding ln B ≳ 5 would indicate robust evidence for the ULDM signal.
There is strong evidence for a temporally-correlated red signal in EPTA DR2, characterized by the same Fourier spectrum of δt in all pulsars [51].This signal may contain contributions from the stochastic gravitationalwave background [59].Because in the 24.7-yr data set this signal does not show significant evidence for interpulsar correlations (unlike in the 10.3-yr data set [41,[49][50][51]), we model it as a spatially uncorrelated red noise process.Individual Fourier components of this broadband signal may contaminate frequencies at which the presence of ULDM is evaluated.Thus, this common red noise signal is included in the null hypothesis, ∅, along with the pulsar-intrinsic noise.The signal hypothesis is based on adding ULDM to the null hypothesis.
Results.-We carry out the search for ULDM in the correlated, pulsar-correlated and uncorrelated limit for φ2 E and φ2 P with the parameters in Table I.The factors φ2 E and φ2 P are drawn from an exponential prior, to correctly model the stochastic nature of the ULDM field [60,61].We find no evidence for a signal in the mass range m ϕ ∼ [10 −24 eV, 10 −22 eV].The largest ln B we find across frequency-mass bins is < 1, i.e. the null hypothesis is favored.Therefore, we calculate the 95% upper limits on the signal amplitude Ψ c and, through Eq. ( 4), on the scalar field density ρ ϕ .The results are shown in Fig. 1.
A highlight of our study is that not only does EPTA DR2 yield more stringent constraints than previous results [38,39], but it also rules out that particles with masses m ϕ ∼ [10 −24 eV, 10 −23.3 eV] can be 100% of the observed local DM density.In particular, the scalar field density is ρ ϕ ≲ 0.15 GeV/cm 3 in the mass range m ϕ ∼ [10 −24 eV, 10 −23. 7eV], while it is constrained to ρ ϕ ≲ 0.30 GeV/cm 3 between m ϕ ∼ [10 −23.7 eV, 10 −23. 4 eV] .Furthermore, the correlated limit in Fig. 1 confirms Lyman-α bounds, which exclude ULDM in this mass range unless it constitutes less than 30% of DM [30].It is worth noticing that the low-frequency end of Fig. 1 extends below the näive expectation f = 1.3 nHz corresponding to the inverse of the observation time T obs = 24.7 yr.In fact, while an ULDM candidate in this mass region does not complete an oscillation cycle during the observation timescale, the signal can still be approximated by a polynomial expansion in (m ϕ t) [47].The sensitivity in this region is limited by the simultaneous fitting to pulsar spin frequency derivatives [62,63].PTAs are only sensitive to the (m ϕ t) 3 term, as the first terms in the expansion are degenerate with the timing model [64].However, since the expected amplitude Ψ c of an ULDM candidate increases as its mass decreases, we can still obtain competitive constraints at low frequency.While, in principle, our analysis could be pushed to even lower masses [65], we choose to focus on the region m ϕ ≳ 10 −24 eV to comply with the aforementioned CMB bounds.We find that the significant improvement in sensitivity to ULDM at low frequencies arises thanks to the larger data span of EPTA DR2, in accordance with the theoretical sensitivity scaling proposed in Eq. ( 13) of Ref. [65].In particular, because of the longer data span, we expect EPTA DR2 limits to be better than NANOGrav [66] ones by a factor of roughly ∼ 3.6, which is in agreement with what observed.At high frequencies, we find that the advantage of the long timing baseline compared to NANOGrav diminishes, also in accordance with the scaling, as pulsar white noise levels become more important.We also performed an identical analysis of the 10-year subset of EPTA DR2 [49,51], as well as of the MeerTime data [67], which yield less stringent upper limits in agreement with the scaling.For comparison, the bounds in both the correlated and uncorrelated limit for the 10-year subset of the EPTA, shown in Fig. 22 in Ref. [41], appear at the level of the pulsar-correlated limit in Fig. 1 of this manuscript.
In the following, we clarify some specific aspects of our results.First, we notice that a similar analysis has been done by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) collaboration [66].TABLE I: Parameters employed for the search along with their respective priors.In the correlated limit, the parameters φ2 E , φ2 P are accounted for by a redefinition of Ψ c , while in the pulsar-correlated regime φ2 E = φ2 P = φ2 is a free parameter.FIG.1: Upper limits on ULDM, and namely on the dimensionless amplitude (Ψ c , left panel ) and the ULDM fraction of the local DM density ρ DM = 0.4GeV/cm 3 (ρ ϕ /ρ DM , right panel ), at 95% credibility.The bottom horizontal axes show the ULDM particle mass, whereas the top horizontal axes show the equivalent oscillation frequency of the scalar field.The upper limits from previous searches [38,39] are shown for comparison.As a reference, we highlight the frequency T −1 obs .In the right panel, we zoom in on the excluded ULDM masses.The horizontal dotted line represents the value of ρ ϕ that would saturate the local DM density.Notice that based on our results ULDM particles with mass −24.0 < log 10 (m ϕ /eV) < −23.7 can only make up at most 30 − 40 % of the total DM energy density, while particles with mass −23.7 < log 10 (m ϕ /eV) < −23.3 can contribute at most up to ∼ 70 %.
There, the upper limits provided in the correlated and uncorrelated scenarios differ at low frequency.This can be understood by noticing that the correlated limit of NANOGrav corresponds to our pulsar-correlated limit.However, in the low mass limit of Fig. 1, the pulsars, the Earth and the stellar and gaseous tracers used for rotation curves estimates lie well within the area spanned by the coherence length; thus, one can only measure the combination Ψ 0 c = Ψ c φ2 , which represents the realization of DM in our Galaxy.Therefore, we remove the φ2 E = φ2 P ≡ φ2 parameter in the correlated limit, as it can be accounted for by a redefinition of Ψ c .Fitting for Ψ c and φ2 separately, instead, introduces an additional uncertainty, which leads our pulsar-correlated analysis to produce a similar mismatch as the one found in Ref. [66], as shown in Fig. 1.Second, Fig. 1 hints at a steep increase in the upper limits at m ϕ ≳ 10 −23.2 eV.In fact, we report the presence of excess signal power on top of the common red noise process, corresponding to a mass of m ϕ ≃ 10 −23 eV and an amplitude of Ψ c ≃ 6 × 10 −14 or, equivalently, a density of ρ ϕ = 90 GeV/cm 3 .At face value, this excess is not compatible with an ULDM candidate, as the corresponding density is outside the local DM measurement uncertainties [42][43][44][45].Moreover, such a mass would be in tension with astrophysical bounds, as extensively discussed in the introduction [27][28][29][30][31][32][33][34].Anyway, the Bayesian odds ratio suggests that it is still consistent with noise (ln B ∼ 0.1).We find a similar excess in the analysis of 10-yr subset of the EPTA DR2 [41].Moreover, the boson mass corresponding to the excess also matches the frequency of the continuous gravitational wave (CGW) candidate studied in [68].This motivates further investigations as part of the International Pulsar Timing Array [69].
Conclusions.-ULDM is a theoretically motivated paradigm that may alleviate the small-scale crisis of structure formation.Here, we focused on the most robust scenario, in which ULDM features only gravitational interactions.These interactions produce a periodic oscillation in the TOAs of the radio beams emitted by pulsars, which can then be collected in PTA telescopes.PTAs stand out as excellent laboratories to test the effects of ULDM in the mass range m ϕ ∼ [10 −24 eV, 10 −22 eV].In this work, we showed that PTAs constrain the presence of ULDM below a few tenths of the observed DM abundance in the mass range m ϕ ∼ [10 −24 eV, 10 −23.3 eV].Therefore, in this range, ULDM cannot constitute 100% of the observed DM.
We wish to thank Diego Blas, Alessio Zicoschi and Paolo Salucci for useful discussions and insights.Moreover, we thank the anonymous referees for improving the quality of our manuscript.

8 log
10 f(Hz) -24.0 -23.8 -23.6 -23.4 -23.2 -23.0 -22.8 -22.6 log 10 m φ (eV) C.S. and E.B. acknowledge support from the European Union's H2020 ERC Consolidator Grant "GRavity from Astrophysical to Microscopic Scales" (Grant No. GRAMS-815673) and the EU Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No. 101007855.B.G. is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN.The European Pulsar Timing Array (EPTA) is a collaboration between European and partner institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), the University of Manchester (UK), the University of Birmingham (UK), the University of East Anglia (UK), the University of Bielefeld (GER), the University of Paris (FRA), the University of Milan-Bicocca (IT), the Foundation for Research and Technology, Hellas (GR), and Peking University (CHN), with the aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves.An Advanced Grant of the European Research Council allowed to implement the Large European Array for Pulsars (LEAP) under Grant Agreement Number 227947 (PI M. Kramer).The EPTA is part of the International Pulsar Timing Array (IPTA); we thank our IPTA colleagues for their support and help with this paper and the external Detection Committee members for their work on the Detection Checklist.The work presented here is a culmination of many years of data analysis as well as software and instrument development.In particular, we thank Drs. N. D'Amico, P. C. C. Freire, R. van Haasteren, C. Jordan, K. Lazaridis, P. Lazarus, L. Lentati, O. Löhmer and R. Smits for their past contributions.We also thank Dr. N. Wex for supporting the calculations of the galactic acceleration as well as the related discussions.The EPTA is also grateful to staff at its observatories and telescopes who have made the continued observations possible.Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany.Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK's Science and Technology Facilities Council (STFC).ICN is also supported by the STFC doctoral training grant ST/T506291/1.The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France.We acknowledge financial support from "Programme National de Cosmologie and Galaxies" (PNCG), and "Programme National Hautes Energies" (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France.We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France.The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO).The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as a National Facility by the National Institute for Astrophysics (INAF).The work is supported by the National SKA programme of China (2020SKA0120100), Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific.This work is also supported as part of the "LEGACY" MPG-CAS collaboration on low-frequency gravitational wave astronomy.J.A. acknowledges support from the European Commission (Grant Agreement number: 101094354).J.A. and S. Chanlaridis were partially supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for