Constraining Fundamental Constant Variations from Ultralight Dark Matter with Pulsar Timing Arrays

Pulsar Timing Arrays (PTAs) are exceptionally sensitive detectors in the frequency band $\text{nHz} \lesssim f \lesssim \mu\text{Hz}$. Ultralight dark matter (ULDM), with mass in the range $10^{-23}\,\text{eV} \lesssim m_\phi \lesssim 10^{-20}\,\text{eV}$, is one class of DM models known to generate signals in this frequency window. While purely gravitational signatures of ULDM have been studied previously, in this work we consider two signals in PTAs which arise in presence of direct couplings between ULDM and ordinary matter. These couplings induce variations in fundamental constants, i.e., particle masses and couplings. These variations can alter the moment of inertia of pulsars, inducing pulsar spin fluctuations via conservation of angular momentum, or induce apparent timing residuals due to reference clock shifts. By using mock data mimicking current PTA datasets, we show that PTA experiments outperform torsion balance and atomic clock constraints for ULDM coupled to electrons, muons, or gluons. In the case of coupling to quarks or photons, we find that PTAs and atomic clocks set similar constraints. Additionally, we discuss how future PTAs can further improve these constraints, and detail the unique properties of these signals relative to the previously studied effects of ULDM on PTAs.


I. INTRODUCTION
Identifying the nature of dark matter (DM) is one of the most important goals of physics. While roughly five times more abundant than visible matter, shockingly little is known about how DM fits in the Standard Model (SM) of particle physics. One class of DM models, known * david.kaplan@jhu.edu † amitri@caltech.edu ‡ ttrickle@caltech.edu as ultralight DM (ULDM), considers DM candidates with tiny particle masses, typically m φ 10 −6 eV. Such light DM candidates must be bosonic, otherwise they would severely violate the Tremaine-Gunn bound [1][2][3][4] on light fermionic DM.
These light DM candidates are well described by classical fields oscillating in time with a frequency 2πf m φ . In the presence of a direct coupling between DM and the SM, these oscillations behave as time-dependent sources for SM fields, and can lead to time-dependent signals with ω ∼ m φ . Therefore Pulsar Timing Arrays (PTAs), thanks to their tremendous ability to detect signals in the nHz − µHz frequency range, offer a unique probe of ULDM candidates in the mass window 10 −23 eV − 10 −20 eV. 1 Models of ULDM in this mass range are theoretically interesting since they are able to suppress structure formation up to length scales much larger than a prototypical WIMP particle [5]. This feature of ULDM models has been proposed as a possible solution of the small-scale challenges faced by the ΛCDM paradigm, such as the missing satellites, and core-cusp problem. See Ref. [6] for a recent review.
CMB observables [15][16][17] rule out ULDM with m φ < ∼ 10 −24 eV while recent studies of the Lyman-α forest [5,24,25], Milky Way (sub)Halo mass functions [22,23], and rotation curves [20,21] assert constraints of 10 −21 eV < ∼ m φ . However non-CMB constraints are susceptible to uncertainties related to determining properties of small scale structures, for which cosmological simulations [19,45], and analytic methods [46] can be insufficient to derive bounds with high precision. Therefore it is important to have complementary probes in this mass region which are not subject to the same uncertainties.
In this work we focus on the effect of ULDM induced fluctuations of fundamental constants. These fluctuations, in the mass window considered here, have previously been searched for before via their effects on atomic clock systems [30][31][32][33]35] and torsion balance experiments [43,44]. Here we detail two experimental signatures in PTAs first discussed in Ref. [36]. The first is due to fluctuations in the PTA reference clock. PTAs are exceptionally sensitive to timing residuals, and fluctuations in the reference clock used can be observed via these residuals. The second is from pulsar spin fluctuations. ULDM induced particle mass fluctuations will cause the pulsar mass to fluctuate. Therefore, by conservation of angular momentum, this must be accompanied by a fluctuation in the pulsar spin. We set constraints on these signals by using mock data closely resembling the current IPTA second data release, and a more futuristic PTA dataset. We find that in the mass range 10 −23 eV < ∼ m φ < ∼ 10 −20 eV current PTA constraints are competitive or stronger than current atomic clocks and fifth force constraints for a wide range of ULDM models.
The outline of this paper is as follows, in Sec. II we begin by reviewing the ULDM models which give rise to fundamental constant fluctuations, and then show how these models generate a PTA signal by inducing pulsar spin fluctuations, Sec. II A, and reference clock shifts, Sec. II B. In Sec. III we describe the analyses performed to set constraints, and the procedure followed to generate the mock data. In Sec. IV we discuss the results of the analyses, and the connection between the signals searched for here and those previously searched for in PTAs (Sec. IV A).

