# Stochastic Ultralight Dark Matter Fluctuations in Pulsar Timing Arrays

Hyungjin Kim  0000-0002-8843-7690 Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany  Andrea Mitridate  0000-0003-2898-5844
Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

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}$.

††preprint: DESY-23-192

###### Contents

  1. I Introduction
  2. II Stochastic ULDM signal
    1. II.1 ULDM Signal Power Spectrum
    2. II.2 Isotropic Distribution
    3. II.3 Anisotropic Distribution
  3. III Analysis
    1. III.1 Noise Modeling and the PTA Likelihood
    2. III.2 Setting Constraints and Projections
    3. III.3 Data Sets
      1. III.3.1 NANOGrav 12.5-year data set
      2. III.3.2 Mock data sets
  4. IV Results
    1. IV.1 Solar system dark matter density
    2. IV.2 Dark matter substructures
  5. V Conclusion
  6. A ULDM fluctuations
  7. B Timing residual
    1. B.1 Time delay
    2. B.2 Spectrum
  8. C Posterior distribution

##  I Introduction

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 ${\rm
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_{\phi}$. These oscillations will cause time-dependent metric
fluctuations at a frequency $\omega=2m_{\phi}$, making ULDM candidates in the
mass window $10^{-23}\,{\rm eV}-10^{-20}\,{\rm 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].

Figure 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
$N_{\rm psr}=66$ pulsar and a total observing time of $T_{\rm obs}=15\,{\rm
yr}$ (dark blue), $N_{\rm psr}=100$ pulsars and $T_{\rm obs}=20\,{\rm yr}$
observing time (medium blue), and $N_{\rm psr}=166$ pulsars and $T_{\rm
obs}=30\,{\rm 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
$\rho_{0}=0.4\,{\rm GeV/cm}^{3}.$

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 $\omega=2m_{\phi}$,
ULDM fields exhibit low-frequency stochastic fluctuations at $\omega\lesssim
m_{\phi}\sigma^{2}$, where $\sigma\simeq 160\,{\rm km/sec}$ is the velocity
dispersion of the virialized dark matter halo. These low-frequency
fluctuations are intricately related to the order-one 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,
$\tau=1/m_{\phi}\sigma^{2}$, and the coherence length,
$\lambda=1/m_{\phi}\sigma$. 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_{\rm eff}=\pi^{3/2}\rho_{\scriptscriptstyle\rm DM}\lambda^{3}$,
with $\rho_{\scriptscriptstyle\rm 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 $\mathcal{O}(\rm{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_{\rm obs}$, larger than the ULDM coherence times, a typical TOA
fluctuation is approximately given by:

$\displaystyle\delta t$ $\displaystyle\simeq
c\,(\bar{a}\tau^{2})\Big{(}\frac{T_{\rm obs}}{\tau}\Big{)}^{\frac{3}{2}}$
$\displaystyle\sim 10^{-10}\,{\rm sec}\left(\frac{10^{-17}\,{\rm
eV}}{m_{\phi}}\right)^{\frac{3}{2}}\left(\frac{\rho_{\oplus}}{0.4\,{\rm
GeV}/{\rm cm}^{3}}\right)$ (1)

where $\bar{a}=Gm_{\rm eff}/\lambda^{2}$ is a typical acceleration on Earth
produced by ULDM quasiparticles, $c\sim{\cal
O}(10^{-2}\,\textrm{--}\,10^{-1})$ is a numerical coefficient, and
$\rho_{\oplus}$ is an average ULDM density around the solar system. This
expression is valid for $T_{\rm obs}/\tau\gg 1$, and we choose $T_{\rm
obs}=30\,{\rm yr}$ for the estimation. In the expression, the quantity in the
first parenthesis ($\bar{a}\tau^{2}$) can be understood as a typical
displacement of Earth due to quasiparticles over one coherence time scale and
the quantity $(T_{\rm obs}/\tau)^{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
$\sigma_{\rm TOA}\sim{\cal O}(10^{-7})\,{\rm sec}$ already suggests that the
current PTA might probe average solar-system ULDM densities of
$\rho_{\oplus}/\rho_{0}\sim 10^{3}$ at $m_{\phi}=10^{-17}\,{\rm eV}$, where
$\rho_{0}\simeq 0.4\;{\rm GeV}/{\rm cm}^{3}$ represents typical local dark
matter density. It is important to note that this local density, $\rho_{0}$,
is measured on volumes of at least $\sim(10^{2}\;{\rm 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.3 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 $\rho_{\oplus}/\rho_{0}\lesssim 4\times 10^{3}$, improving
the existing constraints from solar system ephemerides
$\rho_{\oplus}/\rho_{0}\lesssim 2\times 10^{4}$ [16] in the mass range
$10^{-18}\;{\rm eV}-10^{-16}\;{\rm 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_ , $\delta 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

$ds^{2}=\Big{[}1+2\Phi(t,\bm{x})\Big{]}dt^{2}-\Big{[}1-2\Psi(t,\bm{x})\Big{]}dx^{2}\,.$
(2)

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 $\Phi$ as:

$\frac{\delta\nu_{a}(t)}{\nu_{a}}=\hat{\bm{n}}_{a}\cdot\int^{t}dt^{\prime}\;\nabla\Big{[}\Phi(t^{\prime}-d_{a},\bm{x}_{p})-\Phi(t^{\prime},\bm{x}_{e})\Big{]}\,,$
(3)

where $\nu_{a}$ is the perceived pulsar rotational frequency for the $a$-th
pulsar in the array, $\hat{\bm{n}}_{a}$ is a unit vector pointing from Earth
to pulsar, $d_{a}$ is Earth-pulsar distance, and $\bm{x_{e}}$ and $\bm{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

$\begin{split}\frac{\delta\nu_{a}(t)}{\nu_{a}}=&\Psi(\bm{x}_{e},t)-\Psi(\bm{x}_{p},t-d_{a})\\\
-&\hat{\bm{n}}_{a}\cdot\int_{t-d_{a}}^{t}dt^{\prime}\;\nabla\Big{[}\Phi(t^{\prime},\bm{x}(t^{\prime}))+\Psi(t^{\prime},\bm{x}(t^{\prime}))\Big{]}\end{split}$
(4)

where $\bm{x}(t^{\prime})=\bm{x}_{e}+\hat{\bm{n}}_{a}(t-t^{\prime})$, i.e. the
integral must be performed along the photon path.

Timing residuals can be computed from fluctuations in the arrival frequency as

$\delta
t_{a}(t)=\int^{t}\frac{\delta\nu_{a}(t^{\prime})}{\nu_{a}}dt^{\prime}\,.$ (5)

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.

###  II.1 ULDM Signal Power Spectrum

The goal of this section is to relate the PTA signals discussed in the
previous section with ULDM density fluctuations, $\delta\rho$, defined as

$\delta\rho\equiv\rho-\rho_{\oplus}\,,$ (6)

where, as already mentioned before, $\rho_{\oplus}$ 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

$\langle\widetilde{\delta\rho}(k)\widetilde{\delta\rho}^{*}(k^{\prime})\rangle=(2\pi)^{4}\delta^{(4)}(k-k^{\prime})P_{\delta\rho}(k)\,.$
(7)

where $\tilde{\delta}\rho(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. [17, 9, 18], we can show that the
power spectrum for ULDM density fluctuations can be expressed in terms of the
DM phase space distribution, $f(\bm{p})$, as (see Appendix A for more
details):

$\begin{split}P_{\delta\rho}\approx
m_{\phi}^{2}\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\frac{d^{3}p_{2}}{(2\pi)^{3}}f(\bm{p}_{1})f(\bm{p}_{2})(2\pi)^{4}\delta^{4}(k-p_{1}+p_{2}),\end{split}$
(8)

where we have ignored the modes at $\omega=\pm(\omega_{1}+\omega_{2})\simeq\pm
2m_{\phi}$, which represents the coherently oscillating modes not of interest
for the analysis discussed in this work. Using Poisson equation,
$\nabla^{2}\Phi=4\pi G\delta\rho$,111For the low-frequency stochastic
fluctuations, one can approximate $\Psi\approx\Phi$ [9]. the power spectrum of
ULDM density fluctuations can be easily related to the power spectrum of the
scalar potential fluctuations, $P_{\Phi}$, as

$P_{\Phi}(k)=\frac{(4\pi G)^{2}}{k^{4}}P_{\delta\rho}(k)\,,$ (9)

where $G$ is the gravitational constant, and $P_{\Phi}$, 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:

$\langle\delta t_{ai}\delta
t_{bj}\rangle\supset\Gamma_{ab}^{\scriptscriptstyle\rm DM}\int df\cos[2\pi
f(t_{ai}-t_{bj})]\,S^{\scriptscriptstyle\rm DM}(f)$ (10)

where $\Gamma_{ab}^{\scriptscriptstyle\rm DM}$ is the overlap reduction
function (ORF) characterizing the correlations among timing residual of
different pulsars, $S^{\scriptscriptstyle\rm 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^{\scriptscriptstyle\rm
DM}_{ab}\equiv\Gamma_{ab}^{\scriptscriptstyle\rm DM}S^{\scriptscriptstyle\rm
DM}(f)$ with fluctuations of the gravitational potential:

$S^{\scriptscriptstyle\rm DM}_{ab}(f)=\frac{2}{(2\pi
f)^{4}}\int\frac{d^{3}k}{(2\pi)^{3}}(\hat{n}_{a}\cdot\bm{k})(\hat{n}_{b}\cdot\bm{k})U_{a}U_{b}^{*}P_{\Phi}(f,\bm{k}),$
(11)

where $U_{a}=1-\exp\big{(}2\pi
ifd_{a}[1+(k/\omega)\hat{k}\cdot\hat{n}_{a}]\big{)}$. Given that $2\pi
fd_{a}\gg 1$ for the frequency range of PTA observations and typical Earth-
pulsar distances, we can approximate $U_{a}U_{b}^{*}\approx 1+\delta_{ab}$.

Finally, using the relation in Eq. (9), we can rewrite Eq. (11) as

$S^{\scriptscriptstyle\rm
DM}_{ab}(f)=\frac{2G^{2}}{\pi^{2}f^{4}}\int\frac{d^{3}k}{(2\pi)^{3}}(\hat{n}_{a}\cdot\bm{k})(\hat{n}_{b}\cdot\bm{k})U_{a}U_{b}^{*}\frac{P_{\delta\rho}(f,\bm{k})}{k^{4}}.$
(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.

Figure 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 $\sigma=150\;{\rm km/s}$
and $\sigma=250\;{\rm km/s}$. In this plot we have chosen
$\rho_{\oplus}/\rho_{0}=10^{4}$, and $m_{\phi}=10^{-17}\;{\rm 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.

###  II.2 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

$f(\bm{v})=\frac{\rho_{\oplus}/m_{\phi}}{(2\pi\sigma^{2})^{3/2}}\exp\left[-\frac{v^{2}}{2\sigma^{2}}\right].$
(13)

In this case, we can derive simple analytic expressions for the timing-
residual power spectrum and its ORF (see Appendix B for a derivation):

$\displaystyle\Gamma_{ab}^{\scriptscriptstyle\rm DM}$
$\displaystyle=\frac{1}{2}\big{(}\delta_{ab}+\hat{n}_{a}\cdot\hat{n}_{b}\big{)},$
(14) $\displaystyle S^{\scriptscriptstyle\rm DM}(f)$
$\displaystyle=\frac{\bar{a}^{2}\tau}{(2\pi
f)^{4}}\left[\frac{64}{3\pi}K_{0}(\omega\tau)\right].$ (15)

Here $\tau=1/m_{\phi}\sigma^{2}$ is the coherence time, $\bar{a}=Gm_{\rm
eff}/\lambda^{2}$ is the typical acceleration due to ULDM quasiparticles,
$m_{\rm eff}=\pi^{3/2}\rho_{\oplus}\lambda^{3}$ is the quasiparticle mass,
$\lambda=1/m_{\phi}\sigma^{2}$ 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\lesssim 1/\tau$,
the power spectrum grows like $S^{\scriptscriptstyle\rm DM}\propto 1/f^{4}$;
while for $f\gtrsim 1/\tau$ 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.

###  II.3 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,
$\bm{v}_{0}$:

$f(\bm{v})=\frac{\rho_{\oplus}/m_{\phi}}{(2\pi\sigma^{2})^{3/2}}\exp\left[-\frac{(\bm{v}-\bm{v}_{0})^{2}}{2\sigma^{2}}\right]\,,$
(16)

where $\bm{v}_{0}=v_{0}\hat{v}_{0}$ is the mean DM velocity in the solar
system reference frame. In this analysis, we will use $v_{0}=232\;{\rm
km}/{\rm s}$ and $\hat{v}_{0}=\left[-0.47,0.41,-0.78\right]$ [20, 21].222Here
we ignore the subleading effect induced by Earth’s revolution around the Sun.
The vector $\hat{v}_{0}$ is given in the equatorial coordinate system.

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

$\Gamma_{ab}^{\scriptscriptstyle\rm DM}S^{\scriptscriptstyle\rm
DM}(f)=\Gamma^{\perp}_{ab}S^{\perp}(f)+\Gamma^{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}_{ab}S^{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}(f)\,.$
(17)

The correlation functions are given by

$\displaystyle\Gamma^{\perp}_{ab}$
$\displaystyle=\frac{1}{2}(1+\delta_{ab})\Big{[}\hat{n}_{a}\cdot\hat{n}_{b}-(\hat{n}_{a}\cdot\hat{v}_{0})(\hat{n}_{b}\cdot\hat{v}_{0})\Big{]},$
(18)
$\displaystyle\Gamma^{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}_{ab}$
$\displaystyle=\frac{1}{2}(1+\delta_{ab})(\hat{n}_{a}\cdot\hat{v}_{0})(\hat{n}_{b}\cdot\hat{v}_{0})\,,$
(19)

where the term proportional to $\delta_{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 $\Gamma_{ab}^{A}$
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

$S^{A}(f)=\frac{\bar{a}^{2}\tau}{(2\pi f)^{4}}H^{A}(f,v_{0})\,,$ (20)

for
$A=\perp,\mathrel{\mathchoice{\leavevmode\resizebox{}{6.83331pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{6.83331pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.41666pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}$.
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.

###  III.1 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, time-
correlated stochastic processes (also known as red noise), and the impact of
small errors in the fit to the timing-ephemeris parameters. Specifically, this
phenomenological model parametrizes the timing residuals as:

$\bm{\delta t}=\bm{n}+\bm{F}\,\bm{a}+\bm{M}\,\bm{\epsilon}\,,$ (21)

where $\bm{\delta t}$ is a vector containing all the $N_{\rm TOA}$ measured
timing residuals.

The first term on the right-hand side of Eq. (21), $\bm{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:

$\langle
n_{a,i}n_{b,j}\rangle=\mathcal{F}_{a}^{2}\left[\sigma_{a,i}^{2}+\mathcal{Q}_{a}^{2}\right]\;\delta_{ij}\delta_{ab}\,,$
(22)

where $i$ and $j$ index the TOAs, $a$ and $b$ index pulsars, $\sigma_{a,i}$ is
the TOA uncertainty for the $i$th observation, $\mathcal{F}$ is the Extra
FACtor (EFAC) parameter, and $\mathcal{Q}$ is the Extra QUADrature (EQUAD)
parameter.333For the analysis of the NANOGrav 12.5-year data, which was
derived with narrowband observations, we also include correlated noise across
radio frequencies as  $\displaystyle\langle
n_{a,i}n_{b,j}\rangle\supset\mathcal{J}^{2}\mathcal{U}_{ij}\,$ (23) where
$\mathcal{J}^{2}$ is the ECORR parameter, and $\mathcal{U}_{ij}$ is a block
diagonal matrix with unit entries for TOAs in the same observing epoch and
zeros for all others entries.

The 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}\equiv i/T_{\textrm{obs}}$, where $i$ indexes the
harmonics of the basis, and $T_{\textrm{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, $\bm{F}$. The
Fourier coefficients of this expansion, $\bm{a}$, are assumed to be normally
distributed random variables with zero mean and a covariance matrix,
$\langle\bm{a}\bm{a}^{\textrm{T}}\rangle=\bm{\phi}$, given by

$[\phi]_{(ai)(bj)}=\delta_{ij}\left(\delta_{ab}\varphi_{a,i}+\Gamma_{ab}^{\scriptscriptstyle\rm
GW}\Phi_{i}^{\scriptscriptstyle\rm GW}+\Gamma_{ab}^{\scriptscriptstyle\rm
DM}\Phi_{i}^{\scriptscriptstyle\rm DM}\right)\,,$ (24)

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

$\varphi_{a,i}=\frac{A_{a}^{2}}{12\pi^{2}}\frac{1}{T_{\rm
obs}}\left(\frac{f_{i}}{{\rm yr}^{-1}}\right)^{\gamma_{a}}{\rm yr}^{3}\,.$
(25)

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]:

$\Gamma_{ab}^{\scriptscriptstyle\rm
GW}=\frac{1}{2}\delta_{ab}+\frac{1}{2}-\frac{1}{4}x_{ab}+\frac{3}{2}x_{ab}\ln
x_{ab},$ (26)

where $x_{ab}=(1-\hat{n}_{a}\cdot\hat{n}_{b})/2$. In this analysis, we will
parametrize the GWB common spectrum as

$\Phi_{i}^{\scriptscriptstyle\rm GW}=\frac{A_{\scriptscriptstyle\rm
GW}^{2}}{12\pi^{2}}\frac{1}{T_{\rm obs}}\left(\frac{f_{i}}{{\rm
yr}^{-1}}\right)^{\gamma_{\scriptscriptstyle\rm GW}}{\rm yr}^{3}\,.$ (27)

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

$\Gamma_{ab}^{\scriptscriptstyle\rm DM}\Phi_{i}^{\scriptscriptstyle\rm
DM}=\Gamma_{ab}^{\scriptscriptstyle\rm DM}S^{\scriptscriptstyle\rm
DM}(f_{i})\Delta f\,,$ (28)

where $\Delta f=1/T_{\rm obs}$. In our analysis, we will always consider non-
isotropic DM velocity distribution, such that $S_{ab}^{\scriptscriptstyle\rm
DM}(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 $\delta_{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_ $\bm{M}$ is a
$N_{\rm TOA}\times m$ matrix, which contains the partial derivatives of the
TOAs with respect to the timing-model parameters (evaluated at the best-fit
values), and $\bm{\epsilon}$ 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

$p(\bm{\delta t}|\bm{\eta})=\frac{\exp\left(-\frac{1}{2}\bm{\delta
t}^{T}\bm{K}^{-1}\bm{\delta t}\right)}{\sqrt{{\rm det}(2\pi\bm{K})}}\,,$ (29)

where $\bm{\eta}$ contains all the model parameters (i.e. $\mathcal{F}_{a}$
and $\mathcal{Q}_{a}$ plus all the red noise parameters), and
$\bm{K}=\bm{D}+\bm{F}\bm{\phi}\bm{F}^{T}$. Here, $\bm{D}=\bm{N}+\bm{MEM}^{T}$,
where $\bm{N}$ is a diagonal matrix whose non-zero elements are given by Eq.
(22), and $\bm{E}=\langle\bm{\epsilon\epsilon}^{T}\rangle$ 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 $\bm{\epsilon}$.

Parameter | Description | Prior
---|---|---
$\mathcal{F}_{a}$ |  EFAC∗ | $[0.01,\,10]$
$\log_{10}(\mathcal{Q}_{a}/{\rm s})$ |  EQUAD∗ | $[-8.5,\,-5]$
$\log_{10}A_{a}$ | intrinsic red-noise amplitude | $[-20,\,-11]$
$\gamma_{a}$ | intrinsic red-noise spectral index | $[0,\,7]$
$\log_{10}A_{\scriptscriptstyle\rm GW}$ | GWB signal amplitude | $[-20,\,-11]$
$\gamma_{\scriptscriptstyle\rm GW}$ | GWB signal spectral index | $[0,\,7]$
$\log_{10}(\rho_{\oplus}/\rho_{0})$ | normalized local DM density | $[0,\,12]$
$\log_{10}(m_{\phi}/{\rm eV})$ | DM mass | $[-19,\,-12]$
$\sigma\;[{\rm km/s}]$ | DM velocity dispersion | $[150,\,250]$
Table 1: Prior distributions for the parameters used in this work. The
parameters marked with a ∗ were used only in single pulsar analyses.

###  III.2 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 ($\mathcal{F}_{a}$ and $\mathcal{Q}_{a}$) and the ones contained in
the covariance matrix defined in Eq. (24), which include intrinsic red noise
parameters ($A_{a}$ and $\gamma_{A}$), GWB parameters ($A_{\rm GWB}$ and
$\gamma_{\rm GWB}$), and ULDM parameters ($\rho_{\oplus}$, $m_{\phi}$, and
$\sigma$). Given the PTA likelihood of Eq. (29), we can use Bayes’ theorem to
get

$p(\bm{\eta}|\bm{\delta t})\propto p(\bm{\delta t}|\bm{\eta})p(\bm{\eta})$
(30)

$p(\bm{\eta}|\bm{\delta t})$ is the posterior probability distribution,
$p(\bm{\delta t}|\bm{\eta})$ is the PTA likelihood function given in Eq. (29)
and implemented using the ENTERPRISE [26] and enterprise_extensions [27]
packages, and $p(\bm{\eta})$ are the prior probability distributions for the
noise and ULDM parameters reported in Table 1.

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 $\mathcal{F}_{a}$,
$\mathcal{Q}_{a}$, $A_{a}$, and $\gamma_{a}$. GWB and ULDM parameters are not
included, such that the covariance matrix in Eq. (24) reads
$\bm{\phi}=\bm{\varphi}$.

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 1.

Given the posterior distribution, we marginalize over all parameters except
the local ULDM density, $\rho_{\oplus}$, and the ULDM mass $m_{\phi}$. The
constraints on the local ULDM in each mass bin are then set to the $95^{\rm
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].

###  III.3 Data Sets

####  III.3.1 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 high-cadence 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.

####  III.3.2 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.444The NANOGrav 15-year data set contains
in total of 68 pulsars. Our core catalog is based on the 67 that were used in
the NANOGrav most recent GWB search [1] This 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, $\sigma_{a,i}$, which we derive by sampling from
a normal distribution whose mean and standard deviation are reported in Table
2. 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.1. For each pulsar, the noise parameters are randomly
sampled from the distributions given in Table 2, while the GWB parameters are
set to $\log_{10}A_{\rm GWB}=-14.2$ and $\gamma_{\rm GWB}=3.2$.

Parameter | Description | Distribution
---|---|---
$\sigma_{a,i}$ | TOA uncertainty | $\mathcal{N}(400\,{\rm ns},200\,{\rm ns})$
$\mathcal{F}_{a}$ | EFAC | $\mathcal{N}(1,0.05)$
$\log_{10}\mathcal{Q}_{a}$ | EQUAD | $\mathcal{N}(-8.5,1)$
$\log_{10}A_{a}$ | intrinsic red-noise amplitude | $\mathcal{N}(-16,1)$
$\gamma_{a}$ | intrinsic red-noise spectral index | $\mathcal{U}(1,5)$
Table 2: Probability distributions used to generate values of noise
parameters.

##  IV Results

###  IV.1 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_{\phi}\sim 1/(T_{\rm
obs}\sigma^{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 $\rho_{\oplus}\propto
m_{\phi}^{3/2}$. This is expected as for $m_{\phi}<1/(T_{\rm obs}\sigma^{2})$
the ULDM signal within the PTA band is exponentially suppressed, while for
$m_{\phi}>1/(T_{\rm obs}\sigma^{2})$ the amplitude of the ULDM signal scales
like $S^{\scriptscriptstyle\rm DM}\propto\rho_{\oplus}^{2}/m_{\phi}^{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
$\mathcal{O}(10^{6}\;{\rm 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. $\rho_{0}\sim 0.4\;{\rm 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 $\rho_{\oplus}/\rho_{0}\lesssim
2\times 10^{5}$ at the Earth orbit, and $\rho_{\oplus}/\rho_{0}\lesssim
2\times 10^{4}$ at the Mars/Saturn orbit [16]. Additional constraints arise
from lunar laser ranging and LAGEOS geodetic satellite,
$\rho_{\oplus}/\rho_{0}\lesssim 10^{11}$ [37], and the motion of asteroids in
the solar system, $\rho_{\oplus}/\rho_{0}\lesssim 6\times 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_{\phi}\sim 10^{-18}\;{\rm
eV}-10^{-16}\;{\rm eV}$.

###  IV.2 Dark matter substructures

We perform a similar setup analysis while, this time, fixing the velocity
dispersion to $\sigma=10,\,50,\,164\,{\rm 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_{\phi}\propto 1/\sigma^{2}$ and (ii) the constraint at
that mass scales as $\rho_{\oplus}\propto 1/\sigma$.

Figure 3: Constraints and projections on the local ULDM abundance for cold DM
substructures with $\sigma=10\,{\rm km/sec}$ (blue shaded regions), $50\,{\rm
km/sec}$ (red shaded regions), and $164\,{\rm km/sec}$ (green shaded regions).
The solid lines show the constraints derived by analyzing the NANOGrav
12.5-year data set. The other shaded regions show the constraints derived by
analyzing the 15-year (dark-shaded), 20-year (medium-shaded), and 30-year
(light-shaded) synthetic data sets.

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\pi f\gtrsim m_{\phi}\sigma^{2}$. This suggests that the lowest mass
that can be probed by PTA is $\hat{m}_{\phi}=2\pi/(T_{\rm obs}\sigma^{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^{\scriptscriptstyle\rm
DM}\propto\rho_{\oplus}^{2}/m^{3}\sigma^{4}$. For $m_{\phi}=\hat{m}_{\phi}$,
we then have $S^{\scriptscriptstyle\rm DM}\propto\rho_{\oplus}^{2}\sigma^{2}$,
which explains the $\rho_{\oplus}\propto 1/\sigma$ scaling for the peak
sensitivity.

##  V Conclusion

We have considered the impacts of ULDM low-frequency 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 $\rho_{\oplus}/\rho_{0}\lesssim 4\times
10^{3}$ at $m_{\phi}\simeq 10^{-17}\,{\rm 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.

###### Acknowledgements.

We thank the members of the NANOGrav timing working group for producing the
NANOGrav 12.5-year data set used in this analysis, and the entire NANOGrav
collaboration for useful comments on a preliminary form of this work. We also
thank Xiao Xue for useful discussions. This work was supported by the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC 2121 Quantum
Universe - 390833306. HK was also supported by the Munich Institute for
Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC-2094 – 390783311.

##  Appendix A ULDM fluctuations

In this appendix, following Ref. [9], we provide a more detailed derivation of
the power spectrum for density and pressure fluctuations of a ULDM scalar
field, $\phi$, minimally coupled to gravity and whose action reads

$S=\int
d^{4}x\sqrt{-g}\,\left[\frac{1}{2}(\partial\phi)^{2}-\frac{1}{2}m_{\phi}^{2}\phi^{2}\right]\,,$
(31)

where $g={\rm det}\,g_{\mu\nu}$ is the determinant of the metric $g_{\mu\nu}$.

The ULDM field can be modeled as the sum of many freely propagating plane
waves [17]:555Here we are neglecting the $O(1\%)$ perturbations that the
potential of the Sun can have on the ULDM wave functions.

$\phi(t,\bm{x})=\sum_{\bm{k}}\frac{1}{\sqrt{2m_{\phi}V}}\left[\alpha_{\bm{k}}e^{-ik\cdot
x}+\alpha_{\bm{k}}^{*}e^{+ik\cdot x}\right]$ (32)

where the sum runs over the momenta of a particle quantized in a cubic box of
volume $V$, and $(\alpha_{\bm{k}},\alpha_{\bm{k}}^{*})$ is a set of complex
random number, whose underlying distribution is given by

$p(\alpha_{\bm{k}})=\frac{1}{\pi
f_{\bm{k}}}\exp\left[-\frac{|\alpha_{\bm{k}}|^{2}}{f_{\bm{k}}}\right].$ (33)

Here $f_{\bm{k}}$ is the mean occupation number, which will be later
interpreted as a phase space distribution of dark matter particles. The above
distribution does not depend on the phase of $\alpha_{\bm{k}}$; the phase
follows the uniform distribution. Heuristic arguments for the distribution can
be found in Refs. [46, 47]. A field theoretical derivation is provided in Ref.
[17].

Given this probability distribution, the ensemble average of any operator
built from the $\phi$ field can be written as

$\langle{\cal
O}\rangle=\int\Big{[}\prod_{\bm{k}}d^{2}\alpha_{\bm{k}}\,p(\alpha_{\bm{k}})\Big{]}{\cal
O}\,.$ (34)

We will now apply these results to study the density and pressure fluctuations
of a ULDM field.

The stress-energy tensor for the ULDM field reads

$T_{\mu\nu}=\partial_{\mu}\phi\,\partial_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\left[(\partial\phi)^{2}-m_{\phi}^{2}\phi^{2}\right]\,,$
(35)

from which we can easily derive the ULDM energy density:

$\rho\equiv
T_{00}=\frac{1}{2}\Big{[}(\dot{\phi})^{2}+(\nabla\phi)^{2}+m_{\phi}^{2}\phi^{2}\Big{]}.$
(36)

ULDM density fluctuations are defined as

$\delta\rho=\rho-\langle\rho\rangle\,,$ (37)

with the mean energy density given by

$\langle\rho\rangle\simeq\frac{m_{\phi}}{V}\sum_{\bm{k}}f_{\bm{k}}\simeq
m_{\phi}\int d^{3}{\bm{v}}f(\bm{v})$ (38)

where in the last step we have taken the continuum limit, changed the
integration variable to the DM velocity, and adjusted accordingly the
normalization of the phase space distribution.

Given the above results, we are now ready to derive the power spectrum of ULDM
density fluctuations, $P_{\delta\rho}$, defined as

$\langle\delta\rho(k)\delta\rho^{*}(k^{\prime})\rangle\equiv(2\pi)^{4}\delta^{(4)}(k-k^{\prime})P_{\delta\rho}(k)\,$
(39)

where we have defined

$\delta\rho(k)\equiv\int d^{4}x\,e^{-ik\cdot x}\,\delta\rho(x)\,.$ (40)

Using Eqs. (32)–(34) together with the relations

$\begin{split}\langle\alpha_{\bm{k}}^{\phantom{\ast}}\alpha_{\bm{k}^{\prime}}^{*}\rangle&=f_{\bm{k}}\delta_{\bm{k}\bm{k}^{\prime}}\\\
\langle\alpha_{\bm{k}}^{\phantom{\ast}}\alpha_{\bm{q}}^{\phantom{\ast}}\alpha_{\bm{k}^{\prime}}^{*}\alpha_{\bm{q}^{\prime}}^{*}\rangle&=f_{\bm{k}}^{\phantom{\ast}}f_{\bm{q}}^{\phantom{\ast}}(\delta^{\phantom{\ast}}_{\bm{k}\bm{k}^{\prime}}\delta^{\phantom{\ast}}_{\bm{q}\bm{q}^{\prime}}+\delta^{\phantom{\ast}}_{\bm{k}\bm{q}^{\prime}}\delta^{\phantom{\ast}}_{\bm{q}\bm{k}^{\prime}})\,,\end{split}$
(41)

we find that

$\displaystyle P_{\delta\rho}(k)\approx m_{\phi}^{2}\int
d^{3}v_{1}d^{3}v_{2}\,f(\bm{v}_{1})f(\bm{v}_{1})(2\pi)^{4}\delta(k-p_{1}+p_{2}).$
(42)

Given the DM velocity distribution in the solar system:

$f(\bm{v})=\frac{\rho_{\oplus}/m}{(2\pi\sigma^{2})^{3/2}}\exp\left[-\frac{(\bm{v}-\bm{v}_{0})^{2}}{2\sigma^{2}}\right],$
(43)

with $\bm{v}_{0}$ the mean dark matter velocity, $\rho_{\oplus}$ the mean DM
density near the solar system, and $\sigma$ the velocity dispersion; the power
spectrum for density perturbations takes the form

$\displaystyle
P_{\delta\rho}(k)=\frac{2\pi^{2}\rho_{\oplus}^{2}}{m^{4}\sigma^{5}}\frac{\sigma}{v_{k}}\exp\left[-\frac{v_{k}^{2}}{4\sigma^{2}}-\frac{\sigma^{2}}{v_{k}^{2}}\Big{(}\frac{\bm{v}_{0}\cdot\bm{v}_{k}}{\sigma^{2}}-\omega\tau\Big{)}^{2}\right]$
(44)

where $\bm{v}_{k}=\bm{k}/m$ and $\tau=1/m_{\phi}\sigma^{2}$.

##  Appendix B Timing residual

In this Appendix we provide a detailed derivation of the ULDM signals
discussed in Section II.

###  B.1 Time delay

We start by parametrizing the metric perturbations produced by ULDM
fluctuations in terms of the scalar potentials $\Psi$ and $\Phi$:

$ds^{2}=(1+2\Phi)dt^{2}-(1-2\Psi)dx^{2}.$ (45)

Now, consider a bunch of photons emitted by a pulsar at space-time coordinates
$x_{\rm em}=(t_{\rm em},\bm{x}_{p})$ and propagating to an observer on Earth
that receives them at $x_{\rm obs}=(t_{\rm obs},\bm{x}_{e})$. The geodesic
condition, $ds^{2}=0$, then requires that

$t_{\rm obs}=\Delta+t_{\rm em}-\int_{t_{\rm em}}^{t_{\rm
obs}}dt\Big{[}\Phi(t,\bm{x}(t))+\Psi(t,\bm{x}(t))\Big{]}\,,$ (46)

where $\Delta=|\bm{x}_{p}-\bm{x}_{e}|$. Now consider a second bunch of photons
emitted after a full rotation of the pulsar at $x^{\prime}_{\rm
em}=(t^{\prime}_{\rm em},\bm{x}^{\prime}_{p})$ and received at
$x^{\prime}_{\rm obs}=(t^{\prime}_{\rm obs},\bm{x}^{\prime}_{e})$. As before,
we can use the geodesic condition to get

$t_{\rm obs}^{\prime}=\Delta^{\prime}+t_{\rm
em}^{\prime}-\int_{t^{\prime}_{\rm em}}^{t^{\prime}_{\rm
obs}}dt\Big{[}\Phi(t,\bm{x}^{\prime}(t))+\Psi(t,\bm{x}^{\prime}(t))\Big{]}$
(47)

where $\Delta^{\prime}=|\bm{x}^{\prime}_{p}-\bm{x}^{\prime}_{e}|$. Assuming
that the metric perturbations are small (i.e. $\Phi,\,\Psi\ll 1)$, we can
approximate the integral on the right-hand side of this equation as

$\begin{split}&\int_{t^{\prime}_{\rm em}}^{t^{\prime}_{\rm
obs}}dt\Big{[}\Phi(t,\bm{x}^{\prime}(t))+\Psi(t,\bm{x}^{\prime}(t))\Big{]}\\\
\simeq&\int_{t_{\rm em}+T_{a}}^{t_{\rm
obs}+T_{a}}dt\Big{[}\Phi(t,\bm{x}^{\prime}(t))+\Psi(t,\bm{x}^{\prime}(t))\Big{]}\\\
\simeq&\int_{t_{\rm em}}^{t_{\rm
obs}}dt\Big{[}\Phi(t+T_{a},\bm{x}(t))+\Psi(t+T_{a},\bm{x}(t))\Big{]}\end{split}$
(48)

where, in going from the first to the second line, we have used the
approximation $t_{\rm em,obs}^{\prime}\simeq t_{\rm 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\to t+T_{a}$, and used that for the unperturbed photon trajectory,
$\bm{x}(t)=\bm{x}_{e}+(t_{\rm obs}-t)\hat{n}_{a}$, we have that
$\bm{x}^{\prime}(t+T_{a})\simeq\bm{x}(t)$. Therefore, using Eqs. (46)-(48) we
can write

$t^{\prime}_{\rm obs}-t_{\rm obs}\simeq t^{\prime}_{\rm em}-t_{\rm
em}+\Delta^{\prime}-\Delta+\int_{t_{\rm obs}-d_{a}}^{t_{\rm
obs}}dt\Big{[}\Phi(t,\bm{x}(t))-\Phi(t+T_{a},\bm{x}(t))+\Psi(t,\bm{x}(t))-\Psi(t+T_{a},\bm{x}(t))\Big{]}\,,$
(49)

where in rewriting the integration limits we have used that $t_{\rm em}\simeq
t_{\rm 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

$\Delta\tau\simeq\int_{t_{\rm obs}}^{t_{\rm
obs}^{\prime}}dt\Big{[}1+\Phi(t,\bm{x}_{e})\Big{]}\simeq t^{\prime}_{\rm
obs}-t_{\rm obs}+T_{a}\Phi(t_{\rm obs},\bm{x}_{e})\,$ (50)

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.
(49).

The difference between the two emission times, $t^{\prime}_{\rm em}-t_{\rm
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:

$T_{a}\simeq\int_{t_{\rm em}}^{t_{\rm
em}^{\prime}}dt\Big{[}1+\Phi(t,\bm{x}_{p})\Big{]}\simeq t^{\prime}_{\rm
em}-t_{\rm em}+T_{a}\Phi(t_{\rm em},\bm{x}_{p})\,,$ (51)

where to evaluate the integral we have used that the scalar potential is
approximately constant on the timescale of the pulsar period.

The difference $\Delta^{\prime}-\Delta$ 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:

$\frac{\Delta^{\prime}-\Delta}{T_{a}}\simeq\hat{n}_{a}\cdot\int^{t_{\rm
obs}}dt^{\prime}\,\nabla\Big{[}\Phi(t^{\prime},\bm{x}_{e})-\Phi(t^{\prime}-d_{a},\bm{x}_{p})\Big{]}$
(52)

where we have once again used that the scalar potential is approximately
constant over a pulsar rotation and that $t_{\rm em}\simeq t_{\rm obs}-d_{a}$.

The integral in Eq. (49) can be rewritten by expanding at first order the time
dependence of the scalar potential:

$\begin{split}\Psi(t+T_{a},\bm{x})&\simeq\Psi(t,\bm{x})+T_{a}\,\partial_{t}\Psi(t,\bm{x})\\\
\Phi(t+T_{a},\bm{x})&\simeq\Phi(t,\bm{x})+T_{a}\,\partial_{t}\Phi(t,\bm{x})\,.\end{split}$
(53)

By substituting these expansion into Eq. (49), we can rewrite the integral as

$-T_{a}\int_{t_{\rm obs}-d_{a}}^{t_{\rm
obs}}dt\,\partial_{t}\Big{[}\Phi(t,\bm{x}(t))+\Psi(t,\bm{x}(t))\Big{]}\,.$
(54)

then using that $\partial_{t}=d/dt-\hat{n}_{a}\cdot\nabla$, we get

$-T_{a}\Big{[}\Phi(t_{\rm obs},\bm{x}_{e})-\Phi(t_{\rm
obs}-d_{a},\bm{x}_{p})+\Psi(t_{\rm obs},\bm{x}_{e})-\Psi(t_{\rm
obs}-d_{a},\bm{x}_{p})\Big{]}+T_{a}\hat{n}_{a}\cdot\int_{t_{\rm
obs}-d_{a}}^{t_{\rm
obs}}dt\,\nabla\Big{[}\Phi(t,\bm{x}(t))+\Psi(t,\bm{x}(t))\Big{]}$ (55)

Substituting Eqs. (50)–(55) into Eq. (49), we find that for an observer on
Earth, the interval between the arrival time of the two radiation pulses is
given by

$\frac{\Delta\tau}{T_{a}}\simeq 1+\left(\frac{\Delta
T_{a}}{T_{a}}\right)_{1}+\left(\frac{\Delta
T_{a}}{T_{a}}\right)_{2}+\left(\frac{\Delta T_{a}}{T_{a}}\right)_{3}$ (56)

where the three terms contributing to deviations in the observed pulsar’s
rotation period are given by

$\displaystyle\left(\frac{\Delta T_{a}}{T_{a}}\right)_{1}$
$\displaystyle=\Psi(t-d_{a},\bm{x}_{p})-\Psi(t,\bm{x}_{e})$ (57)
$\displaystyle\left(\frac{\Delta T_{a}}{T_{a}}\right)_{2}$
$\displaystyle=\hat{n}_{a}\cdot\int^{t}_{t-d_{a}}dt^{\prime}\,\nabla\Big{[}\Phi(t^{\prime},\bm{x}(t^{\prime}))+\Psi(t^{\prime},\bm{x}(t^{\prime}))\Big{]}$
(58) $\displaystyle\left(\frac{\Delta T_{a}}{T_{a}}\right)_{3}$
$\displaystyle=\hat{n}_{a}\cdot\int^{t}dt^{\prime}\nabla\Big{[}\Phi(t^{\prime},\bm{x}_{e})-\Phi(t^{\prime}-d_{a},\bm{x}_{p})\Big{]}$
(59)

where in all these expressions we have replaced $t_{\rm 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:

$\displaystyle\left(\frac{\Delta T_{a}}{T_{a}}\right)_{1}$
$\displaystyle=-\int\frac{d^{4}k}{(2\pi)^{4}}\,U(k,d_{a},\hat{n}_{a})e^{-ik\cdot
x_{\rm obs}}\tilde{\Psi}(k)$ (60) $\displaystyle\left(\frac{\Delta
T_{a}}{T_{a}}\right)_{2}$
$\displaystyle=-\int\frac{d^{4}k}{(2\pi)^{4}}\,\frac{\hat{n}_{a}\cdot\bm{k}}{\omega+\hat{k}\cdot\hat{n}_{a}}U(k,d_{a},\hat{n}_{a})e^{-ik\cdot
x_{\rm obs}}$
$\displaystyle\times\big{[}\tilde{\Phi}(k)+\tilde{\Psi}(k)\big{]}$ (61)
$\displaystyle\left(\frac{\Delta T_{a}}{T_{a}}\right)_{3}$
$\displaystyle=-\int\frac{d^{4}k}{(2\pi)^{4}}\,\frac{\hat{n}_{a}\cdot\bm{k}}{\omega}U(k,d_{a},\hat{n}_{a})e^{-ik\cdot
x_{\rm obs}}\tilde{\Phi}(k)$ (62)

where $U(k,d_{a},\hat{n}_{a})=(1-e^{i\omega
d_{a}}e^{i\bm{k}\cdot\hat{n}_{a}d_{a}})$. Given that, as we have shown in
Appendix A, for the dominant Fourier modes we have $\omega\lesssim
m_{\phi}\sigma^{2}$ and $|\bm{k}|\lesssim m_{\phi}\sigma$, the Doppler effect,
$(\Delta T_{a}/T_{a})_{3}$, provides the dominant contribution. Note also that
$\tilde{\Phi}\approx\tilde{\Psi}$ for low-frequency stochastic fluctuations.

It is worthwhile to compare our result with Khmelnitsky and Rubakov [5]. In
their work, they have considered a coherently oscillating mode at
$\omega=2m_{\phi}$. From (60)–(61), it is clear that the coherently
oscillating signal is strongest in $(\Delta T_{a}/T)_{1}$, while the other
term is suppressed by the additional power of DM velocity. Using (60), it is
straightforward to reproduce the results of Khmelnitsky and Rubakov as shown
in the Appendix of Ref. [9].

###  B.2 Spectrum

As discussed in the main text, the timing residual power spectrum induced by
stochastic ULDM fluctuations can be written as

$\displaystyle S^{\rm
DM}_{ab}(f)\approx\frac{2G^{2}}{\pi^{2}f^{4}}\\!\int\\!\frac{d^{3}k}{(2\pi)^{3}}(\bm{k}\cdot\hat{n}_{a})(\bm{k}\cdot\hat{n}_{b})U_{a}U_{b}^{*}\frac{P_{\delta\rho}(f,\bm{k})}{k^{4}}\,,$

where the density fluctuation power spectrum, $P_{\delta\rho}$, is given in
(44). For an explicit computation, we set the coordinate as

$\displaystyle\hat{v}_{0}$ $\displaystyle=(0,0,1),$ (63)
$\displaystyle\hat{k}$
$\displaystyle=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),$ (64)
$\displaystyle\hat{n}_{a}$
$\displaystyle=(\sin\theta_{a}\cos\phi_{a},\sin\theta_{a}\sin\phi_{a},\cos\theta_{a}),$
(65) $\displaystyle\hat{n}_{b}$
$\displaystyle=(\sin\theta_{b}\cos\phi_{b},\sin\theta_{b}\sin\phi_{b},\cos\theta_{b}).$
(66)

After performing an explicit computation, one finds that the timing residual
power spectrum can be written as a sum of two terms:

$\displaystyle
S^{\scriptscriptstyle\textrm{DM}}(f)\Gamma_{ab}^{\scriptscriptstyle\textrm{DM}}(f)=\Gamma^{\perp}_{ab}S^{\perp}(f)+\Gamma^{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}_{ab}S^{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}(f),$
(67)

where the overlap functions are given in (18)–(19). Each spectrum is given as

$\displaystyle S^{A}(f)$ $\displaystyle=\frac{\bar{a}^{2}\tau}{(2\pi
f)^{4}}\left[\frac{64}{\pi}\int\frac{dx}{x}e^{-x^{2}/4}C_{A}(V_{0},\bar{\omega},x)\right]$
(68) $\displaystyle\equiv\frac{\bar{a}^{2}\tau}{(2\pi
f)^{4}}H^{A}(f,V_{0})\,,$ (69)

where we have introduced $V_{0}=v_{0}/\sigma$ and $\bar{\omega}=\omega\tau$
and changed the integration variable to $x=v_{k}/\sigma$. Here the function
$C_{\perp}$ and
$C_{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}$
are given as

$\displaystyle C_{\perp}(V_{0},\bar{\omega},x)$
$\displaystyle=+\frac{1}{8V_{0}^{3}}\Big{[}\big{(}V_{0}-\frac{\bar{\omega}}{x}\big{)}e^{-(V_{0}+\frac{\bar{\omega}}{x})^{2}}-\frac{\sqrt{\pi}}{2}\Big{(}1+2\Big{(}\frac{\bar{\omega}}{x}\Big{)}^{2}-2V_{0}^{2}\Big{)}{\rm
erf}(V_{0}-\frac{\bar{\omega}}{x})\Big{]}+(\bar{\omega}\to-\bar{\omega}),$
(70) $\displaystyle
C_{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}(V_{0},\bar{\omega},x)$
$\displaystyle=-\frac{1}{4V_{0}^{3}}\Big{[}\big{(}V_{0}-\frac{\bar{\omega}}{x}\big{)}e^{-(V_{0}+\frac{\bar{\omega}}{x})^{2}}-\frac{\sqrt{\pi}}{2}\Big{(}1+2\Big{(}\frac{\bar{\omega}}{x}\Big{)}^{2}\Big{)}{\rm
erf}(V_{0}-\frac{\bar{\omega}}{x})\Big{]}+(\bar{\omega}\to-\bar{\omega}),$
(71)

In the limit $V_{0}\to 0$, one finds
$H^{\perp}(f)=H^{\mathrel{\mathchoice{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\displaystyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{4.78333pt}{\raisebox{0.0pt}{$\textstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{3.34833pt}{\raisebox{0.0pt}{$\scriptstyle/\mkern-5.0mu/$}}}{\leavevmode\resizebox{}{2.39166pt}{\raisebox{0.0pt}{$\scriptscriptstyle/\mkern-5.0mu/$}}}}}(f)=(64/3\pi)K_{0}(\omega\tau)$,
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.

Figure 4: The 2D posterior distribution for the analysis of NANOGrav 12.5-year
data is shown. For a detailed description of the analysis, see Section III.
The quoted numbers above each panel correspond to the median and 90% credible
interval of each parameter, while each contour in the 2D histogram contains
$68.3,\,95.4$, and $99.7$% of samples.

##  Appendix C Posterior distribution

In Figure 4, we present the 2D posterior distribution for 5 parameters,

$(\log_{10}m_{\phi},\,\log_{10}(\rho/\rho_{0}),\,\sigma,\,\gamma_{\rm
GW},\,\log_{10}A_{\rm GW}),$

obtained from the analysis of NANOGrav 12.5-year data. For detailed
descriptions of the analysis, see Section III.

## References

  * Agazie  _et al._ [2023a] G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L8 (2023a), arXiv:2306.16213 [astro-ph.HE] .
  * Antoniadis  _et al._ [2023] J. Antoniadis _et al._ , (2023), arXiv:2306.16214 [astro-ph.HE] .
  * Reardon  _et al._ [2023] D. J. Reardon _et al._ , Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Xu  _et al._ [2023] H. Xu _et al._ , Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Khmelnitsky and Rubakov [2014] A. Khmelnitsky and V. Rubakov, JCAP 02, 019 (2014), arXiv:1309.5888 [astro-ph.CO] .
  * Graham  _et al._ [2016] P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, and W. A. Terrano, Phys. Rev. D 93, 075029 (2016), arXiv:1512.06165 [hep-ph] .
  * Kaplan  _et al._ [2022] D. E. Kaplan, A. Mitridate, and T. Trickle, Phys. Rev. D 106, 035032 (2022), arXiv:2205.06817 [hep-ph] .
  * Afzal  _et al._ [2023] A. Afzal _et al._ (NANOGrav), Astrophys. J. Lett. 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Kim [2023] H. Kim, JCAP 12, 018 (2023), arXiv:2306.13348 [hep-ph] .
  * Schive  _et al._ [2014] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Hui  _et al._ [2017] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Bar-Or  _et al._ [2019] B. Bar-Or, J.-B. Fouvry, and S. Tremaine, Astrophys. J. 871, 28 (2019), arXiv:1809.07673 [astro-ph.GA] .
  * Read [2014] J. I. Read, J. Phys. G 41, 063101 (2014), arXiv:1404.1938 [astro-ph.GA] .
  * de Salas and Widmark [2021] P. F. de Salas and A. Widmark, Rept. Prog. Phys. 84, 104901 (2021), arXiv:2012.11477 [astro-ph.GA] .
  * Alam  _et al._ [2021] M. F. Alam _et al._ (NANOGrav), Astrophys. J. Suppl. 252, 4 (2021), arXiv:2005.06490 [astro-ph.HE] .
  * Pitjev and Pitjeva [2013] N. P. Pitjev and E. V. Pitjeva, Astron. Lett. 39, 141 (2013), arXiv:1306.5534 [astro-ph.EP] .
  * Kim and Lenoci [2022] H. Kim and A. Lenoci, Phys. Rev. D 105, 063032 (2022), arXiv:2112.05718 [hep-ph] .
  * Kim  _et al._ [2023] H. Kim, A. Lenoci, G. Perez, and W. Ratzinger, (2023), arXiv:2307.14962 [hep-ph] .
  * Hazboun  _et al._ [2019] J. S. Hazboun, J. D. Romano, and T. L. Smith, Phys. Rev. D 100, 104028 (2019), arXiv:1907.04341 [gr-qc] .
  * Binney and Merrifield [1998] J. Binney and M. Merrifield,  _Galactic Astronomy_ (1998).
  * Schönrich  _et al._ [2010] R. Schönrich, J. Binney, and W. Dehnen, Monthly Notices of the Royal Astronomical Society 403, 1829 (2010).
  * Arzoumanian  _et al._ [2016] Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, _et al._ (NANOGrav), Astrophys. J. 821, 13 (2016), arXiv:1508.03024 [astro-ph.GA] .
  * Arzoumanian  _et al._ [2018] Z. Arzoumanian, P. T. Baker, A. Brazier, _et al._ (NANOGRAV), Astrophys. J. 859, 47 (2018), arXiv:1801.02617 [astro-ph.HE] .
  * Hellings and Downs [1983] R. w. Hellings and G. s. Downs, Astrophys. J. Lett. 265, L39 (1983).
  * Johnson  _et al._ [2023] A. D. Johnson _et al._ (NANOGrav), (2023), arXiv:2306.16223 [astro-ph.HE] .
  * Ellis  _et al._ [2019] J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, “ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE,” Astrophysics Source Code Library, record ascl:1912.015 (2019), ascl:1912.015 .
  * Ellis  _et al._ [2020] J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, “Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite,” Zenodo (2020).
  * Chalumeau  _et al._ [2021] A. Chalumeau _et al._ (EPTA), Mon. Not. Roy. Astron. Soc. 509, 5538 (2021), arXiv:2111.05186 [astro-ph.HE] .
  * Arzoumanian  _et al._ [2020] Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Goncharov  _et al._ [2021] B. Goncharov _et al._ , Astrophys. J. Lett. 917, L19 (2021), arXiv:2107.12112 [astro-ph.HE] .
  * Perera  _et al._ [2019] B. B. P. Perera _et al._ , Mon. Not. Roy. Astron. Soc. 490, 4666 (2019), arXiv:1909.04534 [astro-ph.HE] .
  * Ellis and van Haasteren [2017] J. Ellis and R. van Haasteren, “jellis18/ptmcmcsampler: Official release,”  (2017).
  * Vallisneri [2020] M. Vallisneri, “libstempo: Python wrapper for Tempo2,” Astrophysics Source Code Library, record ascl:2002.017 (2020), ascl:2002.017 .
  * Edwards  _et al._ [2006] R. T. Edwards, G. B. Hobbs, and R. N. Manchester, Mon. Not. Roy. Astron. Soc. 372, 1549 (2006), arXiv:astro-ph/0607664 .
  * Hobbs  _et al._ [2006] G. Hobbs, R. Edwards, and R. Manchester, Mon. Not. Roy. Astron. Soc. 369, 655 (2006), arXiv:astro-ph/0603381 .
  * Agazie  _et al._ [2023b] G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L9 (2023b), arXiv:2306.16217 [astro-ph.HE] .
  * Adler [2008] S. L. Adler, J. Phys. A 41, 412002 (2008), arXiv:0808.0899 [astro-ph] .
  * Tsai  _et al._ [2022] Y.-D. Tsai, J. Eby, J. Arakawa, D. Farnocchia, and M. S. Safronova, (2022), arXiv:2210.03749 [hep-ph] .
  * Read  _et al._ [2008] J. I. Read, G. Lake, O. Agertz, and V. P. Debattista, Mon. Not. Roy. Astron. Soc. 389, 1041 (2008), arXiv:0803.2714 [astro-ph] .
  * Evans  _et al._ [2019] N. W. Evans, C. A. J. O’Hare, and C. McCabe, Phys. Rev. D 99, 023012 (2019), arXiv:1810.11468 [astro-ph.GA] .
  * O’Hare  _et al._ [2018] C. A. J. O’Hare, C. McCabe, N. W. Evans, G. Myeong, and V. Belokurov, Phys. Rev. D 98, 103006 (2018), arXiv:1807.09004 [astro-ph.CO] .
  * O’Hare  _et al._ [2020] C. A. J. O’Hare, N. W. Evans, C. McCabe, G. Myeong, and V. Belokurov, Phys. Rev. D 101, 023006 (2020), arXiv:1909.04684 [astro-ph.GA] .
  * Banerjee  _et al._ [2020a] A. Banerjee, D. Budker, J. Eby, H. Kim, and G. Perez, Commun. Phys. 3, 1 (2020a), arXiv:1902.08212 [hep-ph] .
  * Banerjee  _et al._ [2020b] A. Banerjee, D. Budker, J. Eby, V. V. Flambaum, H. Kim, O. Matsedonskyi, and G. Perez, JHEP 09, 004 (2020b), arXiv:1912.04295 [hep-ph] .
  * Budker  _et al._ [2023] D. Budker, J. Eby, M. Gorghetto, M. Jiang, and G. Perez, (2023), arXiv:2306.12477 [hep-ph] .
  * Derevianko [2018] A. Derevianko, Phys. Rev. A 97, 042506 (2018), arXiv:1605.09717 [physics.atom-ph] .
  * Foster  _et al._ [2018] J. W. Foster, N. L. Rodd, and B. R. Safdi, Phys. Rev. D 97, 123006 (2018), arXiv:1711.10489 [astro-ph.CO] .


$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)\geq
1-c\sqrt{\frac{N^{2}L^{2}\log{NL}}{n}}-c\sqrt{\frac{t}{n}}.$

###  A.11 Proof of Corollary 4.4

To establish Corollary 4.4 we need to find $\mathcal{R}_{n}(\mathcal{F})$ and
$\mathcal{R}_{n}(\mathcal{G})$ in light of Theorem 3.7. Recall that here
$\Phi(\cdot)$ is the feature map, which maps $x\mapsto K(x,\cdot)$. As we have
assumed $\sup_{x}K(x,x)\leq b$, we have for any $h$ in the RKHS generated by
$K$:

$|h(x)|=|\langle
h,K(x,\cdot)\rangle_{\mathcal{H}}|\leq\|h\|_{H}\|K(x,\cdot)\|_{\mathcal{H}}\leq\sqrt{b}\|h\|_{\mathcal{H}}\,.$
(A.8)

In other words, $\|h\|_{\infty}\leq\sqrt{b}\|h\|_{\mathcal{H}}$. Furthermore,
for any $f\in\mathcal{F}$:

$\displaystyle|f(x)|=\left|\left\langle\Phi(x),\mathcal{A}(\Phi)(x)\right\rangle_{\mathcal{H}}\right|$
$\displaystyle=\left|\left\langle\mathcal{A},\Phi(x)\otimes\Phi(x)^{\top}\right\rangle\right|$
$\displaystyle\leq\|\mathcal{A}\|_{\star}\|\Phi(x)\otimes\Phi(x)^{\top}\|_{\mathrm{op}}$
$\displaystyle\leq r\sup_{h\in\mathcal{H}:\|h\|_{\mathcal{H}}\leq
1}\left|\left\langle
h,\Phi(x)\otimes\Phi(x)^{\top}(h)\right\rangle_{\mathcal{H}}\right|$
$\displaystyle\leq r\sup_{h\in\mathcal{H}:\|h\|_{\mathcal{H}}\leq
1}|h^{2}(x)|\leq rb\,.$ (A.9)

The first inequality follows from duality of operator and nuclear norm, second
inequality follows from the definition of operator norm and the last
inequality is from (A.8). Therefore, $\|f\|_{\infty}\leq rb$ for any
$f\in\mathcal{F}$. We now use these bounds to control the Rademacher
complexity:

$\displaystyle\mathcal{R}_{n}(\mathcal{F})$
$\displaystyle=\mathbb{E}\left[\sup_{\|\mathcal{A}\|_{\star}\leq
r}\left|\frac{1}{n}\sum_{i}{\epsilon}_{i}\left\langle\Phi(X_{i}),\mathcal{A}(\Phi)(X_{i})\right\rangle\right|\right]$
$\displaystyle\leq
r\mathbb{E}\left[\left\|\frac{1}{n}\sum_{i}{\epsilon}_{i}\Phi(X_{i})\otimes\Phi(X_{i})^{\top}\right\|_{\mathrm{op}}\right]$
$\displaystyle=r\mathbb{E}\left[\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\frac{1}{n}\sum_{i}{\epsilon}_{i}\langle\Phi(X_{i}),h\rangle^{2}\right|\right]$
$\displaystyle=r\mathbb{E}\left[\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\frac{1}{n}\sum_{i}{\epsilon}_{i}h^{2}(X_{i})\right|\right]$

Next we apply Leduox-Talagrand contraction inequality (see Theorem 7 of Duchi,
(2009)) as $\|h\|_{\infty}\leq\sqrt{b}$ by (A.8) and $\phi(x)=x^{2}$ is
Lipschitz function on $[-\sqrt{b},\sqrt{b}]$ with Lipschitz co-efficient
$2\sqrt{b}$. Therefore we have:

$\displaystyle\mathcal{R}_{n}(\mathcal{F})$ $\displaystyle\leq 4r\sqrt{b}\
\mathbb{E}\left[\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\frac{1}{n}\sum_{i}{\epsilon}_{i}h(X_{i})\right|\right]$
$\displaystyle=4r\sqrt{b}\ \mathbb{E}\left[\sup_{\|h\|_{\mathcal{H}}\leq
1}\left|\langle
h,\frac{1}{n}\sum_{i}{\epsilon}_{i}K(X_{i},\cdot)\rangle\right|\right]$
$\displaystyle\leq 4r\sqrt{b}\
\mathbb{E}\left[\left\|\frac{1}{n}\sum_{i}{\epsilon}_{i}K(X_{i},\cdot)\right\|_{\mathcal{H}}\right]$
$\displaystyle\leq 4r\sqrt{b}\
\sqrt{\mathbb{E}\left[\left\|\frac{1}{n}\sum_{i}{\epsilon}_{i}K(X_{i},\cdot)\right\|^{2}_{\mathcal{H}}\right]}$
$\displaystyle=4r\sqrt{b}\
\sqrt{\mathbb{E}\left[\frac{1}{n^{2}}\sum_{i,j}{\epsilon}_{i}{\epsilon}_{j}K(X_{i},X_{j})\right]}$
$\displaystyle=\frac{4r\sqrt{b}}{n}\ \sqrt{\mathbb{E}\left[{\sf
tr}{K}\right]}=\frac{4r\sqrt{b}}{\sqrt{n}}\sqrt{\mathbb{E}\left[K(X,X)\right]}\,.$

Combining the above inequality with Theorem 3.7, we get the bound of the
obtained bandwidth. The coverage guarantee directly follows from the
observation that
$\mathcal{R}_{n}(\mathcal{G})\lesssim\mathcal{R}_{n}(\mathcal{M})+\mathcal{R}_{n}(\mathcal{F})$.

##  Appendix B Simulations

###  B.1 Simulation results for multivariate input

For the multivariate setups, we generate the inputs $X_{i}$ uniformly from a
unit sphere and define $\beta=[1/\sqrt{3},1/\sqrt{3},-1/\sqrt{3}]$. The
outputs are then generated from the following procedure:

  1. 1.

Setup 1: $Y_{i}=\sqrt{v_{0}(X^{T}_{i}\beta)}\xi_{i}$ where
$\xi_{i}\perp\\!\\!\\!\perp X_{i}$, $\xi_{i}\sim{\rm Unif}([-1,1])$ and
$v_{0}(x)=1+25x^{4}$.

  2. 2.

Setup 2: $Y_{i}=m_{0}(X^{T}_{i}\beta)+\sqrt{v_{0}(X^{T}_{i}\beta)}\xi_{i}$,
where $\xi_{i}\perp\\!\\!\\!\perp X_{i}$ and is generated from ${\rm
Unif}([-1,1])$. The true mean function $m_{0}(x)=1+5x^{3}$ and
$v_{0}(x)=1+25x^{4}$.

  3. 3.

Setup 3: $Y_{i}$ is generated from a conditional Laplace distribution centered
around its mean $(X^{T}_{i}\beta)^{2}+5(X^{T}_{i}\beta)^{4}$, where $Y$ is
conditioned to lie within two times standard deviation
$\sqrt{1+25(X^{T}_{i}\beta)^{4}}$.

Since we are working under the multivariate scenarios, we slightly increase
the number of samples, i.e., we generate $n=3000$ i.i.d samples from the above
data generating process with $n_{pre}=1000$, $n_{opt}=900$, $n_{adj}=100$ and
$n_{t}=1000$. We then employ the same estimation process outlined in Section
5. Our results are shown in the following table.

Table 6: Coverage and Average Width for Multivariate Input | Setup 1 | Setup 2 | Setup 3
---|---|---|---
| Coverage | Width | Coverage | Width | Coverage | Width
UTOPIA | $96.4\%$ | $6.5557$ | $98.7\%$ | $8.6933$ | $93.6\%$ | $17.3982$
LQR | $95.0\%$ | $10.7780$ | $95.5\%$ | $11.2241$ | $92.7\%$ | $28.6803$
SplitCF | $94.5\%$ | $10.0948$ | $96.4\%$ | $11.4834$ | $91.9\%$ | $21.4286$

Here we have not included a comparison with the SDP method proposed by Liang,
(2021). This is primarily due to its considerable computational demands and
substantial memory requirements.