Probing ultralight scalar, vector and tensor dark matter with pulsar timing arrays

Abstract

Pulsar timing arrays (PTAs) are sensitive to oscillations in the gravitational potential along the line-of-sight due to ultralight particle pressure. We calculate the probing power of PTAs for ultralight bosons across all frequencies, from those larger than the inverse observation time to those smaller than the inverse distance to the pulsar. We show that since the signal amplitude grows comparably to the degradation in PTA sensitivity at frequencies smaller than inverse observation time, the discovery potential can be extended towards lower masses by over three decades, maintaining high precision. We demonstrate that, in the mass range $10^{-26}-10^{-23}$ eV, existing 15-year PTA data can robustly detect or rule out an ultralight component down to $O(1-10)\%$ of the total dark matter. Non-detection, together with other bounds in different mass ranges, will imply that ultralight scalar/axion can comprise at most $1-10\%$ of dark matter in the $10^{-30}\!-\!10^{-17}$ eV range. With 30 years of observation, current PTAs can extend the reach down to $0.1-1\%$, while next-generation PTAs such as SKA can attain the $0.01-0.1\%$ precision. We generalize and derive predictions for ultralight spin-1 vector (i.e. dark photon) and spin-2 tensor dark components.

Ultralight particles are intriguing dark matter candidates, motivated by high energy theories Svrcek:2006yi . As they have de-Broglie wavelengths comparably to galactic or even larger cluster scales, depending on their mass, they would leave marks on cosmological and astrophysical observables Arvanitaki:2009fg . One such effect comes from the pressure of the ultralight field which could cause the energy-momentum tensor and equivalently spacetime to fluctuate in a monochromatic way Khmelnitsky:2013lxt . These fluctuations would modify the observed arrival time of pulses from pulsars and can be used to extract the properties, e.g. mass and energy density, of ultralight scalar/axion-like (spin-0), vector (dark photon/spin-1) and tensor (spin-2) bosons, even if they form a fraction Blum:2021oxj ; Ye:2021iwa of dark matter.

Ultralight bosonic degrees of freedom have undergone intense exploration in recent years (see Ref. Ferreira:2020fam for a review). The mass parameter space has already been constrained at several distinct scales by a number of probes. Cosmic microwave background experiments (CMB) Hlozek:2014lca ; Poulin:2018dzj combined with large-scale structure surveys Lague:2021frh have derived bounds on the extremely light range, roughly $10^{-32}\!-\!10^{-25}\,{\rm eV}$. The ultraviolet luminosity function and reionization constrain the $10^{-23}\!-\!10^{-22}\,{\rm eV}$ range Bozek:2014uqa . Lyman-$\alpha$ experiments reaching large wavenumbers of the matter power spectrum probe approximately the range $10^{-23}\!-\!10^{-20.5}\,{\rm eV}$ Kobayashi:2017jcf ; Irsic:2017yje ; Armengaud:2017nkf ; Rogers:2020ltq ; Rogers:2020cup , and the shape of Eridanus-II the $10^{-20.5}-10^{-19}\,{\rm eV}$ rangeMarsh:2018zyw . Galaxy rotation curves probe the $10^{-23.7}-10^{-22}\,{\rm eV}$ range Bar:2021kti . Secular variations in binary-pulsar orbital parameters can test the mass range $10^{-23}\!-\!10^{-18}\,{\rm eV}$ Blas:2019hxz ; Armaleo:2019gil ; LopezNacir:2018epg ; Blas:2016ddr . Finally, non-observation of superradiance in rapidly spinning supermassive black holes (SMBHs) can explore the $10^{-20.5}\!-\!10^{-16.7}\,{\rm eV}$ range Unal:2020jiy . With all of the above, the mass range $10^{-26}\!-\!10^{-23}\,{\rm eV}$ still remains relatively unexplored.

