Stochastic Ultralight Dark Matter Fluctuations in Pulsar Timing Arrays

Metric perturbations induced by ultralight dark matter (ULDM) fields have long been identified as a potential target for pulsar timing array (PTA) observations. Previous works have focused on the coherent oscillation of metric perturbations at the characteristic frequency set by the ULDM mass. In this work, we show that ULDM fields source low-frequency stochastic metric fluctuations and that these low-frequency fluctuations can produce distinctive detectable signals in PTA data. Using the NANOGrav 12.5-year data set and synthetic data sets mimicking present and future PTA capabilities, we show that the current and future PTA observations provide the strongest probe of ULDM density within the solar system for masses in the range of $10^{-18}\;{\rm eV}-10^{-16}\;{\rm eV}$.

Millisecond pulsars are one of nature's most precise clocks, with a frequency stability that rivals one of human-made atomic clocks.By monitoring the ticking of a net of these cosmic clocks, pulsar timing arrays (PTAs) have achieved unparalleled sensitivity to spacetime perturbations in the nHz band.By leveraging this remarkable sensitivity, several PTA collaborations have recently found evidence for the presence of a gravitational wave background permeating our galaxy [1][2][3][4].
While the primary objective of PTAs is the detection of gravitational waves (GWs), we can also leverage their sensitivity to search for spacetime fluctuations produced by different sources.A notable example is the metric fluctuations induced by ultralight dark matter-a bosonic dark matter candidate with a mass smaller than a few eV.Such light bosonic dark matter candidates behave like a classical wave, coherently oscillating at a frequency set by their mass, m ϕ .These oscillations will cause time-dependent metric fluctuations at a frequency ω = 2m ϕ , making ULDM candidates in the mass window 10 −23 eV − 10 −20 eV potentially detectable by PTAs, as first demonstrated by Khmelnitsky and Rubakov [5].In addition, if ULDM couples to the Standard Model particles non-gravitationally, it will affect the spin of pulsars as well as the frequency of atomic clocks used for pulsar timing measurements, leading to another monochromatic signal in the timing observation [6][7][8].
However, ULDM density fluctuations exhibit much richer temporal and spectral behaviors than a simple coherent oscillation.In a recent work [9], it was shown that, in addition to the coherent oscillation at ω = 2m ϕ , ULDM fields exhibit low-frequency stochastic fluctuations at ω ≲ m ϕ σ 2 , where σ ≃ 160 km/sec is the velocity dispersion of the virialized dark matter halo.These lowfrequency fluctuations are intricately related to the orderone density fluctuations within the ULDM halo observed in numerical simulations (e.g.Ref. [10]) over characteristic time and length scales given by the coherence time, τ = 1/m ϕ σ 2 , and the coherence length, λ = 1/m ϕ σ.Implications of these low-frequency stochastic fluctuations in PTA observations have not yet been investigated and will be the primary focus of this work.
We can intuitively understand these low-frequency fluctuations in terms of quasiparticles, whose size is set by the coherence length and whose mass is given by m eff = π 3/2 ρ DM λ 3 , with ρ DM being the ULDM density [11,12].In Ref. [9], it was shown that the gravitational interaction between these quasiparticles and test masses in GW interferometers inevitably introduces stochastic signals, potentially detectable by existing and planned future GW detectors.It was also found that interferometers with longer baselines are better suited to look for such stochastic signals-since they are sensitive to larger quasiparticles which can impart larger accelerations on the test masses.In the same study, it was also highlighted that more than two detectors would be required to distinguish stochastic ULDM signals from detector noise.For these reasons, PTAs arise as a natural candidate to look for such stochastic ULDM fluctuations.Thanks to the O(kpc) distance between pulsars and the Earth, PTAs can be perturbed by larger quasiparticles compared to test masses in Earth-based interferometers.Moreover, correlations of the ULDM signal across the pulsars in the array provide a natural means to distinguish ULDM signals from detector-specific noise.
In the quasiparticle picture, one can easily estimate the perturbation to the TOAs induced by stochastic ULDM fluctuations.For PTAs with an observing baseline, T obs , larger than the ULDM coherence times, a typical TOA fluctuation is approximately given by: where ā = Gm eff /λ 2 is a typical acceleration on Earth produced by ULDM quasiparticles, c ∼ O(10 −2 -10 −1 ) is a numerical coefficient, and ρ ⊕ is an average ULDM density around the solar system.This expression is valid for T obs /τ ≫ 1, and we choose T obs = 30 yr for the estimation.In the expression, the quantity in the first parenthesis (āτ 2 ) can be understood as a typical displacement of Earth due to quasiparticles over one coherence time scale and the quantity (T obs /τ ) 3/2 can be understood as a result of random walk motion of the relative velocity of Earth-pulsar system induced by ULDM quasiparticles.
A naive comparison of this pulsar timing fluctuation with typical TOA errors σ TOA ∼ O(10 −7 ) sec already suggests that the current PTA might probe average solar-system ULDM densities of ρ ⊕ /ρ 0 ∼ 10 3 at m ϕ = 10 −17 eV, where ρ 0 ≃ 0.4 GeV/cm 3 represents typical local dark matter density.It is important to note that this local density, ρ 0 , is measured on volumes of at least ∼ (10 2 pc) 3 surrounding the solar system (see for example Refs.[13,14]), while there is currently no direct measurement of dark matter density within and around the solar system.This naive estimate suggests that PTAs might provide one of the strongest probes of ULDM density near the solar system.The main objective of this work is to sharpen the above estimate by carefully characterizing the effects of stochastic ULDM fluctuations in current and future PTA observations.We summarize our main results in Figure 1 where we show the constraints and projections on the ULDM density near the solar system derived by analyzing the NANOGrav 12.5-year data set [15], and a collection of simulated data sets mimicking the capabilities of future observations (see Section III C for details on the generation of the mock data).A null result from the analysis of NANOGrav 12.5-year data sets the strongest upper limit on ULDM density near the solar system as ρ ⊕ /ρ 0 ≲ 4 × 10 3 , improving the existing constraints from solar system ephemerides ρ ⊕ /ρ 0 ≲ 2 × 10 4 [16] in the mass range 10 −18 eV − 10 −16 eV.The analyses with simulated data sets show that current constraints could be improved by an order of magnitude or more with a larger number of pulsars with a longer observation period, which could be easily achieved in future PTAs.
The paper is organized as follows.In section II, we discuss the statistical property of the ULDM field.We compute the timing residual power spectrum induced by stochastic ULDM fluctuations and derive their correlation pattern.In section III, we detail our analysis: we discuss the data analysis tools used in our analysis, and how we generated the mock data sets.In section IV, we present the result from our analysis of the mock pulsar timing data in the context of ultralight dark matter.We conclude in section V.