II. PTA SIGNALS FROM ULDM INDUCED FUNDAMENTAL CONSTANT VARIATIONS
As a model for ULDM we consider a singlet scalar field, φ, which constitutes all the cosmic DM. Owing to their high phase space density, cold (v φ ∼ 10 −3 ) ULDM candidates are well described by a classical field oscillating in time with a frequency ω = m φ , whereφ( x) is a random amplitude with zero mean and unit variance, ρ φ = 0.4 GeV cm −3 is the local ULDM density, and γ( x) is the phase of the field. This ULDM field can directly couple to ordinary matter in a plethora of ways. Following the notation in Ref. [8], we parameterize the couplings to the QED sector as where Λ = M Pl / √ 4π, M Pl is the Planck mass, and d's are dimensionless coupling coefficients. 2 Due to these couplings, fluctuations in the ULDM background induce variations in the fundamental constants of the SM. Specifically, the couplings in Eq. (2) drive fluctuations in the electromagnetic coupling constant, α, and the electron and muon masses, m e,µ , ULDM can also couple to the QCD sector, however, more care must be taken relative to the QED case due to the running of masses and couplings (implicitly ignored for QED). Still following the notation of Ref. [8], we parameterize the couplings to the QCD sector as where β 3 is the QCD beta function, and γ q are the light quark anomalous dimensions. This specific parameterization of the couplings is useful since it makes the φ induced shift to fundamental constants independent of β 3 and γ q . Specifically, the light quark mass shifts are given by while the shifts of the nucleon masses are δm p,n m p,n where C n = 0.048 [8], and we have defined the symmetric combination of the quark mass couplings as The first term in Eq. (6) comes from the shift to the QCD scale, δΛ QCD /Λ QCD = d g (φ/Λ). The second, subleading, term comes from the shift to the light quark masses, Eq. (5). The contribution from a strange quark mass shift is expected to be smaller by a factor of ∼ 4 [34], although there are uncertainties in deriving its precise value [8]. The effect of the other d's are even smaller [8], and neglected here for simplicity, though in principle these searches can put constraints on those parameters as well.
With the effects of ULDM on fundamental constants defined, we now describe how these can affect PTA observables. Pulsars are rotating, highly magnetized neutron stars that emit beams of electromagnetic radiation from their magnetic poles. Given a misalignment between the rotation and magnetic axes, the pulsar rotation can cause the radiation beam to sweep across Earth. If this happens, a pulsar will appear to an Earth observer as a periodic emitter. Thanks to pulsars' extremely stable rotation periods, the time of arrival (TOA) of these radiation pulses can be predicted with great accuracy. PTAs accurately measure the TOAs by looking for deviations from these predictions, a quantity referred to as timing residuals. In the following sections we will discuss how the fundamental constant variations can source these timing residuals and therefore be detected, or constrained, by PTAs.
Specifically, in Sec. II A we describe how mass fluctuations can induce pulsars spin variations which naturally source timing residuals. In Sec. II B we describe how shifts in fundamental constants can induce shifts in the energy levels of the atomic clocks used by PTAs, and lead to apparent aberrations in the pulsars TOAs.