In this Letter , we show that current data from the International Pulsar-Timing Array (IPTA) Manchester:2013ndt ; Verbiest:2016vem , a collaboration between NANOGRAV NANOGrav:2015aud , PPTA Manchester:2012za and EPTA Lentati:2015qwp can be used to probe the ultralight boson energy density fraction ${\cal F}_{\rm ULDM}$ in the $10^{-26}\!-\!10^{-23}\,{\rm eV}$ mass range with about 1-10% precision. The bounds remarkably extend to the smaller masses (lower frequencies) in this window, since the strength of the ultralight particle signal changes comparable to the PTA sensitivity curve. Therefore, smaller mass particles can be efficiently constrained by current IPTA data down to frequencies corresponding to $1/D_{\rm pulsar}\!\sim\!1/{\rm kpc}$ (where $D_{\rm pulsar}$ is the distance to the pulsar, and we set the speed of light $c=1$). Extending the observation period up to 30 years can improve these constraints by one order of magnitude, and with future PTA experiments, such as SKA Weltman:2018zrl , by two more orders of magnitude. Combined with CMB/LSS, Lyman-$\alpha$ and SMBH superradiance constraints, the discovery potential of PTAs implies that ultralight bosons can be probed continuously throughout the mass range $10^{-30}\!-\!10^{-17}$ eV with $\sim{\cal O}(1)\%$ precision (see Hlozek:2016lzm ; Safarzadeh:2019sre ; Dentler:2021zij ; Flitter:2022pzf ; Libanore:2022ntl ; Sun:2021yra ).

We start with the signature of ultralight dark matter in PTAs. Free scalar/axion-like degrees of freedom, $\phi$, with a canonical kinetic term have a Lagrangian density ${\cal L}=\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{1}{2}m^{2}\phi^{2}$. If such particles are non-relativistic, then the field configuration is given by a plane wave of nearly single frequency with corrections up to the kinetic term, which is about $10^{-6}$ times smaller compared to the rest mass. If we neglect the expansion of the background, $\phi\!=\!{\cal A}(x)\,\cos(mt+\beta)$. The energy density is nearly time independent, $\rho=\dot{\phi}^{2}/2+V\simeq\frac{1}{2}m^{2}{\cal A}^{2}$, and the pressure oscillates with angular frequency twice the mass of the particle, i.e. $p=\dot{\phi}^{2}/2-V\simeq\rho\cdot\cos(2\pi ft+2\beta)$, where $f$ is the oscillation frequency that can be computed using the corresponding mass as $2\pi f=w=2m$, and $\beta$ is a phase. Hence,

$$f=5\cdot 10^{-9}{\rm Hz}\left(\frac{m}{10^{-23}\,{\rm eV}}\right).$$ (1)

Then, solving the 00 and ii components of Einstein’s equations at first order, with the following metric perturbations: $\delta g_{00}=-2\Phi$ and $\delta g_{ij}=2\psi\delta_{ij}$, one ends up with $\psi=\psi_{0}(x)+\frac{\pi}{2}G_{N}{\cal A}^{2}\cos(2mt+2\beta)=\psi_{0}(x)+\pi G_{N}\rho/m^{2}\cos(2mt+2\beta)$, the second term being the oscillating piece defined as $\psi_{c}\!=\!\pi G_{N}\rho/m^{2}$. The time residual can be calculated via the fractional frequency shift as

$$\delta\Delta t=\int_{t_{p}}^{t}\frac{\nu^{\prime}-{\bar{\nu}}}{{\bar{\nu}}}dt,$$ (2)

where $t$ and $t_{p}$ are the time at Earth and pulsar emission, and the fractional change in the frequency is given by

$$\frac{\nu^{\prime}-{\bar{\nu}}}{{\bar{\nu}}}=\psi(x,t)-\psi(x_{p},t_{p})-\int_{t_{p}}^{t}n_{i}\partial_{i}(\Phi+\psi)dt^{\prime},$$ (3)

where the first term is the difference in the gravitational potential and the second term is the change of its gradient during propagation, which is suppressed by a factor $k/m$. Now, plugging Eq. (3) into Eq. (2), we get the expression

$$\delta\Delta t = \frac{\psi_{c}}{m}\sin(mD_{\rm pulsar}+\beta_{e}-\beta_{p})\times$$ (4)

$$\times\cos(2mt-mD_{\rm pulsar}+\beta_{e}+\beta_{p}).$$

where $\beta_{e}$ and $\beta_{p}$ corresponds to phase values at the Earth and at the individual pulsar. The typical amplitude for the residual oscillating signal is the root-mean-square (rms) path average which is given by

$$\sqrt{\langle(\delta\Delta t)^{2}\rangle}=\sqrt{\frac{1}{L}\int_{0}^{L}dl\,(\delta\Delta t)^{2}}={\cal P}\cdot\psi_{c}/m,$$ (5)