II. STOCHASTIC ULDM SIGNAL
Pulsar timing arrays track the TOAs of radio pulses emitted by a collection of galactic millisecond pulsars.These TOAs are then compared with the predictions of a pulsar timing model.The discrepancies between these predictions and the observed TOAs define the timing residuals, δt.In this section, we will discuss the imprint of stochastic ULDM fluctuations on these timing residuals.
We parametrize the metric perturbations induced by fluctuations of the ULDM field as ( These metric perturbations can influence the TOAs and leave an imprint in the timing residuals in several different ways: • Doppler effect: The force exerted by ULDM fluctuations induces an acceleration of both Earth and pulsars.This results in perturbations of the perceived pulsar period due to the induced relative motion between the observer (Earth) and the emitters (pulsars).We can write these perturbations in terms of the scalar potential Φ as: where ν a is the perceived pulsar rotational frequency for the a-th pulsar in the array, na is a unit vector pointing from Earth to pulsar, d a is Earth-pulsar distance, and x e and x p represent Earth and pulsar position respectively.
• Other effects: ULDM fluctuations source metric fluctuation between Earth and pulsars, thereby perturbing the light travel time along the line-of-sight.This effect is often called as Shapiro delay.Additionally, metric fluctuations at the position of the Earth and pulsar change the measurements of TOAs and pulsar period as both are measured by the proper time of observer at the Earth and pulsar, respectively.This effect is called Einstein delay.Combined, they induce a shift in the perceived rotational pulsar period as where x(t ′ ) = x e + na (t − t ′ ), i.e. the integral must be performed along the photon path.
Timing residuals can be computed from fluctuations in the arrival frequency as In this work, we focus on the Doppler signal given in Eq. (3) as the Shapiro and Einstein effects are suppressed by an additional factor of dark matter velocity.See Appendix B for more details.

A. ULDM Signal Power Spectrum
The goal of this section is to relate the PTA signals discussed in the previous section with ULDM density fluctuations, δρ, defined as where, as already mentioned before, ρ ⊕ is the average ULDM density within the solar system.Specifically, since we are interested in stochastic ULDM fluctuations we want to relate the power spectrum of ULDM density fluctuations to the power spectrum of PTA timing residuals.The first of these two quantities is defined as where δρ(k) denotes the Fourier components of the ULDM density fluctuations, and the angle brackets denote the ensemble average over all possible realizations of the ULDM field.
Following the procedure outlined in Refs.[9,17,18], we can show that the power spectrum for ULDM density fluctuations can be expressed in terms of the DM phase space distribution, f (p), as (see Appendix A for more details): where we have ignored the modes at ω = ±(ω 1 + ω 2 ) ≃ ±2m ϕ , which represents the coherently oscillating modes not of interest for the analysis discussed in this work.Using Poisson equation, ∇ 2 Φ = 4πGδρ,1 the power spectrum of ULDM density fluctuations can be easily related to the power spectrum of the scalar potential fluctuations, P Φ , as where G is the gravitational constant, and P Φ , is defined analogously to the density power spectrum in Eq. ( 7).
The power spectrum for scalar perturbations can also be related to the timing residuals covariance matrix, defined as: (10) where Γ DM ab is the overlap reduction function (ORF) characterizing the correlations among timing residual of different pulsars, S DM (f ) is the timing residual power spectrum, and the limits of the integration are determined by the inverse of the observational cadence and the total observation period.Here i and j index TOAs, and a and b index pulsars.By using Eq. ( 3) together with Eqs. ( 5) and (10), we can relate S DM ab ≡ Γ DM ab S DM (f ) with fluctuations of the gravitational potential: Given that 2πf d a ≫ 1 for the frequency range of PTA observations and typical Earth-pulsar distances, we can approximate U a U * b ≈ 1 + δ ab .Finally, using the relation in Eq. ( 9), we can rewrite Eq. ( 11) as (12) In the two following sections, we will use this expression to derive the timing residual signal for two specific forms of the DM phase space distribution.

B. Isotropic Distribution
If we ignore the motion of the solar system with respect to the galactic center, the DM has an isotropic velocity distribution that can be written as In this case, we can derive simple analytic expressions for the timing-residual power spectrum and its ORF (see Appendix B for a derivation): Here τ = 1/m ϕ σ2 is the coherence time, ā = Gm eff /λ 2 is the typical acceleration due to ULDM quasiparticles, is the wavelength, and K n (x) are the modified Bessel functions of the second kind.The spectral shape for the ULDM signal in the isotropic limit is shown in blue in the left panel of Figure 2. At frequencies smaller than f ≲ 1/τ , the power spectrum grows like S DM ∝ 1/f 4 ; while for f ≳ 1/τ the power spectrum decays exponentially.In the right panel of the same figure, we compare the dipole structure of the isotropic ULDM signal with the Hellings & Downs correlations expected for a gravitational wave signal.

C. Anisotropic Distribution
The solar system revolution around the galactic center induces a dark matter wind oriented in the opposite direction of the solar system velocity.Because of this effect, the DM velocity distribution in the solar system reference frame is not isotropic and is characterized by a non-zero mean velocity, v 0 : where v 0 = v 0 v0 is the mean DM velocity in the solar system reference frame.In this analysis, we will use v 0 = 232 km/s and v0 = [−0.47,0.41, −0.78] [20,21]. 2  Similarly, the orbital motion of pulsars will induce a DM wind in the pulsars' reference frame.Due to differences in the peculiar velocity, each pulsar will experience its own DM wind velocity.However, given that pulsars used in PTA observation are within a few kpc from the solar system over which rotational velocities are not expected to change substantially, we may assume that they all share the same wind velocity as long as the peculiar velocity of each pulsar is smaller than the DM wind velocity.This approximation will be used to derive the full PTA signal in the anisotropic limit, although it will not affect our final results as those will be derived only using the Earth term and neglecting the contribution from the pulsar term.
Assuming that the DM velocity distribution is given by Eq. ( 16) both in the Earth and pulsars' reference frame, the spectrum and the overlap reduction function are most conveniently written as The correlation functions are given by where the term proportional to δ ab is the pulsar term, while the remaining one gives the Earth term.In the anisotropic limit, inter-pulsar correlations for the ULDM signal are no longer a simple function of the pulsars' angular separation but depend on the positions of the pulsars in the sky.In the right panel of Figure 2, we report the value of Γ A ab for all the pulsar pairs in our 15-year mock data set.It is clear from this figure that, for a fixed angular separation, there is a spread in the level of correlations due to different locations of the pulsars with respect to the orientation of the DM wind.The spectral functions appearing in Eq. ( 17) can be written as for A =⊥, // .In the left panel of Figure 2, we compare the functions S A (f ) (light and dark green lines) with the power spectrum in the isotropic limit.Detailed expressions and derivations for the H A functions are provided in Appendix B.

III. ANALYSIS
In this section, we discuss noise and signal modeling, the data sets used in this analysis, and the statistical tools used to analyze them.

A. Noise Modeling and the PTA Likelihood
The statistical tools needed to model PTA timing residuals have been extensively discussed in the literature (see, e.g., Refs.[22,23]).In this section, we provide a short overview of these tools paying particular attention to how we can use them to model the ULDM signals derived in the previous section.
The pulsar timing residuals receive contributions from several noise and astrophysical sources.All these sources are usually described by using a phenomenological model containing three main components: white noise, timecorrelated stochastic processes (also known as red noise), and the impact of small errors in the fit to the timingephemeris parameters.Specifically, this phenomenological model parametrizes the timing residuals as: where δt is a vector containing all the N TOA measured timing residuals.
The first term on the right-hand side of Eq. ( 21), n, describes the white noise that is assumed to be left in each of the timing residuals after subtracting all known systematics.White noise is assumed to be a zero mean normal random variable, fully characterized by its covariance.Assuming that a single instrument measures all the TOAs using a wideband approach, the white-noise covariance matrix reads: where i and j index the TOAs, a and b index pulsars, σ a,i is the TOA uncertainty for the ith observation, F is the Extra FACtor (EFAC) parameter, and Q is the Extra QUADrature (EQUAD) parameter. 3he second term on the right-hand side of Eq. ( 21) describes time-correlated stochastic processes, which include pulsar-intrinsic red noise, and any stationary stochastic signal.These processes are modeled using a Fourier basis of frequencies f i ≡ i/T obs , where i indexes the harmonics of the basis, and T obs is the timing baseline, extending from the first to the last recorded TOA in the complete PTA data set.Since we are generally interested in processes that exhibit long timescale correlations, this expansion is truncated after, N f , frequency bins.In this work we will set N f = 30 for intrinsic red noise processes, and N f = 14 for GWB and ULDM signals in the mock and NANOGrav 12.5-year data sets.This set of N f sine-cosine pairs evaluated at the different observation times are contained in the Fourier design matrix, F .The Fourier coefficients of this expansion, a, are assumed to be normally distributed random variables with zero mean and a covariance matrix, ⟨aa T ⟩ = ϕ, given by where a and b index the pulsars, and i and j index the frequency harmonics.The first term in Eq. ( 24) describes intrinsic pulsar noise, and is parametrized as The second term in Eq (24) describes the contribution to the timing residuals induced by a GWB.The ORF for this contribution is given by the well-known Hellings-Downs (HD) function [24]: where x ab = (1 − na • nb )/2.In this analysis, we will parametrize the GWB common spectrum as Finally, the last term in Eq. ( 24) describes the stochastic ULDM contribution to the timing residual and is related to the timing residuals' PSD by where ∆f = 1/T obs .In our analysis, we will always consider non-isotropic DM velocity distribution, such that S DM ab (f ) is given by Eq. ( 17).However, since we are interested in probing the ULDM abundance within the solar system, we will neglect the pulsar term for the ULDM signal, which corresponds to neglecting the δ ab terms in the ORFs given in Eqs.(18) and (19).
The third and last term in Eq. ( 21) describes the impact that linear deviations from the initial best-fit values for the m timing-model parameters have on the timing residuals.The design matrix M is a N TOA × m matrix, which contains the partial derivatives of the TOAs with respect to the timing-model parameters (evaluated at the best-fit values), and ϵ is a vector containing the linear offset from the best-fit values of the timing model parameters.
Given this parametrization for the timing residuals, we can use the two-step marginalization of the timing and noise parameters described in Ref. [25] to obtain the PTA likelihood function  where η contains all the model parameters (i.e.F a and Q a plus all the red noise parameters), and K = D + F ϕF T .Here, D = N +M EM T , where N is a diagonal matrix whose non-zero elements are given by Eq. ( 22), and E = ⟨ϵϵ T ⟩ is set to be a diagonal matrix of very large values (10 40 ), which effectively means that we assume flat priors for the parameters in ϵ.

B. Setting Constraints and Projections
The goal of our analysis is to set constraints and projections on the local ULDM density as a function of the ULDM mass.To do this, we make use of Bayes inference, a technique that uses Bayes' rule of conditional probabilities to derive probability distribution for the parameters of the statistical model used to describe the data.In our case, these parameters are white noise parameters (F a and Q a ) and the ones contained in the covariance matrix defined in Eq. ( 24), which include intrinsic red noise parameters (A a and γ A ), GWB parameters (A GWB and γ GWB ), and ULDM parameters (ρ ⊕ , m ϕ , and σ).Given the PTA likelihood of Eq. ( 29), we can use Bayes' theorem to get p(η|δt) is the posterior probability distribution, p(δt|η) is the PTA likelihood function given in Eq. ( 29) and implemented using the ENTERPRISE [26] and enterprise extensions [27] packages, and p(η) are the prior probability distributions for the noise and ULDM parameters reported in Table I.Following the standard practice [28][29][30][31], we first perform a single-pulsar noise analysis to derive the maximum likelihood values for each pulsars' white noise parameters.In this single-pulsar analysis, we model the timing residuals only using white and intrinsic red noise, such that the PTA likelihood will only depend on the parameters F a , Q a , A a , and γ a .GWB and ULDM parameters are not included, such that the covariance matrix in Eq. ( 24) reads ϕ = φ.
After this single-pulsar analysis, we performed a full PTA analysis in which the white noise parameters are fixed to their maximum likelihood value extracted from the single-pulsar analysis, while all the other parameters are varied within the prior ranges summarized in Table I.
Given the posterior distribution, we marginalize over all parameters except the local ULDM density, ρ ⊕ , and the ULDM mass m ϕ .The constraints on the local ULDM in each mass bin are then set to the 95 th percentile of the ULDM density in that bin.
Practically, we derive marginalized posterior distributions by using Markov Chain Monte Carlo (MCMC) techniques to sample randomly from the posterior distributions.To assess convergence of our MCMC runs we use the Gelman-Rubin statistic, R, and require R < 1.02 for all parameters including pulsar intrinsic red noise parameters.All the MCMC runs performed in this analysis used the PTMCMC sampler [32].

NANOGrav 12.5-year data set
The NANOGrav 12.5-year data set [15] consists of observations of 47 millisecond pulsars made between July 2004 and June 2017 by the Arecibo Observatory and the Green Bank Telescope.Most pulsars in this data set were observed approximately once per month, except six pulsars that were observed once per week as part of a highcadence campaign carried out at the Green Bank Telescope since 2013 and at the Arecibo Observatory since 2015.In our analysis, we will use the data from the 45 pulsars of this data set that have an observation baseline longer than 3 years.

Mock data sets
We create our mock data using libstempo [33], a Python wrapper of TEMPO2 [34,35].We build the data set starting from a core catalog consisting of the 67 pulsars contained in the NANOGrav 15-year data set [36], for which we assume an observation baseline of 15 years. 4his core catalog is expanded by adding 33 pulsars every 5 years until a total observing time of 30 years is reached.Shorter observing baselines are derived by slicing this data set into smaller catalogs.
The final data set consists of 166 pulsars placed in random sky locations.Of these pulsars, 67 have an observing  baseline of 30 years, while the remaining 99 are divided into three blocks which have a total observing baseline of 15, 10, and 5 years.For all pulsars, we assume an observing cadence of 3 weeks, with small random fluctuations added on top.Each TOA in the mock data set is associated with a TOA error, σ a,i , which we derive by sampling from a normal distribution whose mean and standard deviation are reported in Table II.The synthetic TOAs are injected with withe and red noise, plus a GWB signal.All these processes are described using the statistical models discussed in Section III A. For each pulsar, the noise parameters are randomly sampled from the distributions given in Table II, while the GWB parameters are set to log 10 A GWB = −14.2and γ GWB = 3.2.

A. Solar system dark matter density
The results of our analysis are summarized in Fig. 1 where we show the constraints and projections on the average ULDM density obtained by analyzing the NANOGrav 12.5-year data set and a simulated data set with increasing observing baseline and number of pulsars.Full posterior distributions for the ULDM and GWB parameters are reported in Appendix C.
The maximum sensitivity is achieved at the ULDM masse, m ϕ ∼ 1/(T obs σ 2 ), for which the ULDM coherence time approximately match the total observing baseline.At lower masses, the sensitivity weakens exponentially, while at higher masses they decrease like ρ ⊕ ∝ m 3/2 ϕ .This is expected as for m ϕ < 1/(T obs σ 2 ) the ULDM signal within the PTA band is exponentially suppressed, while for m ϕ > 1/(T obs σ 2 ) the amplitude of the ULDM signal scales like S DM ∝ ρ 2 ⊕ /m 3 ϕ , while the ULDM signals are entirely within the PTA frequency band.
It is important to stress that PTAs are sensitive to the DM abundance near the solar system, a quantity for which we have no direct experimental measure.Indeed, typical measurements of the DM abundance are derived using large-scale properties of the Milky Way, and probe volumes of order O(10 6 pc 3 ) or larger (for a review on measurements of the local DM abundance see for example Refs.[13,14]).In other words, the result of these measurements, i.e. ρ 0 ∼ 0.4 GeV/cm 3 , do not exclude the possibility of DM overdensities at much smaller scales.
To the best of our knowledge, the strongest constraints on the dark matter density within the solar system come from solar system ephemerides, which constraint the DM density to be smaller than ρ ⊕ /ρ 0 ≲ 2 × 10 5 at the Earth orbit, and ρ ⊕ /ρ 0 ≲ 2×10 4 at the Mars/Saturn orbit [16].Additional constraints arise from lunar laser ranging and LAGEOS geodetic satellite, ρ ⊕ /ρ 0 ≲ 10 11 [37], and the motion of asteroids in the solar system, ρ ⊕ /ρ 0 ≲ 6 × 10 6 [38].Our result shows that PTAs provide the strongest probe on ULDM densities around the solar system in the mass range of m ϕ ∼ 10 −18 eV − 10 −16 eV.

B. Dark matter substructures
We perform a similar setup analysis while, this time, fixing the velocity dispersion to σ = 10, 50, 164 km/sec and assuming vanishing mean velocity v 0 = 0.This analysis is motivated by the potential existence of dark matter substructures near the solar system with distinct kinematic properties.This may include a hypothetical dark disc [39], stream-like structures accompanying stellar streams [40][41][42], or dense ULDM structures bound to the solar system via capture processes [43][44][45].This particular analysis examines the capacity of PTA observations to probe cold dark matter substructures with small velocity dispersion.
The result is shown in Figure 3.We observe two interesting features: (i) as the velocity dispersion decreases, the mass at which the strongest constraint appears shifts as m ϕ ∝ 1/σ 2 and (ii) the constraint at that mass scales as ρ ⊕ ∝ 1/σ.The first feature can be easily explained.From the analytic expression for the power spectrum (15), one finds that the ULDM signal drops exponentially for 2πf ≳ m ϕ σ 2 .This suggests that the lowest mass that can be probed by PTA is mϕ = 2π/(T obs σ 2 ), which explains the first feature shown in the figure .The second feature can be explained similarly.The timing residual power spectrum induced by ULDM fluctuations is proportional to S DM ∝ ρ 2 ⊕ /m 3 σ 4 .For m ϕ = mϕ , we then have S DM ∝ ρ 2 ⊕ σ 2 , which explains the ρ ⊕ ∝ 1/σ scaling for the peak sensitivity.

V. CONCLUSION
We have considered the impacts of ULDM lowfrequency stochastic fluctuations on pulsar timing observations and derived the overlap reduction function and the timing residual power spectrum induced by these fluctuations.We have found that these stochastic fluctuations allow us to probe ULDM in a mass range approximately six orders of magnitude higher than usual searches based on coherent ULDM oscillations.To fully assess the prospect of PTA searches for this kind of fluctuations, we have analyzed (i) the NANOGrav 12.5-year data set and (ii) synthetic data with injected stochastic GWB and noise characteristics resembling actual PTA data sets.From the NANOGrav 12.5-year analysis, we have not found any signals of stochastic ULDM fluctuations, from which we place an upper bound on ULDM density near the solar system as ρ ⊕ /ρ 0 ≲ 4 × 10 3 at m ϕ ≃ 10 −17 eV.From the analyses of the simulated data set, we show that the sensitivity of future PTA observations could be improved up to an order of magnitude with the analysis of upcoming data sets.We have extended this analysis to probe cold dark matter substructures near the solar system and found that PTA will provide a strong probe of such objects for a wide range of mass.

Time delay
We start by parametrizing the metric perturbations produced by ULDM fluctuations in terms of the scalar potentials Ψ and Φ: Now, consider a bunch of photons emitted by a pulsar at space-time coordinates x em = (t em , x p ) and propagating to an observer on Earth that receives them at x obs = (t obs , x e ).The geodesic condition, ds 2 = 0, then requires that where Assuming that the metric perturbations are small (i.e.Φ, Ψ ≪ 1), we can approximate the integral on the right-hand side of this equation as where, in going from the first to the second line, we have used the approximation t ′ em,obs ≃ t em,obs + T a , where T a is the rotational period of the a-th pulsar; while, in going from the second to the last line, we have redefined the integration variable as t → t + T a , and used that for the unperturbed photon trajectory, x(t) = x e + (t obs − t)n a , we have that x ′ (t + T a ) ≃ x(t).Therefore, using Eqs.(B2)-(B4) we can write where in rewriting the integration limits we have used that t em ≃ t obs − d a , where d a is the nominal distance between the observer and the a-th pulsar.The left-hand side of this equation is related to the proper time interval between the detection of the two photon bunches for an observer on Earth as ∆τ ≃ t ′ obs t obs dt 1 + Φ(t, x e ) ≃ t ′ obs − t obs + T a Φ(t obs , x e ) (B6) where to evaluate the integral we have used that the scalar potential is approximately constant on the timescale of the pulsar period.In the rest of this section, we will discuss each contribution on the right-hand side of Eq. (B5).
The difference between the two emission times, t ′ em − t em , can be rewritten in terms of the pulsar period by requiring that the proper time interval between two emissions in the pulsar reference frame must equal the pulsar period: B7) where to evaluate the integral we have used that the scalar potential is approximately constant on the timescale of the pulsar period.
The difference ∆ ′ − ∆ represents the change in path lengths the two consecutive photon bunches have to traverse.Indeed, fluctuations in the scalar field potential can accelerate the pulsar or the Earth and change their relative position: (B8) where we have once again used that the scalar potential is approximately constant over a pulsar rotation and that t em ≃ t obs − d a .
The integral in Eq. (B5) can be rewritten by expanding at first order the time dependence of the scalar potential: Substituting Eqs.(B6)-(B11) into Eq.(B5), we find that for an observer on Earth, the interval between the arrival time of the two radiation pulses is given by ∆τ where the three terms contributing to deviations in the observed pulsar's rotation period are given by where in all these expressions we have replaced t obs with t.The first two contributions come from the combination of Shapiro delay and proper time correction, while the last term is due to the Doppler shift of the arrival times caused by the Earth-pulsar relative motion.It is useful to rewrite these quantities in Fourier space: It is worthwhile to compare our result with Khmelnitsky and Rubakov [5].In their work, they have considered a coherently oscillating mode at ω = 2m ϕ .From (B16)-(B17), it is clear that the coherently oscillating signal is strongest in (∆T a /T ) 1 , while the other term is suppressed by the additional power of DM velocity.Using (B16), it is straightforward to reproduce the results of Khmelnitsky and Rubakov as shown in the Appendix of Ref. [9].