A. Pulsar Spin Fluctuations
We begin by studying pulsar spin fluctuations generated by changes in the pulsar moment of inertia, I. By conservation of angular momentum, these moment of inertia fluctuations will induce corresponding fluctuations in the pulsar spin frequency, ω, Fluctuations in the pulsar constituent particle (dominantly neutrons, protons, electrons, and muons) masses can induce fluctuations in I. This dependence comes implicitly via the pulsar mass, M , as well as explicitly via the neutron mass which controls the balance of the Fermi degeneracy and gravitational pressures. In the simplest model of a spherically symmetric, non-rotating neutron star consisting of only non-relativistic neutrons one can show that R ∝ M −1/3 m −8/3 n , and therefore, . This is obviously a simplistic description, and a more realistic model for the mass dependence of I would require numerically solving the Tolman-Oppenheimer-Volkov equations [47], and include relativistic corrections. These corrections will induce O(1) deviations from the naive scaling; because of this we introduce the η and δη, pulsar dependent, parameters such that In the following, for simplicity, we will assume that η = 1/3, δη = −16/3, the values computed from the simplest model. These fluctuations in the pulsar spin frequency induce corresponding fluctuations, h, in the pulsar timing residuals via, From here it is clear how to understand the effects of ULDM models on pulsar timing; one just has to relate the particle mass fluctuations, δm f /m f , to pulsar mass fluctuations, δM/M 0 , where Y f ≡ N f /(N n + N p ) is the relative number of f particles to protons and neutrons in the pulsar. While these values will vary from pulsar to pulsar, we take Y n ∼ 0.9, Y p ∼ 0.1, Y e = Y p and Y µ ∼ 0.05 as fiducial values following Ref. [48].
We can now substitute Eqs. (3), (6) in to Eq. (11) to finally obtain the ULDM induced timing residuals in terms of the DM background field: where x P is the position of the pulsar, and y · d ≡ i y i d i where i ∈ {e, µ, γ, g,m} and we have introduced the sensitivity parameters of this search: One noteworthy feature is that while muons are less abundant than electrons, due to their larger mass the sensitivity is greater.

B. Reference Clock Shifts
The TOAs measured by PTAs are referenced to the Temps Atomique International (TAI) realization of the Terrestrial Time (TT). The TAI estimates TT using measurements from an ensemble of more than 400 atomic clocks (most of them based on the hyperfine transition of the ground state of the Cesium atom). A shift in the frequency of these clocks will induce an apparent shift in the pulsar spin frequency, which translates to pulsar timing residuals according to the first equality in Eq. (10). In this section, we will discuss how fundamental constant variations can induce such shifts.
Using a similar notation to Ref. [34], the scaling of the frequency of an atomic clock is given by where F rel (Zα) is the relativistic correction to the energy levels of an atom with nuclear charge Z, µ is the nuclear magnetic moment, and ζ = 1 (ζ = 0) for clocks using hyperfine (optical) transitions. It is clear then that fluctuations in fundamental constants will induce fluctuations in atomic clock frequencies. Specifically, for a clock using an atom A, where δF rel /F rel = K A δα/α, and δµ/µ = C A δm q /m q . For the case of Cesium atoms, K A = 0.83, and C A = 0.110 [34]. Finally, by substituting Eqs. (3), (6) in to Eq. (15) we obtain the induced timing residuals in terms of the ULDM background: where x E is the position of the Earth, and the sensitivity parameters of this search are given by: where While the signals in Eq. (12) and Eq. (16) appear similar, they have a few distinguishing features. First, notice the difference in the sensitivity parameters, y in Eqs. (17), (13). The pulsar spin fluctuation signal is dominant for muon-philic models, while the reference clock shift is dominant for models coupling to photons or electrons. The two signals also differ in the correlations between pulsars. Assuming the distance between pulsars is much larger than the DM coherence length, the pulsar spin fluctuations will be uncorrelated between pulsars, δ IJ . Conversely, the reference clock shift affects all of the pulsars in an identical way; this will leave an imprint of monopole correlations across the array, or h I h J ∝ 1. Note that both of these signals are distinct from those of gravitational waves or Doppler shifts which will leave quadrupole, and dipole signatures, respectively. A summary of the difference between these effects and other ULDM effects in PTAs is given in Table III in Sec. IV A.