where $L\equiv D_{\rm pulsar}$ and ${\cal P}$ is defined as

$$\!\!\!\!{\cal P}=\frac{1}{\sqrt{2}}\left(1-\frac{\sin(2\,m\,L)}{2\,m\,L}\right)^{\frac{1}{2}}\!\!\,\,$$ (6)

In the case where $m\cdot D_{\rm pulsar}>1$ we have ${\cal P}\simeq 1$. However, in the limit $m\cdot D_{\rm pulsar}\ll 1$, we have extra suppression factor in the signal, of approximately $m\,D_{\rm pulsar}$.

The oscillations in the energy-momentum tensor will cause oscillations in the gravitational potential which correspond to an equivalent characteristic strain $h_{c}=2\sqrt{3}\psi_{c}$, stressing that this is not GW signal. A similar analysis can be done for vector and tensor fields. For the massive spin-1 case, i.e. a massive dark photon, the field has 2 transverse and 1 longitudinal polarizations. For the spin-2 case, i.e. tensor, the field has 1 scalar, 2 transverse and 2 transverse-traceless polarizations (for subluminal gravitational waves Qin:2020hfy ). These polarizations are of similar order-of-magnitude in strength. Although there will be accompanying polarizations, the distinctive modes will be transverse polarization for a vector field and transverse-traceless polarization for a tensor field, and so we will conduct our analysis for the detectability of these specific polarizations. The signal from a single pulsar corresponding to the scalar Khmelnitsky:2013lxt , vector Nomura:2019cvc and tensor Armaleo:2020yml dark matter scenarios can be estimated as

$$\!\!\!\!\!\!\!\!h_{c,\,{\rm scalar}}\simeq 2\cdot 10^{-15}\cdot{\cal P}\cdot{\cal F}_{\rm ULDM}\left(\frac{5{\rm nHz}}{f}\right)^{2}\,$$

$$\!\!\!\!\!\!\!\!h_{c,\,{\rm vector}}\simeq 6\cdot 10^{-15}\cdot{\cal P}\cdot{\cal F}_{\rm ULDM}\left(\frac{5{\rm nHz}}{f}\right)^{2}\,$$

$$\!\!\!\!\!\!\!\!h_{c,{\rm tensor},\,\alpha}\simeq 5\cdot 10^{-15}\cdot{\cal P}\cdot{\cal F}^{1/2}_{\rm ULDM}\,\,\left(\frac{5{\rm nHz}}{f}\right)\left(\frac{\alpha}{10^{-7}}\right)\;\;\;$$ (7)

where ${\cal F}_{\rm ULDM}\equiv(\frac{\rho_{\rm ULDM}}{0.4\,GeV/cm^{3}})$ is the fraction of ultralight dark matter (ULDM) in the local Universe and $\alpha$ is the spin-2 universal coupling to matter.

A number of comments are in order. First, the signal suppression from propagation is important, since with the corresponding modification the signal strain grows more slowly than the red-noise strain and as a result the signal-to-noise diminishes for frequencies smaller than inverse pulsar distance, i.e. $f<1/D_{\rm pulsar}$, see Fig. 1.

Secondly, the tensor case is special in the sense that the only known viable massive spin-2 theory, bigravity, includes a universal direct coupling term to matter with strength $\alpha$ that cannot be tuned away. Because of this, the signal from the massive tensor field is proportional to $f^{-1}$ instead of $f^{-2}$ as in the scalar and vector cases Armaleo:2020yml . Therefore, in the tensor case the constraint is effectively on the combination $\alpha^{2}{\cal F}_{\rm ULDM}$ ${}^{,}\!\!$ .

Figure 1: Characteristic strain, $h_{c}$, for the signal in the case of scalar ultralight particles with ${\cal F}_{\rm ULDM}=1$, and for the noise with current PTAs (assuming 15 and 30 years of observation time Chen:2021rqp , for 60 and 90 pulsars) and with SKA (15 years with 5K pulsars), respectively. The effective noise curve is defined here as $h_{\rm eff}\equiv h_{\rm noise}(2{\cal R}_{\beta}^{2})^{-1/4}$. The signal, the PTA sensitivity and the detectability in three regimes separated by vertical dashed lines are discussed in detail in the text. Note that the signal grows comparable to the the noise in the regime $1/D_{\rm pulsar}<f<1/T_{\rm obs}$, allowing for precision probing, and slower for $f>1/T_{\rm obs}$ and $f<1/D_{\rm pulsar}$.