Spectrum
As discussed in the main text, the timing residual power spectrum induced by stochastic ULDM fluctuations can be written as where the density fluctuation power spectrum, P δρ , is given in (A14).For an explicit computation, we set the coordinate as v0 = (0, 0, 1), (B19) k = (sin θ cos ϕ, sin θ sin ϕ, cos θ), (B20) na = (sin θ a cos ϕ a , sin θ a sin ϕ a , cos θ a ), (B21) nb = (sin θ b cos ϕ b , sin θ b sin ϕ b , cos θ b ). (B22) After performing an explicit computation, one finds that the timing residual power spectrum can be written as a sum of two terms: ab S / / (f ), (B23) where the overlap functions are given in ( 18)- (19).Each spectrum is given as where we have introduced V 0 = v 0 /σ and ω = ωτ and changed the integration variable to x = v k /σ.Here the function C ⊥ and C// are given as In the limit V 0 → 0, one finds H ⊥ (f ) = H / / (f ) = (64/3π)K 0 (ωτ ), reproducing the results obtained in the isotropic limit.For the Bayesian analysis in the main text, we tabulate each function H A and use them for the numerical evaluations.

10 − 19 12 m
FIG.1.The red line shows the 95% constraints on the dark matter density near the solar system derived by analyzing the NANOGrav 12.5-year data set.The blue lines show projections derived by analyzing mock data sets with Npsr = 66 pulsar and a total observing time of T obs = 15 yr (dark blue), Npsr = 100 pulsars and T obs = 20 yr observing time (medium blue), and Npsr = 166 pulsars and T obs = 30 yr observing time (light blue).These constraints are compared with the ones derived from solar system ephemerides (black lines) at Earth, Mars, and Saturn orbit.All the constraints are normalized to ρ0 = 0.4 GeV/cm 3 .

FIG. 2 .
FIG. 2. Left:Comparison between the sensitivity curve for the NANOGrav 11-year data set[19] (red), the characteristic strain for the ULDM signal in the isotropic limit (blue), and the two components of the ULDM strain in the anisotropic limit (light and dark green).The shaded bands are derived by varying the value of the velocity dispersion between σ = 150 km/s and σ = 250 km/s.In this plot we have chosen ρ⊕/ρ0 = 10 4 , and m ϕ = 10 −17 eV.Right: The blue line shows the overlap reduction function for the ULDM signal in the isotropic limit.Light and dark green dots represent the correlations amongst the pulsars in our 15-year mock data set for the ULDM signal in the anisotropic limit.The dashed black line represents the Hellings & Downs correlation pattern.

TABLE I .
Prior distributions for the parameters used in this work.The parameters marked with a * were used only in single pulsar analyses.

TABLE II .
Probability distributions used to generate values of noise parameters.