III. DATA ANALYSIS
In this section, we discuss the data analysis techniques used to constrain ULDM couplings. The analysis will closely follow the standard procedure adopted by PTA collaborations. This section aims to summarize how this procedure can be applied in our context. See Ref. [49] for a more detailed discussion on Bayesian inference with PTA data, as well as Ref. [50] for a more pedagogical discussion.

A. Noise Modeling and the PTA Likelihood
The main observable in a PTA experiment is the timing residuals, δt, which measure the discrepancy between the observed times of arrival (TOAs) and the ones predicted by the pulsar timing model every cadence, ∆t. Generally, there are three main contributions to these timing residuals: white noise, red noise, and small errors in the fit to the timing-ephemeris parameters. Specifically, we can model the timing residuals as: We will now discuss each of these three terms in more detail. For each of the N TOAs, the white noise is assumed to be a normally distributed random variable, with zero mean and variance, where σ i is the TOA uncertainty for the i−th observation, E µ is the Extra FACtor (EFAC) parameter for the receiver-backend system µ, and Q µ is the Extra QUADreature (EQUAD) parameter. 3 Red noise is modeled in a Fourier basis for frequencies k/T , where k indexes the harmonics of the basis and T is the observation time. Since red noise is a low-frequency process, the summation over k is truncated to some reasonable frequency, N f . For example, in the most recent NANOGrav search for a gravitational wave background (GWB) this cutoff was set to N f = 30. We will use the same cutoff in our analyses. This set of 2N f sine-cosine pairs evaluated at the different observation times is contained in the Fourier density matrix, F. The Fourier coefficients a are assumed to be normally distributed random variables, with zero mean and covariance matrix a a T = φ. This covariance matrix will contain all the possible sources of low-frequency achromatic noise. For our analysis we will consider two possible sources: pulsar intrinsic red noise, and GWB. Therefore, φ takes the form where a and b index pulsars, k and j index frequency harmonics, and Γ ab is the GWB overlap function for the pulsar pair (a, b). In this expression the terms ρ k (κ ak ) are related to the power spectral density (PSD) of the timing residuals, S(f ), induced by the GWB (pulsar intrinsic red noise) such that ρ(f ) = S(f )∆f with ∆f = 1/T (and similarly for κ a (f )). Following the assumptions commonly made in PTA analysis, we model the PSD for the two red noise contributions as Finally, possible small deviations from the initial bestfit values of the m timing-ephemeris parameters are accounted for by the term M . The design matrix, M , is an N × m matrix containing the partial derivatives of the TOAs with respect to each timing-ephemeris parameter (evaluated at the initial best-fit value), and is a vector containing the linear offset from these best-fit parameters.
Since in most analyses, ours included, we do not care about the specific realization of the noise but are instead interested in its statistical properties, we can analytically marginalize over all the possible noise realizations (i.e., integrate over all the possible values of a and ). This leaves us with a marginalized likelihood that depends only on the hyper-parameters, η = (A GWB , γ, . . .), characterizing the statistical properties of the noise [53,54]: where C = N +T BT T . Here N is the covariance matrix of the white noise, B = diag(∞, φ), and T = [M , F ].