Having described the three different types of signals we are after, we now turn to a discussion of the sensitivity of current and future PTAs. Although all ultralight particles generate pressure which leads to perturbations in the arrival time of pulsar signals, the type of perturbations are distinct from each other. Accordingly, PTAs have different response functions for scalar, vector and tensor fields. The variables typically used in the literature to indicate the strength of the signal: the characteristic strain $h_{c}$, the power spectral density $S_{h}$ and the energy density $\Omega_{h}$, are connected to each other as follows

$$H_{0}^{2}\,\Omega(f)=\frac{2\pi^{2}}{3}f^{3}S_{h}(f)=\frac{2\pi^{2}}{3}f^{2}h_{c}^{2}(f).$$ (10)

To estimate the detectability of the monochromatic signal, we start from the formula

$$\!\!\!\!\!\!\!{\rm SNR}^{2}=2\;(f\cdot T_{\rm obs})^{2}\sum_{I=1}^{N}\sum_{J>I}^{N}\;r_{(\beta)\;IJ}^{2}\;\left(\frac{S_{h,{\rm signal}}}{S_{h,{\rm noise}}}\right)^{2},\,\,\,\,\,$$ (11)

where $T_{\rm obs}$ is the observation time, and $I,J$ indicate two pulsars in each pair correlation, and $r_{(\beta){IJ}}^{2}\equiv\frac{1}{4\pi}\int d{\Omega}\;\chi^{2}_{IJ}$, $\beta$=Scalar (S), Vector(V), Tensor(T), and $N$ is the total number of pulsars, assumed to have identical properties to derive an ensemble average. The explicit expressions for $\chi_{IJ}$, the correlation coefficient between each pair of pulsars separated by an angle $\zeta$, for each $\beta$ are Jenet:2014bea

$${\rm S}:\;\chi_{IJ}=1,~{}~{}~{}~{}~{}~{}~{}{\rm V}:\chi_{IJ}=\frac{1}{3}\cos{\zeta}$$

$${\rm T}:\chi_{IJ}=\frac{1}{2}-\frac{1}{4}\left(\frac{1-\cos{\zeta}}{2}\right)\left(1+6\ln{\left(\frac{1-\cos{\zeta}}{2}\right)}\right)\,$$

Defining ${\cal R}^{2}_{\beta}\equiv\sum_{I=1}^{N}\sum_{J>I}^{N}r^{2}_{(\beta)IJ}$, this yields

$$\!\!\!\!\!\!{\cal R}_{\rm S}^{2}=\frac{N\,(N-1)}{2},\;{\cal R}_{\rm V}^{2}=\frac{N\,(N-1)}{2\cdot 27},\;{\cal R}_{\rm T}^{2}=\frac{N\,(N-1)}{2\cdot 48}.$$ (12)

The PTA sensitivity curve can be described as

$$S_{n}=12\pi^{2}f^{2}\,\frac{{\cal N}}{{\cal T}}=h_{n}^{2}/f$$ (13)

where ${\cal T}=\left(\frac{(f\cdot T_{\rm obs})^{3}}{1+(f\cdot T_{\rm obs})^{3}}\right)^{2}$ is the transmission function, ${\cal N}=2\sigma^{2}\Delta t+P_{r}\,f^{-\gamma}$ (white noise plus red noise which can be relevant for small frequencies). $S_{n}$ varies between the different frequency ranges (see Fig. 1):

Figure 2: Bounds on the fraction of a scalar ultralight boson component out of the total dark matter. We show existing constraints from CMB Hlozek:2014lca ; Poulin:2018dzj ; combined with BOSS Lague:2021frh ; reionization Bozek:2014uqa ; Lyman-$\alpha$ Kobayashi:2017jcf ; Irsic:2017yje ; Armengaud:2017nkf ; Rogers:2020ltq ; Rogers:2020cup ; Eridanus-II Marsh:2018zyw ; galaxy rotation curves Bar:2021kti ; and SMBH superradiance Unal:2020jiy . Our forecasts for current PTA (future; 30 year PTA, and 15 year SKA) are shown in solid (dashed) blue.

• I) $f>1/T_{\rm obs}$ : The high frequency part of the sensitivity regime is controlled by white timing noise, whose expression is given as $S_{n}(f)=12\pi^{2}f^{2}(2\Delta t\sigma^{2})=h_{n}^{2}/f$. Here $\Delta t$ is the timing period, and $\sigma$ is the rms error in timing residuals. This frequency regime has scaling $h_{n}\propto f^{3/2}$, or equivalently $\Omega_{n}\propto f^{5}$ Thrane:2013oya ; Moore:2014lga ; Romano:2016dpx ; Vigeland:2016nmm .

• II) $f<1/T_{\rm obs}$: Here $S_{n}$ is set by three factors: (i) the transmission function, which accounts for the information absorbed by the timing model fit (we assume a quadratic spin-down model which fits for the pulsar phase offset, spin period and period derivative). In this regime, the transmission function scales as $f^{-6}$ and limits the detection capability  Hazboun:2019vhv . (ii) Measurement white noise, which is frequency-independent; (iii) Pulsar specific red noise Goncharov:2021oub , which is more effective at low frequencies Lam:2016iie (see also NANOGrav:2020bcs ; Goncharov:2020krd ). For simplicity, in making our forecast we will limit ourselves to a subset of pulsars for which, in the frequency regime we focus on, the red noise is subdominant compared to white noise. We assume that this holds for 1/4th of the total dataset (indeed, there are many such pulsars  Cordes:2010fh ; Lam:2016iie ; NANOGrav:2020bcs ; Goncharov:2020krd ; NANOGrav:2020qll ).

  When the frequency of the signal is less than $1/T_{obs}$, there are two ways to do the analysis: i) One can keep the white noise fixed and expand the signal in a Taylor series. Then, the first two expansion terms are absorbed by the period and period derivative terms due to the lack of an independent measurement of these parameters, hence we are only left at the next order with the cubic term. Keeping the noise fixed, we therefore have an extra $m^{3}/T^{3}\sim f^{3}/T^{3}$ suppression. ii) One can keep the signal the same, but then the response will include again the ambiguity in period and period derivative, hence one can modify the definition of the noise via the transfer function Hazboun:2019vhv . The strain transfer function is unity for frequencies larger than inverse observation time ($f>1/T_{obs}$), on the other hand, and it is $f^{-3}$ at low frequencies ($f>1/T_{obs}$). The transfer function for power spectral amplitude is quadratic in the strain, hence they have $f^{-6}$. Red noise makes the sensitivity worse on top of this white noise curve.

The noise curve for a generic single pulsar is parametrized as in Ref. Unal:2020mts , consistent with both simulated curves Hazboun:2019vhv and data NANOGrav:2020qll

$$\!\!\!\!\!\!\!\!\!\!\!\!h_{c,{\rm noise}} = \sqrt{f\cdot S_{n}}=\sqrt{12\pi^{2}f^{3}\,({\cal N}/{\cal T})}\simeq 10^{-14}\times$$ (14)

$$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\times \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\;\sqrt{\frac{\Delta t_{14d}\;\sigma^{2}_{\mu s}}{{T_{\rm obs,15yr}^{3}}}}\;\left(\xi_{\rm frac}\cdot\left(f\cdot T_{\rm obs}\right)^{-3/2}+\left(f\cdot T_{\rm obs}\right)^{3/2}\right)$$

where $T_{\rm obs,15yr}=T_{\rm obs}/15{\rm yr}$ is scaled observation time, $\Delta t_{14d}=\frac{\Delta t}{14{\rm days}}$ is cadence, $\sigma_{\mu s}=\frac{\sigma}{{\rm\mu sec}}$ is rms signal error , $\xi_{\rm frac}$ is the inverse of the square-root of the fraction of pulsars where red noise is subdominant compared to white noise. We take 1/4th of the pulsars to be white-noise dominated in the regime we focus on, so $\xi_{frac}=2$ in this study. We emphasize that Eq. (14) is a generic result for PTAs, hence plugging in distinct timing errors or observation periods, one can produce approximate sensitivity curves. The scalings correspond to current PTA experiments, namely the rms timing error is normalized for ${\rm\mu s}$. The next generation SKA experiment is expected to improve this result by an order of magnitude, i.e. to $30\,{\rm ns}$. Moreover, the number of observed pulsars with SKA will be about two orders of magnitude larger than the current PTA pulsar number Smits:2008cf . Combined, these two effects will improve the sensitivity of SKA by roughly two orders of magnitude, as shown in Fig. 1.