B. Signal Parameterization
The marginalized likelihood in Eq. (23) can be easily generalized to take into account deterministic signals in the timing residuals. In presence of a deterministic signal, h( θ), which depends on the parameters θ, we just need to make the replacement δt → δt − h( θ) in Eq. (23). In our case the deterministic signal is given by the ULDM induced fluctuations described in Sec. II.
Generally, the ULDM signals for the I th pulsar in the array will have the form where A i = d i 2ρ φ /(m φ Λ) here, and we have defined the Earth and pulsars phases as where k ∼ m φ v φ is the characteristic DM momentum.
We will now briefly discuss how we treat the parameters appearing in Eq. (24). Since the observation times of interest here are T ∼ 10 year, the motion of celestial bodies is of order O(10 −3 pc) and we can safely ignore the variation in k · x due to the Earth/pulsar motion and take γ E,P to be constants. Moreover, since typical values of | x P − x E | range between 0.1 and several kpc, the term m φ | x E − x P | is never negligible for the DM masses of interest here. Therefore, γ P and γ E need to be treated as independent parameters. The DM momentum sets the coherence length, l c , of the DM background, The DM background amplitude,φ( x), is correlated within a coherence length, and uncorrelated outside of it.
Since the typical pulsar-pulsar and pulsar-Earth separations are comparable with the ULDM coherence length, following the same practice of [39], we perform our analysis in two limits: • Correlated: in this case we will have oneφ parameter for both the Earth and all the pulsars terms: φ E =φ P,I for all I • Uncorrelated: in this caseφ E and all theφ P,I will be uncorrelated free parameters The correlated analysis is expected to provide reliable results for the low-mass region considered in this work, and vice-versa for the uncorrelated analysis.

C. Setting Constraints
To set constraints, we derive the Bayesian posteriors distributions for the noise and DM parameters by using the Markov Chain Monte Carlo (MCMC) techniques implemented in the PTMCMCSampler package [55], together with the marginalized likelihood given in Eq. (23) (implemented by using the enterprise [56] and enterprise extensions [57] packages), and the priors distributions given in Table I. We then marginalize over all parameters except the DM signal amplitude, A i , and the DM mass, m φ . The constraints on the signal amplitude in each mass bin are then set to the 95 th percentile of the amplitude in that bin. These constraints can then be easily translated to constraints on the DM couplings d i .

D. Mock Data
While future analyses will be performed using NANOGrav and IPTA data, in this work we will set constraints by using mock data that should closely resemble the current IPTA DR2 dataset [58] and future PTAs data [59][60][61]. The mock data are generated by using the pulsar parameters contained in the .par files of the International Pulsar Timing Array (IPTA) second Mock Data Challenge (MDC2) [62], the python wrapper libstempo [63] and the pulsar timing package TEMPO2 [64,65].
For each pulsar we generate a list of TOAs by using the PTA parameters (i.e. cadence, ∆t, TOAs uncertainties, σ i , number of pulsars, N P , and observation time, T ) given in Table II. The TOAs are injected with the noise sources described in the previous sections: white noise, pulsar intrinsic red noise, and GWB. For each pulsar, the values of the noise parameters are randomly sampled from the priors given in Table I.