To proceed, we plug Eq. (7) into Eq. (11), using our source is monochromatic, with $f_{*}$ given in Eq. (1)

$${\rm SNR}^{2} \simeq 2\,(f_{*}\cdot T_{\rm obs})^{2}\;{\cal R}^{2}_{\beta}\;\left(\frac{h_{c,{\rm signal}}}{h_{c,{\rm noise}}}\right)^{4}\bigg{|}_{f=f_{*}}\;,$$ (15)

where ${\cal R}^{2}_{\beta}$ is the response given to different spin fields, i.e. scalar, vector and tensor, given in Eq. (12). With this result in hand, we are equipped with all the required items to compute the sensitivity of current and future PTA experiments to large portions of the parameter space of ULDM particles of each distinct nature (scalar, vector and tensor). Our results are shown in Figs. 2, and 3.

Figure 3: Zoomed in for scalar, vector (dark photon) and tensor (spin-2) (setting $\alpha\!=\!10^{-7},N_{p}^{PTA}=60,N_{p}^{SKA}=5K$, and 15 years of observation time for both).

To provide intuition for the results, we now discuss in detail the capabilities of PTAs to probe ULDM by focusing on the scaling properties of the signal and noise in each separate frequency regime in Figs. 1, 2 and 3 by assuming conservative noise curve.

• a) $f>1/T_{\rm obs}$ : For scalar and vectors we have $h_{\rm signal}\propto f^{-2}$, while for tensor(spin-2) $h_{\rm signal}\propto f^{-1}$; and $h_{\rm noise}\propto f^{3/2}$.

  $$S,V : \;{\rm SNR}^{2}\propto(f\cdot T_{obs})^{2}\,\frac{h_{c,s}^{4}}{h_{c,n}^{4}}\propto{\cal F}_{\rm ULDM}^{\,4}\;(f\cdot T_{\rm obs})^{-12}$$

  $$T : \;{\rm SNR}^{2}\propto(f\cdot T_{obs})^{2}\,\frac{h_{c,s}^{4}}{h_{c,n}^{4}}\propto{\cal F}_{\rm ULDM}^{\,2}\;(f\cdot T_{\rm obs})^{-8}\;\;\;\;\;\;\;\;\;\;$$ (16)

  The fraction of dark matter we can probe scales as ${\cal F}_{\rm ULDM}\propto(f\cdot T_{\rm obs})^{3}$ for scalar and vector, and ${\cal F}_{\rm ULDM}\propto(f\cdot T_{\rm obs})^{4}$ for tensor. This is the regime for masses larger than $10^{-23}\,{\rm eV}$.

• b) $1/D_{\rm pulsar}<f<1/T_{\rm obs}$: For scalar and vectors we have $h_{\rm signal}\propto f^{-2}$, while for tensor $h_{\rm signal}\propto f^{-1}$; and $h_{\rm noise}\propto f^{-3/2}$.

  $$S,V : \;{\rm SNR}^{2}\propto(f\cdot T_{obs})^{2}\,\frac{h_{c,s}^{4}}{h_{c,n}^{4}}\propto{\cal F}_{\rm ULDM}^{\,4}\;(f\cdot T_{\rm obs})^{0}$$

  $$T : \;{\rm SNR}^{2}\propto(f\cdot T_{obs})^{2}\,\frac{h_{c,s}^{4}}{h_{c,n}^{4}}\propto{\cal F}_{\rm ULDM}^{\,2}\;(f\cdot T_{\rm obs})^{4}\;\;\;\;\;\;\;\;\;\;$$ (17)

  For scalars and vectors, in the $f<1/T_{\rm obs}$ regime, the signal grows slightly faster than noise, which results in ${\cal F}_{\rm ULDM}\propto(f\cdot T_{\rm obs})^{0}$. In the tensor case, we have ${\cal F}_{\rm ULDM}\propto(f\cdot T_{\rm obs})^{-2}$. Depending on $T_{\rm obs}$ of the PTA, this regime typically lies between $m\sim 10^{-26}\,{\rm eV}$ and $m\sim 10^{-23}\,{\rm eV}$.