IV. RESULTS
Applying the analysis detailed in Sec. III to the signals discussed in Sec. II, we place projected constraints on the  Table II. Dashed (solid) lines correspond to searches for an (un)correlated signal (labelled in "unco.", "corr.", respectively), as discussed in Sec. III. In deriving these plots we have set η = 1/3, δη = −16/3, and assumed that PTA clocks use the transition between the two hyperfine levels of the 133 Cs ground state. The lower gray shaded regions correspond to regions of parameter space where the signal amplitude is less than the purely gravitational signal. Current constraints "Rb/Cs atomic clocks" (purple) are from Ref. [31], "Al/Hg atomic clocks" (turquoise) are from Ref. [35], "MICROSCOPE" (teal) are from Ref. [44], "H/Si clock shift" (orange) are from Ref. [30], and "NS binary system" are from Refs. [41,42].
five scalar ULDM coupling constants: dm, d e , d µ , d γ and d g . The results are shown in Figs. 1, where we compare them to constraints from atomic clocks [30,31,35], torsion balances [44], and decay of neutron star binaries orbits [41,42]. While other constraints apply to these models (see, e.g., Refs. [30,44,66]) to simplify Fig. 1 we only show those which have the strongest constraints in the mass range of interest here: 10 −24 eV < ∼ m φ < ∼ 10 −20 eV. We find that, in the mass range under consideration, an IPTA-like array can already be competitive with the most precise atomic clock experiments [30,31,35] available today when ULDM couples to the QCD sector via dm, and strictly outperforms them when the coupling happens through d g . In the QED sector, atomic clock experiments are slightly better at constraining d γ . 4 However, they are typically not able to constrain d e since relative frequency shifts between clocks, using the same type of atomic transition, cause the d e dependence to drop out. To avoid this cancellation one needs to compare different types of transitions as done in [30]. At the moment, the observation time of this experiment limits the constraints to 10 −21 eV < ∼ m φ . PTAs can also beat torsion balance constraints [43,44] by several orders of magnitude. Lastly, we note that PTAs are exceptionally sensitive to d µ . This is due to the large number of muons within neutron stars, sustained through β equilibrium [67,68]. This allows PTAs to improve by more than three orders of magnitude the next best projected constraints from the orbital decay of a neutron star binary system [41,42].
All of the previous conclusions are amplified for future PTA experiments. With better timing parameters lead to O(1) differences in the projections [14]. For simplicity, we directly reproduce the constraints from their respective paper and do not rescale them. and more pulsars, a future PTA will be a strict improvement over the IPTA, as can be seen in Fig. 1. 5 Similar improvements could also be achieved by the IPTA collaborations, as the observation time and number of pulsars keep increasing. To gain intuition for the dependence of these constraints on the PTA parameters we will momentarily resort to a frequentist signal analysis approach, see Refs. [69][70][71][72] for more details. Given a deterministic signal in each pulsar, h I = A sin (m φ t + γ) for simplicity, and assuming N P identical pulsars with only white noise, the signal SNR can be shown to scale as, This explains the scaling of the constraints in Fig. 1 at high masses; SNR ∝ A ∝ d/m 2 φ and therefore the constraints on d are proportional to m 2 φ . The scaling of Eq. (28) is appropriate when the signal is known in all the pulsars, i.e., for an Earth term signal, since this SNR is derived using a matched filter approach. A pulsar dependent signal would have a more complicated, weaker, dependence on N P , with all other dependencies being equal. With the scaling in Eq. (28) future PTAs are expected to perform better than IPTA-like ones by a factor, For the specific future PTA this means a predicted improvement of a factor of ∼ 5, in good agreement with Fig. 1.
At low masses, m φ < ∼ 1/T , this picture breaks down. In this limit, the signal no longer oscillates over the observing time and looks like a polynomial expanded in m φ t. PTAs are not sensitive to the first few terms, up to O(t 2 ), in this polynomial since they are degenerate with the timing model [69]. The first term they are sensitive to is ∝ A(m φ t) 3 , and therefore at low masses, SNR ∝ A(m φ T ) 3 and the constraints weaken as d ∝ 1/m φ . In addition to these subtraction effects, red noise can also deteriorate the signal significance [73] at low masses. This explains why the constraints begin to flatten before turning to the d ∝ 1/m φ scaling.
Lastly, we note that the reduced sensitivity around f ∼ year −1 is due to the fact that the relative position of the pulsar, which oscillates yearly, is not known to high precision [64,65] and must be fit away, therefore reducing sensitivity. 5 The increased dimensionality of the parameter space for the future PTA setup makes sampling computationally demanding. Because of this, some constraints are less converged, specifically the correlated constraints on dm, dµ, and dg. In Fig. 1, these constraints have been smoothed by averaging the constraints from neighboring mass points.

A. Connections to Other ULDM Effects in PTAs
We will conclude this section by discussing the specific differences in the ULDM induced signals presented here versus those considered previously [26-29, 38, 39]. A summary of our discussion is given in Table III, previous effects are labeled "Shapiro" and "Doppler", and the new effects are labeled "Pulsar Spin Fluctuations" and "Reference Clock Shift". The primary purpose of most PTAs is to measure the stochastic gravitational wave background. These ripples in space-time affect the photons' geodesic when traveling from the pulsar to the Earth. ULDM can induce a similar effect; the local ULDM density sources fluctuations to the metric which can affect light travel time in near complete analogy with the gravitational waves PTAs are searching for. This effect only depends on the background ULDM density and is model independent; any ULDM model which constitutes all of the DM will produce this effect. The amplitude of the signal is given by, see Ref. [26] for more details of its derivation. To justify ignoring it in our previous analysis we require that the amplitude of the effects considered here is greater than the amplitude in Eq. (30). This happens when and explains the gray shaded region in the bottom left corner of Figs. 1. In this region the purely gravitational signal is stronger and a more thorough analysis should include all effects. However these signals are not degenerate for a few reasons. First, the purely gravitational signal has a frequency of ω = 2m φ whereas the model dependent ULDM signals have a frequency of ω = m φ . In addition, the pulsar spin fluctuation signal is uncorrelated between pulsars, whereas the purely gravitational one is perfectly correlated giving another method of disentangling the signals.
Recently, the Doppler effect from vector ULDM has been studied using the Parkes PTA dataset [39]. This effect can be understood by thinking about the vector DM sourcing a dark "electric" field from ∂ t φ. If the pulsar or Earth are charged under this dark force they will accelerate in the presence of this field. A similar effect is present for scalar ULDM, however the acceleration induced will be due to the spatial derivative of the field, a ∼ ∇φ, leading to a suppression since v φ ∼ 10 −3 . Therefore this effect will be subdominant to those discussed here which are not velocity suppressed. However the signal correlations are dipolar between pulsars which could potentially be used to search for this smaller signal.
Eq. (16) TABLE III. Summary of effects in PTA timing residuals from scalar ULDM candidates. A is the signal amplitude, P/E denotes whether there are "Pulsar" or "Earth" terms in the signal, assuming the signal in each pulsar is written as h = AP sin(ω + γP ) + AE sin(ω + γE). Note that the sensitivity parameters, y, will depend on whether the term comes from a pulsar or Earth effect. ω is the signal frequency, and hI hJ represents the signal correlation between pulsars. Effects which are velocity suppressed are highlighted in red.

V. CONCLUSIONS
While the primary goal of PTAs is to detect the stochastic gravitational wave background (SGWB), their remarkable sensitivity can be leveraged for new physics searches. From searching for novel contributions to the GWB produced via cosmic strings [74,75], inflationary fluctuations [76][77][78], and cosmological phase transitions [79][80][81][82][83][84] to signals produced from passing DM substructure [69,73,85,86] PTAs have been shown to set quite powerful constraints. In this work we have discussed two types of signals generically expected from scalar ULDM models with direct couplings to the SM. These couplings will generate apparent fluctuations in the SM fundamental constants, e.g., particle masses and couplings, which are then realized as new contributions to the timing residuals.
The two signals we discussed were: "pulsar spin fluctuations", Sec. II A, generated via changes to the moment of inertia and conservation of angular momentum, and "reference clock shifts", Sec. II B which change the tick frequency of the reference clocks used for the PTA. Notably, neither of these effects suffer from a velocity suppression as a Doppler shift would. We find that for a wide range of ULDM models the IPTA, and future PTA, can compete or outperform other constraints in the 10 −23 eV < ∼ m φ < ∼ 10 −20 eV mass window such as other atomic clock experiments, torsion balances, and the orbital decay of neutron star binaries.
Despite the fact that the results derived in this paper were obtained using mock data, the formalism and code employed are ready to be deployed on real PTA datasets. Having proved the constraining power of PTAs for these kinds of ULDM models, we plan to apply the analyses performed in this paper to NANOGrav and IPTA data in the near future.