• c) $f<1/D_{\rm pulsar}$: For scalar and vectors $h_{\rm signal}\propto f^{-1}$, while for tensor $h_{\rm signal}\propto f^{0}$; and $h_{\rm noise}\propto f^{-3/2}$.

  $$S,V : \;{\rm SNR}^{2}\propto(f\cdot T_{obs})^{2}\,\frac{h_{c,s}^{4}}{h_{c,n}^{4}}\propto{\cal F}_{\rm ULDM}^{\,4}\;(f\cdot D_{\rm pulsar})^{4}$$

  $$T : \;{\rm SNR}^{2}\propto(f\cdot T_{obs})^{2}\,\frac{h_{c,s}^{4}}{h_{c,n}^{4}}\propto{\cal F}_{\rm ULDM}^{\,2}\;(f\cdot D_{\rm pulsar})^{8}\;\;\;\;\;\;\;\;$$ (18)

  In the regime $f<1/D_{pulsar}$, the ${\rm SNR}$ decreases for scalar, vector and tensor. The fraction of dark matter that can be probed scales as ${\cal F}_{\rm ULDM}\propto(f\cdot D_{\rm pulsar})^{-1}$ for scalar and vector; and ${\cal F}_{\rm ULDM}\propto(f\cdot D_{\rm pulsar})^{-4}$ for tensor. Depending on the pulsar distance $D_{\rm pulsar}$, this is typically the regime where $m\lesssim 10^{-26}$eV.

To conclude, our results indicate that current as well as future PTAs can probe ultralight bosons of scalar, vector (dark photon) and tensor (spin-2) types with excellent precision. The mass range $10^{-26}-10^{-23}\,{\rm eV}$ is especially interesting since CMB and large-scale structure experiments have less sensitivity in that regime, and in contrast PTAs have the most sensitivity there, so that we can close this gap, as shown in Fig. 2 and 3. Crucially, this mass range roughly corresponds to the frequency regime $1/D_{\rm pulsar}\!<\!f\!<\!1/T_{\rm obs}$, where the signal stays strong compared to the noise. Therefore, smaller mass particles (in the scalar and vector scenarios) can be tightly constrained until frequencies $f\!\sim\!1/D_{\rm pulsar}$. For tensor particles, the potential oscillations scale with $1/f$, hence the best regime to probe them is around $f\!\sim\!1/T_{\rm obs}$.

Our calculations show that with current PTA data, the abundance of ultralight bosons in the mass range $10^{-26}-10^{-23}\,{\rm eV}$ can be probed down to ${\cal O}(1-10)\%$ of the total dark matter energy density. We also found that with 30 year PTA data, the precision improves by about one order of magnitude to $1\%$, and with SKA to $0.1\%$.

This work implies that combining PTAs with current constraints from CMB, large-scale structure, Lyman-$\alpha$ and superradiance (see Fig. 2), ultralight scalar dark matter can be constrained throughout the mass range $10^{-30}\!-\!10^{-17}\,{\rm eV}$ to less than $\mathcal{O}(10\%)$ of the dark matter.

Acknowledgements. We thank Jeff Hazboun, Sarah Libanore, David Marsh, Andrea Mitridate, Joseph Romano, Debanjan Sarkar, Marc Kamionkowski, Michael Lam, Sarah Vigeland, Stephen Taylor for discussions and feedback on the manuscript, and Jordan Flitter for the SPARC bounds. We thank especially Tristan Smith for the detailed explanations on low frequency signal, sensitivity and feedback on the draft. We also thank the anonymous referees for suggesting important corrections and improvements. CÜ thanks his family for their support during this work, and dedicates this work to Ece Ceyda Güdemek and Rümeysa Berin Şen. CÜ is supported by the Kreitman fellowship of BGU, and the Excellence fellowship of the Israeli Academy of Sciences and Humanities, and the Council for Higher Education. FU is supported by the European Regional Development Fund (ESIF/ERDF) and the Czech Ministry of Education, Youth and Sports (MEYS) through Project CoGraDS-CZ.02.1.01/0.0/0.0/15_003/0000437. EDK is supported by an Azrieli Foundation Faculty Fellowship.