© The Authors 2024

1. Introduction

Binary pulsars, namely astrophysical binary systems in which at least one of the two orbiting bodies is a pulsar, are outstandingly precise clocks, and as such they have been used in a wide variety of applications–readers can refer to Lorimer (2008) for a comprehensive review. Among those applications, binary pulsars can act as potential detectors for dark matter whose mass lies in the lower end of the ultra-light dark matter (ULDM) range: 10⁻²³ eV ≲ m ≲ 10⁻¹⁸ eV. In the late Universe, ULDM with such a mass behaves like a coherent superposition of oscillators all with the same inverse frequency of about 1.15 hours × (10⁻¹⁸ eV/m), which roughly corresponds to the span of periods of observed binary pulsars. Two interesting phenomena take place when the orbital period of the system, or any of its harmonics, is in resonance with the ULDM oscillations. First, because the ULDM oscillations perturb the gravitational potential in which the binary system is immersed, they cause secular drifts in the orbital motion, which however are unlikely to be observed. Second, if the ULDM carries a (fifth) force that couples directly to the bodies of the binary system, it brings about further secular drifts in the orbital parameters that can be potentially detected¹. Such resonant effects have been studied for ULDM of spin 0, spin 1 and spin 2 Blas et al. (2017, 2020), López Nacir & Urban (2018), Armaleo et al. (2020a), and in all cases binary pulsars have been shown to be promising detectors for such kinds of ULDM

These results however are limiting in three respects. First, resonances appear only for a very narrow range of ULDM mass m peaked at a given harmonic of the orbital period; for the linear direct coupling we consider here m ≈ kω_b where $k\in \mathbb{N}$ is the harmonic number and ω_b = 2π/P_b is the orbital frequency and P_b the orbital period. This means that each binary system can probe a very small parameter space in ULDM mass, or, in other words, that only very few, if any, binary pulsars are sensitive to ULDM. Second, for the most-studied case of a universally-coupled spin 0, the resonant secular drift on the period P_b–which is often the best-measured parameter–vanishes for systems with small eccentricities e → 0, which are the majority (e.g. more than 80% in the ATNF catalogue Manchester et al. (2005) have e < 0.001)². Hence, in this case we are limited to an even smaller sub-sample of all observed binary pulsars. Third, so far the sensitivity of binary pulsars to ULDM has only been estimated by requiring that the ULDM effects be within the errors of the measured orbital parameters, thus ignoring the time behaviour of the ULDM perturbations on the orbital parameters

In this work we develop a new method to estimate the sensitivity of binary pulsars as detectors of ULDM. For practical purposes here we limit ourselves to the universally-coupled linear scalar case, but our results are readily generalised to other spins and types of coupling. Our proposal draws from existing methods used to determine the sensitivity of pulsar-timing arrays to continuous gravitational waves, and takes advantage of the full time evolution of the ULDM perturbations to the binary system. In this way we are able to improve on existing results in three different ways. First, for a given system, we are able to compute semi-analytically the sensitivity to ULDM for all values of the ULDM mass, that is, also away from the resonance. Second, since we are not limited to resonances, the effects of ULDM do not vanish for e → 0, only the resonances do; hence, even for a universally-coupled scalar ULDM, all binary pulsars can act as ULDM detectors. Third, with our method we can straightforwardly combine multiple binary pulsars, each of which spans the entire mass range of interest. This vastly improves the sensitivity to ULDM and can tighten the constraint on the ULDM parameter space should no signal be found

The rest of this paper is organised as follows: in Sect. 2 we introduce the post-Keplerian, osculating-orbits formalism and derive the effects of ULDM on the six orbital parameters. In Sect. 3 we describe our method to estimate the sensitivity. We present our results in Sect. 4 and conclude with an outlook for future work in Sect. 5. In the Appendix we collect several useful formulas and additional tests and cross-checks of our method. Throughout the paper we use units for which c = 1 and ℏ = 1

2. The perturbed orbital motion

2.1. Timing models

The phase N of the pulsar at the moment at which the pulse is emitted is related to the proper time as measured by a hypothetical clock on the pulsar T by

$$\begin{array}{c}\hfill N=N_0+\nu T+\frac{1}{2}\frac{\mathrm{d}\nu }{\mathrm{d}T}{T}^2+\frac{1}{6}\frac{{\mathrm{d}}^2\nu }{\mathrm{d}{T}^2}{T}^3,\end{array}$$(1)

where ν is the spin frequency of the pulsar and N₀ an arbitrary constant–readers can refer to Lorimer & Kramer (2004). The pulsar’s proper time T is related to the time of arrival of each pulse at the detector on Earth through a timing model. In this work we focus on two timing models, the BT model of Blandford & Teukolsky (1976) and the ELL1 model of Lange et al. (2001), which better describes near-circular orbits

For non-relativistic orbits the relation between proper and Earth’s times can be written as

$$\begin{array}{c}\hfill t=T+{\alpha }_b(cosE-e)+{\beta }_{b}sinE,\end{array}$$(2)

where we defined

$$\begin{array}{cc}\hfill {\alpha }_b& \equiv xsin\omega ,{\beta }_b\equiv {(1-e^2)}^{1/2}xcos\omega ,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill x& \equiv a_{1}s,s\equiv sin\iota ,a_1\equiv \frac{M_2}{M_1+M_2}a,\hfill \end{array}$$(4)

where a is the semi-major axis of the binary, e its eccentricity, ι its inclination, ω the argument of periapsis (not to be confused with the orbital frequency ω_b), and M₁ and M₂ are the masses of the pulsars and its companion, respectively–the Euler angles ι and ω, together with the longitude of the ascending node ω define the rotation of the (x,y,z) system with respect to the observer’s system (X,Y,Z) (in which the observer is at $Z=-\infty$), we refer to Fig. 1³. The eccentric anomaly E satisfies Kepler’s equation

$$\begin{array}{c}\hfill E-esinE={\int }_{T_0}^T{\omega }_b\mathrm{d}T\prime \equiv \mathrm{\Theta },\end{array}$$(5)

where T₀ is the proper time of periapsis passage⁴. In order to make contact with the literature, we then define

$$\begin{array}{c}\hfill {\eta }_b\equiv {\beta }_b+{\gamma }_b,\end{array}$$(6)

where ${\gamma }_b$ is the Einstein delay; in practice this term is presently indistinguishable from β_b, which means that ${\eta }_b\approx {\beta }_b$. Lastly, in this work we always ignore the Shapiro delay term in the timing model

With the timing model at hand we can perturbatively solve for N as a function of t and the timing-model parameters to obtain

$$\begin{array}{cc}\hfill N& =N_0+\nu t+\frac{\dot{\nu }}{2}{t}^2+\frac{\ddot{\nu }}{6}{t}^3-\nu f(E\prime )-\dot{\nu }tf(E\prime )\hfill \\ \hfill & +\frac{\nu {\omega }_b}{1-ecosE\prime }f(z)\frac{df(z)}{\mathit{dz}}{|}_{z=E\prime },\hfill \end{array}$$(7)

where we defined the function

$$\begin{array}{c}\hfill f(z)\equiv {\alpha }_b(cosz-e)+{\eta }_{b}sinz,\end{array}$$(8)

and a new eccentric anomaly E′ as

$$\begin{array}{c}\hfill E\prime -esinE\prime ={\int }_{T_0}^t{\omega }_b\mathrm{d}t\prime \equiv \mathrm{\Theta }\prime .\end{array}$$(9)

Eq. (7) defines the BT timing model, for which we can equivalently use the parameter sets $\{a,e,\omega ,s,T_0\}$ or $\{a,e,{\alpha }_b,{\eta }_b,\mathrm{\Theta }\prime \}$. On top of these, there are the parameters $N_0$, ν, $\dot{\nu }$ and $\ddot{\nu }$ that are related to the pulsar’s own rotation (and are common to all timing models)

If the orbit is nearly circular, e → 0, T₀ and ω are ill-defined because there is no periapsis. In this case we can define a different set of variables, the Laplace-Lagrange parameters η (not to be confused with parameter ${\eta }_b$ defined in Eq. (6)) and κ, as

$$\begin{array}{cc}\hfill \eta & \equiv esin\omega ,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill \kappa & \equiv ecos\omega ,\hfill \end{array}$$(11)

and we can replace T₀ with T_asc, which stands for the time of ascending node. In analogy with the definition (5) we write

$$\begin{array}{c}\hfill \mathrm{\Psi }\equiv {\int }_{T_{\mathrm{asc}}}^T{\omega }_b\mathrm{d}T\prime =E-esinE+\omega ,\end{array}$$(12)

with $\omega ={\displaystyle {\int }_{T_{\mathrm{asc}}}^{T_0}}{\omega }_b\mathrm{d}T\prime$⁵. In order to make contact with the new eccentric anomaly E′ we can finally define a new variable Ψ′ through

$$\begin{array}{c}\hfill \mathrm{\Psi }\prime \equiv {\int }_{T_{\mathrm{asc}}}^t{\omega }_b\mathrm{d}t\prime =E\prime -esinE\prime +\omega .\end{array}$$(13)

With these parameters, the relationship (7) becomes

$$\begin{array}{cc}\hfill N& =N_0+\nu t+\frac{\dot{\nu }}{2}{t}^2+\frac{\ddot{\nu }}{6}{t}^3-\dot{\nu }xg(\mathrm{\Psi }\prime )t\hfill \\ \hfill & -\nu g(\mathrm{\Psi }\prime )+\nu {\omega }_{b}g(z)\frac{dg(z)}{\mathit{dz}}{|}_{z=\mathrm{\Psi }\prime },\hfill \end{array}$$(14)

where now

$$\begin{array}{c}\hfill g(z)\equiv -x\frac{\eta }{2}cos2z-x\frac{3\eta }{2}+x\frac{\kappa }{2}sin2z+xsinz.\end{array}$$(15)

Eq. (14) defines the ELL1 timing model Lange et al. (2001), that is, the set $\{x,\eta ,\kappa ,s,T_{\mathrm{asc}}\}$ or with $T_{\mathrm{asc}}\to \mathrm{\Psi }\prime$

2.2. Ultra-light dark matter as a perturbing force

The unperturbed orbital motion of two bodies in the binary is completely characterised by six orbital parameters: a, e, ι, ω, T₀ and the longitude of the ascending node Ω. The ULDM field $\mathrm{\Phi }(t)$, owing to its interaction with the binary system, perturbs the system’s orbital motion and modifies the timing models of (7) or (14) that describe it. To illustrate our method, here we consider an effective, universal, direct linear coupling between the ULDM field and the two bodies in the binary system. The coupling is characterised by a single parameter α for which we can perturbatively write

$$\begin{array}{c}\hfill M_A^{(\alpha )}(\mathrm{\Phi })=M_A(1+\alpha \mathrm{\Phi }),\end{array}$$(16)

where A ∈ {1, 2}, with M₁ being the pulsar’s mass and M₂ the mass of the companion⁶

The ULDM field can be considered homogeneous within its de Broglie wavelength (or coherence length), as determined by

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{dB}}\sim 1.3\times {10}^{12}\mathrm{km}\left(\frac{{10}^{-3}}{v_0}\right)\left(\frac{{10}^{-18}\mathrm{eV}}{m}\right),\end{array}$$(17)

where $v_0$ is the typical ULDM velocity in the Milky Way halo. However, its amplitude and phase are expected to show large oscillations from binary system to binary system. Therefore, we describe the ULDM field as Foster et al. (2018)

$$\begin{array}{c}\hfill \mathrm{\Phi }={\mathrm{\Phi }}_0\varrho cos(mt+\mathrm{{\rm Y}}),\end{array}$$(18)

where, assuming that each binary system is in a different coherence patch, Υ is a random phase with a flat distribution in [0, 2π) ,

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_0\equiv \frac{\sqrt{2{\rho }_{\mathrm{DM}}}}{m},\end{array}$$(19)

the local ULDM density is taken to be ρ_DM = 0.3 GeV/cm³ and 𝜚 ≥ 0 is a random variable drawn from the Rayleigh distribution

$$\begin{array}{c}\hfill P(\varrho )=2\varrho e^{-{\varrho }^2}.\end{array}$$(20)

The stochastic character of the ULDM field comes about because pulsar-timing observation campaigns typically last up to a few decades, which is a much shorter time than the ULDM coherence time

$$\begin{array}{c}\hfill t_{\mathrm{coh}}\sim 65\mathrm{y}{\left(\frac{{10}^{-3}}{v_0}\right)}^2\left(\frac{{10}^{-18}\mathrm{eV}}{m}\right);\end{array}$$(21)

indeed, only by sampling the field during several $t_{\mathrm{coh}}$ does the amplitude converge to $\varrho \to 1$

The perturbed, non-relativistic Keplerian equations of motion of a binary system interacting with ULDM are given by

$$\begin{array}{cc}\hfill \ddot{\mathit{R}}& =-\alpha \dot{\mathrm{\Phi }}\dot{\mathit{R}},\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill \ddot{\mathit{r}}& =-(1+\alpha \mathrm{\Phi })\frac{G{M}_T\mathit{r}}{r^3}-\alpha \dot{\mathrm{\Phi }}\dot{\mathit{r}}\hfill \\ \hfill & =-\frac{G{M}_T\mathit{r}}{r^3}+\mathit{F},\hfill \end{array}$$(23)

where a dot denotes a derivative with respect to (coordinate) time, M_T ≡ M₁ + M₂, R is the binary barycentre position defined using the unperturbed masses, r is the relative coordinate with origin in the companion body with mass M₂ (readers can refer to Blas et al. 2020) that defines the polar coordinate system $(\widehat{r},\widehat{\theta },\widehat{z}$) of the orbital motion. In the last equality we defined F which we refer to as the ‘perturbing force’

2.3. Osculating orbits

The perturbing force F defined in Eq. (23) leads to deviations from the Kepler motion that can be studied using the osculating orbits formalism as described in Danby (1970). After some algebra (more in Appendix A), choosing {a , e , ω , ϵ₁ , Ω , s} as the six independent orbital parameters, we obtain

$$\begin{array}{cc}\hfill \frac{\dot{a}}{a}& =-\frac{2\alpha \mathrm{\Phi }e{\omega }_b}{\sqrt{1-e^2}}\frac{a^2}{r^2}sin\theta \hfill \\ \hfill & -\frac{2\alpha \dot{\mathrm{\Phi }}(1+e^2+2ecos\theta )}{1-e^2},\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill \dot{e}& =-\alpha \mathrm{\Phi }\sqrt{1-e^2}{\omega }_b\frac{a^2}{r^2}sin\theta \hfill \\ \hfill & -2\alpha \dot{\mathrm{\Phi }}(cos\theta +e),\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill \dot{\omega }& =\frac{\alpha \mathrm{\Phi }{\omega }_b\sqrt{1-e^2}}{e}\frac{a^2}{r^2}cos\theta \hfill \\ \hfill & -\frac{2\alpha \dot{\mathrm{\Phi }}}{e}sin\theta ,\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill \dot{{ϵ}_1}& =\frac{2\alpha \mathrm{\Phi }{\omega }_b}{1-e^2}(1+ecos\theta )\hfill \\ \hfill & +\frac{2\alpha \dot{\mathrm{\Phi }}e}{\sqrt{1-e^2}}\frac{r}{a}sin\theta \hfill \\ \hfill & +[1-\sqrt{1-e^2}]\dot{\omega },\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill \dot{\mathrm{\Omega }}& =\dot{s}=0,\hfill \end{array}$$(28)

where the parameter ϵ₁ is given by

$$\begin{array}{c}\hfill {ϵ}_1\equiv {\omega }_b(t-T_0)+\omega +\mathrm{\Omega }-{\int }_{T_0}^t\mathrm{d}t\prime {\omega }_b,\end{array}$$(29)

and is related to Θ′=Ψ′−ω by $\dot{{ϵ}_1}=\dot{\mathrm{\Omega }}+\dot{\omega }+\dot{\mathrm{\Theta }}\prime -{\omega }_b$.⁷ We note that scalar ULDM does not have any effect on the parameters ω and s, which we can henceforth discard. Integrating $\dot{\mathrm{\Theta }}\prime$, using ${\omega }_b^2=G{M}_T/a^3$ and allowing for a possible error and secular evolution of the semi-major axis $a\to a+\dot{a}(t-T_{\mathrm{asc}})$, we find

$$\begin{array}{cc}\hfill \delta \mathrm{\Theta }\prime & =\delta \mathrm{\Theta }{\prime }_0\hfill \\ \hfill & +{\int }_{T_0}^t[\sqrt{\frac{G{M}_T(1+\alpha \mathrm{\Phi })}{{(a+\delta a)}^3}}-\sqrt{\frac{G{M}_T}{a^3}}]\mathrm{d}t\prime \hfill \\ \hfill & +2\alpha {\omega }_b{\int }_{T_0}^t\mathrm{\Phi }\mathrm{d}t\prime \hfill \\ \hfill & +\alpha \sum _{n=1}^{\infty }{\int }_{T_0}^t[A_n\mathrm{\Phi }cos(n{\omega }_b\tilde{t}\prime )+B_n\dot{\mathrm{\Phi }}sin(n{\omega }_b\tilde{t}\prime )]\mathrm{d}t\prime ,\hfill \end{array}$$(30)

where δΘ′₀ is the error in the determination of Θ′ at T₀, $\stackrel{\sim }{t}\prime \equiv t\prime -T_0$, with T₀ = T_asc + ω/ω_b, and

$$\begin{array}{cc}\hfill A_n& \equiv 4{\omega }_b{J}_n(ne)-2n{\omega }_b\frac{1-e^2}{e}{J}_n\prime (ne),\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill B_n& \equiv \frac{4}{n}{J}_n(ne)+4\frac{1-e^2}{e}{J}_n\prime (ne),\hfill \end{array}$$(32)

where $J_n(z)$ is a Bessel function of order n and $J_n\prime (z)$ its derivative with respect to its argument; in obtaining these expressions we have made use of the Fourier decomposition of the orbit, more in Appendix A

For orbits with small eccentricity that are described by the ELL1 model, at leading order in e → 0 the osculating orbit equations become

$$\begin{array}{cc}\hfill \frac{\dot{x}}{x}& =-2\alpha \dot{\mathrm{\Phi }},\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill \frac{\dot{\eta }}{\eta }=\frac{\dot{\kappa }}{\kappa }& =-\frac{\alpha }{e}[\mathrm{\Phi }{\omega }_{b}sin{\omega }_b\tilde{t}+2\dot{\mathrm{\Phi }}cos{\omega }_b\tilde{t}]\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill \dot{\mathrm{\Psi }\prime }-{\omega }_b& =2\alpha \mathrm{\Phi }{\omega }_b.\hfill \end{array}$$(35)

Proceeding as with δΘ′, we find

$$\begin{array}{cc}\hfill \delta \mathrm{\Psi }\prime & =\delta {\mathrm{\Psi }}_{\mathrm{asc}}\prime \hfill \\ \hfill & +{\int }_{T_{\mathrm{asc}}}^t[\sqrt{\frac{G{M}_T(1+\alpha \mathrm{\Phi })}{{(a+\delta a)}^3}}-\sqrt{\frac{G{M}_T}{a^3}}]\mathrm{d}t\prime \hfill \\ \hfill & +2\alpha {\omega }_b{\int }_{T_{\mathrm{asc}}}^t\mathrm{\Phi }\mathrm{d}t\prime .\hfill \end{array}$$(36)

With Φ from Eq. (18), upon integration we find the expressions for the variations (denoted by δ) of the orbital parameters, which we collect in Appendix C

3. A two-step approach to estimate the sensitivity

Step 1: variances. In order to estimate the sensitivity to ULDM we follow the procedure outlined in Blandford & Teukolsky (1976), which we adapt to our needs by incorporating the Bayesian analysis of Moore et al. (2015). From Eq. (7) we see that the time-dependent function N is parametrised by N₀, ν, $\dot{\nu }$, $\ddot{\nu }$ and the orbital parameters, namely N = N(t, N₀, ν, …). The first step is to assume that, for each binary system, we have a reasonable ‘first guess’ as to what the values of the orbital parameters are, which we denote {N₀⁽¹⁾, ν⁽¹⁾, …}. This can be for instance obtained by fitting time-of-arrival data with the help of tools such as TEMPO2 Hobbs et al. (2006) or PINT Luo et al. (2021). We can then define the time residuals R(t) as

$$\begin{array}{c}\hfill -{\nu }^{(1)}R(t)\equiv N(t,N_0,\nu ,\dots )-N(t,N_0^{(1)},{\nu }^{(1)},\dots ),\end{array}$$(37)

where N(t, N₀, ν, …) corresponds to the observed time series (or a synthetic one with model parameters $\{N_0,\nu ,\dots \}$). Explicitly, the function R(t) reads

$$\begin{array}{cc}\hfill R(t)& =\delta K-\frac{\partial N}{\partial a}{|}_1\frac{\delta a}{\nu }\hfill \\ \hfill & -\frac{\partial N}{\partial e}{|}_1\frac{\delta e}{\nu }-\frac{\partial N}{\partial {\alpha }_b}{|}_1\frac{\delta {\alpha }_b}{\nu }\hfill \\ \hfill & -\frac{\partial N}{\partial {\eta }_b}{|}_1\frac{\delta {\eta }_b}{\nu }-\frac{\partial N}{\partial \mathrm{\Theta }\prime }{|}_1\frac{\delta \mathrm{\Theta }\prime }{\nu },\text{(BT)}\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill R(t)& =\delta K-\frac{\partial N}{\partial x}{|}_1\frac{\delta x}{\nu }\hfill \\ \hfill & -\frac{\partial N}{\partial \eta }{|}_1\frac{\delta \eta }{\nu }-\frac{\partial N}{\partial \kappa }{|}_1\frac{\delta \kappa }{\nu }\hfill \\ \hfill & -\frac{\partial N}{\partial s}{|}_1\frac{\delta s}{\nu }-\frac{\partial N}{\partial \mathrm{\Psi }\prime }{|}_1\frac{\delta \mathrm{\Psi }\prime }{\nu },\text{(ELL1)}\hfill \end{array}$$(39)

where

$$\begin{array}{c}\hfill \delta K\equiv -\frac{\delta N_0}{\nu }-\frac{\delta \nu }{\nu }t-\frac{\delta \dot{\nu }}{2\nu }{t}^2-\frac{\delta \ddot{\nu }}{6\nu }{t}^3.\end{array}$$(40)

In order to simplify the notation below we define a vector $\delta S$ containing the model parameters of the residuals that can in principle be measured for both the BT and ELL1 timing models, namely $\delta S_{\mathrm{BT}}\equiv \{\delta K,\delta a,\delta e,\delta {\alpha }_b,\delta {\eta }_b,\delta \mathrm{\Theta }\prime \}$ and $\delta S_{\mathrm{ELL}1}\equiv \{\delta K,\delta s,\delta x,\delta \eta ,\delta \kappa ,\delta \mathrm{\Psi }\prime \}$, respectively. In what follows we neglect the contribution of the term with $\dot{\nu }$ everywhere in Eq. (38) except in δK

We now assume that the observations are made at regular time intervals at a rate n_c for a duration d such that P_b ≪ d ≪ T_obs where T_obs is the total time of observation. We also assume that each observation has a constant variance ϵ² and zero covariance. In each interval of length d we minimise the χ² for the n_cd observations made in that interval and obtain an equation of the form

$$\begin{array}{c}\hfill \frac{1}{{ϵ}^2}\sum _{i=1}^{n_{c}d}{\mathsf{M}}^i\delta \mathsf{S}=\mathsf{D},\end{array}$$(41)

where the summation is over the number of observations, and we defined the vector D which contains the data. We then approximate the summation with its average times the number of observations as

$$\begin{array}{c}\hfill \sum _{i=1}^{n_{c}d}{\mathsf{M}}^i\simeq \frac{n_{c}d}{2\pi }{\int }_0^{2\pi }{\mathsf{M}}^i\mathrm{d}\mathrm{\Theta }\prime =n_{c}d\overline{\mathsf{M}},\end{array}$$(42)

where dΘ′ should be replaced by dΨ′ in the ELL1 model. Finally, the covariance matrix is obtained as

$$\begin{array}{c}\hfill \mathsf{C}=\frac{{ϵ}^2}{n_{c}d}{\overline{\mathsf{M}}}^{-1},\end{array}$$(43)

from where we can read off the variances for all the parameters that define the system–we collect these variances in Appendix B

Step 2: time-dependence. The second step consists of fitting the variances obtained in the first step to a time-dependent function which captures the evolution of the perturbations to the binary system. In Blandford & Teukolsky (1976) the time-evolution of the orbital parameters is assumed to be a polynomial up to the quadratic order. However, since we are seeking to detect the sinusoidal perturbation exerted by ULDM, we employ a Bayesian approach that mirrors the one for gravitational waves Moore et al. (2015)

We assume that the number of measurements 𝒩 = T_obs/d of the parameter set $\delta \mathsf{S}$ is large and that they are regularly made within the observation time. For any of the parameters S for which a variance has been estimated in step 1, σ_S² = var(δS) where $\{\delta S(t)\}\equiv \delta \mathsf{S}$, we define

$$\begin{array}{cc}\hfill \delta S(t)& \equiv \mathbf{m}(\mathbf{\Xi },t)+\mathbf{h}(t)+\mathbf{n},\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill \mathbf{h}(t)& \equiv \alpha [{\mathbf{A}}_X(t)X+{\mathbf{A}}_Y(t)Y],\hfill \end{array}$$(45)

$$\begin{array}{cc}\hfill \mathbf{m}(\mathbf{\Xi },t)& \equiv {\mathbf{\Xi }}^T·\mathit{N}\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill \mathrm{with}{\mathbf{\Xi }}^T& \equiv \{{\mathrm{\Xi }}_0,{\mathrm{\Xi }}_1,{\mathrm{\Xi }}_2\},{\mathit{N}}^T\equiv \{1,t,t^2\};\hfill \end{array}$$(47)

here m and h are 𝒩-dimensional vectors that describe the model–we note that these vectors are not of the same dimensionality as the parameter sets $\delta \mathsf{S}$. The parameters Ξ₀, Ξ₁ and Ξ₂ are the constants that define the quadratic fit of Blandford & Teukolsky (1976) (we note that some of the variables have Ξ₂ = 0); n is the noise of the measurement, which is assumed to be white, Gaussian and uncorrelated between pulsars; $\mathbf{h}$ is the signal caused by the presence of ULDM. For future convenience, here we have also separated the sinusoidal ULDM contribution into its $A_X(t)$ and $A_Y(t)$ contributions where $X\equiv \sqrt{2}\varrho cos\mathrm{{\rm Y}}$ and $Y\equiv \sqrt{2}\varrho sin\mathrm{{\rm Y}}$⁸. For instance, in the ELL1 model, from Eq. (33) we find

$$\begin{array}{cc}\hfill \delta x(t)& =\delta x_{\mathrm{asc}}+\dot{x}(t-T_{\mathrm{asc}})\hfill \\ \hfill & -2\alpha x{\mathrm{\Phi }}_0\varrho [cos(mt+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]\hfill \\ \hfill & ={\mathrm{\Xi }}_0^x+{\mathrm{\Xi }}_1^{x}t+h^x(t),\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill h^x(t)& =-2\alpha x{\mathrm{\Phi }}_0\varrho [cos(mt+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})],\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill A_X(t)& =-\sqrt{2}x{\mathrm{\Phi }}_0[cos(mt)-cos(m{T}_{\mathrm{asc}})],\hfill \\ \hfill A_Y(t)& =\sqrt{2}x{\mathrm{\Phi }}_0[sin(mt)-sin(m{T}_{\mathrm{asc}})],\hfill \end{array}$$(50)

where ${\mathrm{\Xi }}_2^x=0$ in this case as there is no quadratic term–further explicit expressions for the variations $\delta S$ can be found in Appendix C

The time-model for a given orbital parameter variation needs to be tested against an alternative model wherein there is no ULDM, namely $\mathbf{h}(t)=0$. We call these two competing models the signal and noise hypotheses, respectively. In order to account for our ignorance in the quadratic fitting parameters as well as the specific realisation of the local ULDM at the location of the pulsar, we marginalise over Ξ_j (which include possible secular variations of the orbital parameters) assuming uniform prior. Moreover, unless otherwise specified, we also use a uniform prior for the ULDM phase Υ, whereas for 𝜚 we use Eq. (20); thence, the probability distribution functions for the X and Y variables is given by

$$\begin{array}{c}\hfill P(X,Y)=\frac{1}{2\pi }{e}^{-\frac{X^2+Y^2}{2}}.\end{array}$$(51)

Lastly, the noise distribution can be written as

$$\begin{array}{c}\hfill P(\mathbf{n})d\mathbf{n}=\frac{exp(-\frac{1}{2}{\mathbf{n}}^T{{\mathrm{\Sigma }}_{\mathbf{n}}}^{-1}\mathbf{n})d\mathbf{n}}{\sqrt{{(2\pi )}^{\mathcal{N}}det({\mathrm{\Sigma }}_n)}},\end{array}$$(52)

where Σ_n is a diagonal matrix with all diagonal elements identical to σ_S², calculated in the first step

To compare the two hypotheses, the signal $\mathbf{h}(t)\ne 0$ and the noise $\mathbf{h}(t)=0$ for a given parameter S, we compute the ratio of the evidences 𝒪_h and 𝒪_n for the signal and the noise hypotheses, respectively. The evidences are defined as the likelihood marginalised over all the parameters of the model. The ratio of the evidences of the two models defines the Bayes factor ℬ as:

$$\begin{array}{c}\hfill \mathcal{B}\equiv \frac{{\mathcal{O}}_h}{{\mathcal{O}}_n}.\end{array}$$(53)

Explicitly, log-likelihoods for the signal and noise hypotheses are

$$\begin{array}{c}\hfill log\mathcal{L}(\mathbf{\Xi },X,Y,\alpha ,m)\propto -\frac{1}{2}{(\mathbf{d}-\mathbf{m}-\mathbf{h})}^T{{\mathrm{\Sigma }}_{\mathbf{n}}}^{-1}(\mathbf{d}-\mathbf{m}-\mathbf{h}),\end{array}$$(54)

$$\begin{array}{c}\hfill log\mathcal{L}(\mathbf{\Xi })\propto -\frac{1}{2}{(\mathbf{d}-\mathbf{m})}^T{{\mathrm{\Sigma }}_{\mathbf{n}}}^{-1}(\mathbf{d}-\mathbf{m}),\end{array}$$(55)

where d is the 𝒩−dimensional vector that contains the data. Following Moore et al. (2015) we further define the matrices

$$\begin{array}{c}\hfill \mathbf{m}\equiv \mathbf{M}\mathbf{\Xi },\mathbf{M}\equiv \mathbf{U}\mathbf{T}{\mathbf{V}}^{†},\mathbf{U}\equiv (\mathbf{F},\mathbf{G}),\end{array}$$(56)

where M is the so-called design matrix which admits a singular value decomposition into the matrices U (of dimension 𝒩 × 𝒩), T (of dimension 𝒩 × 3) and V (of dimension 3 × 3)–U is conveniently written in terms of F and G with G an 𝒩 × 𝒩 − 3 matrix⁹. From this decomposition the evidence for the noise hypothesis results in

$$\begin{array}{cc}\hfill {\mathcal{O}}_n& =\int \mathcal{L}(\mathrm{\Xi })\mathrm{d}\mathbf{\Xi }\hfill \\ \hfill & =\frac{exp(-\frac{1}{2}{\mathbf{d}}^T\mathbf{G}{({\mathbf{G}}^T{\mathrm{\Sigma }}_{\mathbf{n}}\mathbf{G})}^{-1}{\mathbf{G}}^T\mathbf{d})}{\sqrt{{(2\pi )}^{\mathcal{N}-3}det({\mathbf{G}}^T{\mathrm{\Sigma }}_n\mathbf{G})}}.\hfill \end{array}$$(57)

In order to obtain the evidence for the signal hypothesis we need to assume priors for α and the mass m. If we adopt a stringent detection threshold (namely, a Bayes factor of 1000) we can assume that the posterior distributions for α and m are sufficiently peaked at the true values, and that the choice of any reasonable prior for them would not affect this. Therefore, since we are interested in evaluating the sensitivity curve for α as a function of m, we use a delta-function prior for α and m centred at their fiducial values α_f and m_f

$$\begin{array}{c}\hfill P(\alpha ,m)=\delta (\alpha -{\alpha }_f)\delta (m-m_f).\end{array}$$(58)

Hence, the evidence for the signal hypothesis is given by

$$\begin{array}{cc}\hfill {\mathcal{O}}_h& =\int \mathrm{d}\mathbf{\Xi }\mathrm{d}X\mathrm{d}Y\mathrm{d}\alpha \mathrm{d}mP(X,Y)P(\alpha ,m)\hfill \\ \hfill & \times \mathcal{L}(\mathrm{\Xi },X,Y,\alpha ,m).\hfill \end{array}$$(59)

After performing the two trivial integrations over α and m, we are left with a number of gaussian integrals which can be computed straightforwardly, more in Appendix C

Finally, we average the Bayes factor over many realisations of the noise $\overline{\mathcal{B}}\equiv {\displaystyle \int \mathrm{d}}\mathbf{n}P(\mathbf{n})\mathcal{B}$. The result is

$$\begin{array}{cc}\hfill \overline{\mathcal{B}}& =exp\left\{{\alpha }_f^2\frac{(u_X^2+u_Y^2)}{2}\right\},\hfill \end{array}$$(60)

$$\begin{array}{cc}\hfill u_X& \equiv {\mathbf{A}}_X^T\mathbf{G}{({\mathbf{G}}^T{\mathrm{\Sigma }}_{\mathbf{n}}\mathbf{G})}^{-1}{\mathbf{G}}^T\mathbf{h},\hfill \end{array}$$(61)

$$\begin{array}{cc}\hfill u_Y& \equiv {\mathbf{A}}_Y^T\mathbf{G}{({\mathbf{G}}^T{\mathrm{\Sigma }}_{\mathbf{n}}\mathbf{G})}^{-1}{\mathbf{G}}^T\mathbf{h},\hfill \end{array}$$(62)

where in these quantities (their definition is in Appendix C) all the parameters are replaced by their fiducial values (i.e. α → α_f, m → m_f, 𝜚 → 𝜚_f, and Υ → Υ_f). If instead of using P(X, Y) given in Eq. (51) we use delta-function priors also for X and Y (which is equivalent to assuming that we know the exact ULDM configuration for a given pulsar), we find $\overline{\mathcal{B}}=exp\left\{u\right\}$ where $u\equiv {\mathbf{h}}^T\mathbf{G}{({\mathbf{G}}^T{\mathrm{\Sigma }}_{\mathbf{n}}\mathbf{G})}^{-1}{\mathbf{G}}^T\mathbf{h}$

At last, we can combine all binary pulsars by multiplying the Bayes factors as

$$\begin{array}{c}\hfill {\overline{\mathcal{B}}}^C=\prod _{p=1}^{N_p}{\overline{\mathcal{B}}}_p,\end{array}$$(63)

where N_p is the total number of binary pulsars, ${\overline{\mathcal{B}}}_p$ corresponds to the p-th binary system for which the Bayes factors individually are given in Eq. (60). The constraint ${\alpha }_f(m_f)$, which we call the sensitivity, is obtained by choosing a threshold value for detection, customarily set to be

$$\begin{array}{c}\hfill {\overline{\mathcal{B}}}^C={\overline{\mathcal{B}}}_{\mathrm{thr}}^C=1000,\end{array}$$(64)

readers can refer to Moore et al. (2015). For simplicity, in what follows, and in particular in Sects. 4 and 5, we drop the subscript ‘f’ for the fiducial values of α and m

4. Results

4.1. Near-circular orbits

As a first application of our method, we apply it to a set of 19 binary pulsars for which the orbital parameters were obtained by the NANOGrav collaboration Gabriella et al. (2023) using the ELL1 model and 3 binary pulsars modelled with the ELL1H (a variant on ELL1 which differs only in the parametrisation of the relativistic Shapiro time delay which we do not consider here Hobbs et al. (2006), Luo et al. (2021)) – we list the properties of all the binary pulsars that we use in this work in Tables E.1 and E.2. For each binary system we use the NANOGrav values of T_obs (MJD), ϵ and n_c, and assume ${\rho }_{\mathrm{DM}}=0.3\mathrm{GeV}/{\mathrm{cm}}^3$. We note that, as we are forecasting the constraining power of data, and not running an actual data analysis, we need to assume the fiducial values for ϱ and Υ, or, equivalently, X and Y; in other words, we need to simulate our data vectors and with them the unknown values of the ULDM parameters; we do this by assigning separately to each binary system a random realisation of such parameters as taken from their corresponding theoretical distributions

In Fig. 2 we show the sensitivity obtained by applying Eq. (64) to the parameters $\delta x$ and $\delta \mathrm{\Psi }\prime$ for ten different realisations of the fiducial ULDM model. In this figure we also show the current bounds obtained by Solar System constraints Bertotti et al. (2003) and Armstrong et al. (2003), and the bounds obtained with pulsar-timing arrays (PTAs) Porayko et al. (2018), Smarra et al. (2023). We see that, first of all, for a wide range of ULDM masses the sensitivity of binary pulsars supersedes the current Solar System constraints. This is mostly coming from the $\delta \mathrm{\Psi }\prime$ orbital parameter, which is by far more sensitive for small masses. The $\delta x$ parameter, closely related to the semi-major axis a, is the more sensitive of the two at large ULDM masses, but it currently does not give a competitive constraint. However, current PTA searches provide the most stringent limits on α in the low-frequency regime, which are two to three orders of magnitude better than our limits

Thanks to the near-analyticity of our model we can readily predict how much this result would improve for any given set of parameters of future observations. In Fig. 3 we combine 111 systems taken from the ATNF pulsar catalogue Manchester et al. (2005) and, based on the analysis presented in Liu et al. (2011) for a new generation of radio telescopes (e.g. the Square Kilometre Array, SKA), as an illustration of forecasting for such experiments we assume T_obs = 10 years, ϵ = 0.1 μs, n_c = 1/ day. As we can see in the figure, the sensitivity in this case improves by well over an order of magnitude, bringing the binary pulsar constraint below the Solar System one for a larger range of ULDM masses. In particular, the sensitivity obtained with our method extends to masses beyond the Nyquist frequency of PTAs, and can beat Solar System tests in that range by up to one order of magnitude

In order to show how our results depend on the choice of priors for the ULDM parameters X and Y (corresponding to choices of ϱ and Υ), in Figs. 2 and 3 we also show the sensitivity curves for delta-function priors $\delta (X_{\mathrm{fid}}-X)$ and $\delta (Y_{\mathrm{fid}}-Y)$. We observe that, while the choice of prior is not entirely negligible, nonetheless the effect is limited to a factor of order unity. Lastly, we note that the variations of the Laplace-Lagrange orbital parameters η and κ do exhibit resonances (as indicated in Eqs. (C.20) and (C.21)); however, beyond the resonant peak these parameters are less sensitive than x and $\mathrm{\Psi }\prime$

4.2. Resonances

When a system possesses a significant eccentricity e its orbital parameters experience resonant secular drifts Blas et al. (2017, 2020), and obtaining the sensitivity becomes much more computationally intensive. Therefore, here we limit ourselves to the effects of ULDM on P_b, and not the full variation of the parameter $\mathrm{\Theta }\prime$ defined in Eqs. (9) and (30) and refer to Appendix D for a more detailed discussion on this simplification, which does not appreciably alter our results. We call this contribution δΘ′_Pb, defined as

$$\begin{array}{cc}\hfill \delta \mathrm{\Theta }{\prime }_{P_b}& =\delta \mathrm{\Theta }{\prime }_0-{\int }_{T_0}^t{\omega }_b\frac{\delta P_b}{P_b}\mathrm{d}t\prime .\hfill \end{array}$$(65)

With this proviso, in Figs. 4–6 we show the sensitivity obtained from this simplified $\delta \mathrm{\Theta }\prime$ for the three binary pulsars J1946+3417, J2234+0611 and J1903+0327, respectively (we collect their parameters in Table E.2); as in Sect. 4.1 we run ten realisations of the ULDM parameters for each binary system

The combined sensitivity from all three systems, which beats each of them separately, as expected, is shown in Fig. 7. In these figures we limit ourselves to a smaller mass range to focus on the resonances, and we can see how these systems are in effect less constraining than the near-circular ones, their combined sensitivity approaching that of Solar System tests only for the most peaked resonances at low masses. Finally, the results of this section are directly comparable to Blas et al. (2020); we see that, while we reproduce these results very closely near the resonances, the off-resonance sensitivity as obtained from our method is in fact better than originally estimated

5. Conclusion and outlook

Summary. In this paper we computed the sensitivity of binary pulsars for the detection of scalar ULDM for the entire mass range ${10}^{-23}\mathrm{eV}\lesssim m\lesssim {10}^{-18}\mathrm{eV}$, and for both eccentric and near-circular systems. We improve on previous results which only considered resonances between the systems’ orbital periods and ULDM oscillations by developing a method that goes beyond resonances. With the new method it is also possible to quantify our ignorance on the local DM field (i.e. the local amplitude determined by the parameter 𝜚 and the phase Υ). We find that, while our method is not competitive in the range of frequencies surveyed by PTAs, beyond the Nyquist frequency binary pulsars can be the most sensitive detectors for ULDM. For instance, we expect that observations made with the future SKA would be sensitive to ULDM couplings $\alpha \approx \mathrm{few}\times {10}^{-22}{\mathrm{GeV}}^{-1}$ for $m\approx \mathrm{few}\times {10}^{-22}\mathrm{eV}$, better than Solar System constraints. In this regard, we notice that the most sensitive orbital parameter is $\mathrm{\Psi }\prime$ (for ELL1 binaries) or $\mathrm{\Theta }\prime$ (for BT binaries), highlighting the relevance of orbital parameters other than the most-studied (and often best-measured) orbital period $P_b(t)$

Relevance. Our results highlight why it is important to go beyond resonances. Firstly, resonances, being narrow resonances, are limited to a very narrow range of ULDM masses in correspondence with integer multiples of the binary system inverse period: going beyond resonances allows us to extend the sensitivity of each binary system to the entire mass range of interest, only limited by the binary pulsar back-reactions on the field Φ on the higher end Blas et al. (2020) and cosmological constraints on ${\rho }_{\mathrm{DM}}$ on the lower end (readers can refer, for instance, to Shevchuk et al. (2023) and references therein). Secondly, in the scalar ULDM case resonances disappear for near-circular orbits; since near-circular orbits constitute the majority of binary pulsars in nature, for scalar ULDM searches we are limited to a minority of systems. Going beyond resonances means that each observed binary system contributes to the sensitivity for ULDM. Thirdly, each binary system contributes separately to the sensitivity to ULDM, as obtained from the secular orbital drifts caused by resonances. With our new method all binary pulsars jointly contribute to the sensitivity, thereby improving it across all masses

Context. Our findings should be compared with existing limits on scalar ULDM. We have already shown the limit coming from the Doppler tracking of the Cassini spacecraft, which results in $\alpha \lesssim {10}^{-21}{\mathrm{GeV}}^{-1}$ across all the masses we consider here Bertotti et al. (2003). Because a putative ULDM direct coupling α is equivalent to a scalar gravitational wave, the Cassini spacecraft also sets a constrain on the mass range $5\times {10}^{-21}\mathrm{eV}\lesssim m\lesssim {10}^{-18}\mathrm{eV}$ that is about an order of magnitude stronger Armstrong et al. (2003). In the lowest frequency range considered in this work, current searches of the stochastic gravitational wave background by pulsar-timing arrays set the limits on ULDM around $\alpha \lesssim {10}^{-26}{\mathrm{GeV}}^{-1}$ for a decade in mass around ${10}^{-23}\mathrm{eV}\lesssim m\lesssim {10}^{-22}\mathrm{eV}$Porayko et al. (2018), Smarra et al. (2023). Furthermore the dynamical friction within the solitonic core of the ULDM halo changes the spectral slope of the stochastic background of gravitational waves which is at odds with what has been measured by pulsar-timing arrays, possibly implying the exclusion of ULDM masses within the range $1.3\times {10}^{-21}\mathrm{eV}\lesssim m\lesssim 1.4\times {10}^{-20}\mathrm{eV}$Aghaie et al. (2023). Existing limits also come from superradiance. By measuring the spin and mass of known black holes and other astrophysical objects, and assuming no interaction between Φ and standard matter, the mass ranges $2.5\times {10}^{-21}\mathrm{eV}\lesssim m\lesssim 1.2\times {10}^{-20}\mathrm{eV}$ and $5.5\times {10}^{-20}\mathrm{eV}\lesssim m\lesssim 1.3\times {10}^{-16}\mathrm{eV}$ are excluded, or else these black holes and celestial objects would not be there Stott (2020). Therefore, the limits we obtain from binary pulsars at the same time constitute an independent test of some of the parameter space that is probed by superradiance, as well as a test of regions not accessible by it

Prospects. In this work we have concerned ourselves with the simplest case of a universally-coupled linear scalar ULDM field only. In the future we plan to investigate the spin-1 and spin-2 ULDM cases, for which only the resonant effects have been studied López Nacir & Urban (2018), Armaleo et al. (2020a). The main difference with respect to the scalar case is that, if ULDM has spin 1 or spin 2, even circular orbits experience resonant secular orbital drifts. Moreover, since the geometry of the ULDM perturbations on the orbital motion is different, we expect that it is possible to differentiate between different ULDM spins. The spin-2 case is especially interesting because ULDM with spin 2 mimics a continuous gravitational wave, albeit with five polarisations, readers can refer to Armaleo et al. (2021). Important extension to this work are the cases of a non-universal coupling as well as a quadratic coupling, which present distinctive phenomenology compared to the linear case–readers can refer, for instance, to Blas et al. (2020). Furthermore, a more detailed analysis would include other sources of noise (such as pulsar red noise) as well as more complex and complete timing models, which are optimally tailored to each binary pulsar system. Finally, we plan to empirically test our results by directly injecting a model ULDM signal to simulated binary pulsar time-of-arrival data, and then using different methods (applied to simulated data) in order to recover the injected ULDM signal within the predicted level of sensitivity

With the advent of the next generation of radio surveys such as the SKA we expect to observe an order-of-magnitude more binary pulsars than today. With our method it is immediate to combine all the systems, each with very different properties, into a single sensitivity curve across all relevant masses, and thus assess the chances for either ULDM detection, or the tightest constraints on its parameter space as new systems are discovered

---

Acknowledgments

The Authors wish to thank Stéphane Ilić for useful discussion. DLN acknowledges support from ESIF and MEYS (Project ‘FZU researchers, technical and administrative staff mobility’ – CZ.02.2.69/0.0/0.0/18_053/0016627), UBA and CONICET. PK and FU acknowledge support from MEYS through the INTER-EXCELLENCE II, INTER-COST grant LUC23115. This article is based upon work from COST Action COSMIC WISPers CA21106, supported by COST (European Cooperation in Science and Technology)

References

Appendix A: Post-Keplerian orbits

The orbital parameters of a perturbed binary system evolve with time according to the equations Turner (1979), Danby (1970),

$$\begin{array}{cc}\hfill \dot{a}& =\frac{2}{{\omega }_b}\{\frac{F_{r}e}{\sqrt{1-e^2}}sin\theta +\frac{F_{\theta }}{r}a\sqrt{1-e^2}\},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \dot{e}& =\frac{\sqrt{1-e^2}}{a{\omega }_b}\{F_{r}sin\theta +F_{\theta }(cos\theta +\frac{1}{e}-\frac{r}{\mathit{ae}})\},\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill \dot{\mathrm{\Omega }}& =\frac{F_{z}rsin(\theta +\omega )}{a^2{\omega }_b\sqrt{1-e^2}sin\iota },\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill \dot{\iota }& =\frac{F_{z}rcos(\theta +\omega )}{a^2{\omega }_b\sqrt{1-e^2}},\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill \dot{\varpi }& =\frac{\sqrt{1-e^2}}{ae{\omega }_b}\{-F_{r}cos\theta +F_{\theta }sin\theta [1+\frac{r}{a(1-e^2)}]\}+2{sin}^2\left(\frac{\iota }{2}\right)\dot{\mathrm{\Omega }},\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\dot{ϵ}}_1& =-\frac{2r{F}_r}{a^2{\omega }_b}+[1-\sqrt{1-e^2}]\dot{\varpi }+2\sqrt{1-e^2}{sin}^2\left(\frac{\iota }{2}\right)\dot{\mathrm{\Omega }},\hfill \end{array}$$(A.6)

where

$$\begin{array}{c}\hfill {\omega }_b=\sqrt{\frac{G{M}_T}{a^3}},\end{array}$$(A.7)

is the binary orbital frequency and we have decomposed the perturbing force in the cylindric coordinates,

$$\begin{array}{c}\hfill \mathit{F}=F_r\widehat{r}+F_{\theta }\widehat{\theta }+F_z\widehat{z}.\end{array}$$(A.8)

In Eqs. (A.5) and (A.6) we have introduced the argument of the pericentre,

$$\begin{array}{c}\hfill \varpi =\omega +\mathrm{\Omega },\end{array}$$(A.9)

and the parameter

$$\begin{array}{c}\hfill {ϵ}_1={\omega }_b(t-T_0)+\varpi -{\int }_{T_0}^t\mathrm{d}{t}^{\prime }{\omega }_b.\end{array}$$(A.10)

The orbital motion for elliptic orbits is not expressible in closed form, but can be represented as Fourier decomposition of combinations of the unperturbed Keplerian orbits, which amounts to a harmonic expansion in eccentricities Watson (1995):

$$\begin{array}{cc}& \frac{r}{a}=1+\frac{e^2}{2}-2e\sum _{n=1}^{\infty }\frac{J_n^{\prime }(ne)}{n}cosn{\omega }_b\tilde{t},\hfill \end{array}$$(A.11)

$$\begin{array}{cc}& \frac{r^2}{a^2}=1+\frac{3}{2}{e}^2-4\sum _{n=1}^{\infty }\frac{J_n(ne)}{n^2}cosn{\omega }_b\tilde{t},\hfill \end{array}$$(A.12)

$$\begin{array}{cc}& cos\theta =-e+\frac{2(1-e^2)}{e}\sum _{n=1}^{\infty }{J}_n(ne)cosn{\omega }_b\tilde{t},\hfill \end{array}$$(A.13)

$$\begin{array}{cc}& sin\theta =2\sqrt{1-e^2}\sum _{n=1}^{\infty }{J}_n^{\prime }(ne)sinn{\omega }_b\tilde{t},\hfill \end{array}$$(A.14)

$$\begin{array}{cc}& \frac{r}{a}cos\theta =-\frac{3e}{2}+2\sum _{n=1}^{\infty }\frac{J_n^{\prime }(ne)}{n}cosn{\omega }_b\tilde{t},\hfill \end{array}$$(A.15)

$$\begin{array}{cc}& \frac{r}{a}sin\theta =\frac{2\sqrt{1-e^2}}{e}\sum _{n=1}^{\infty }\frac{J_n(ne)}{n}sinn{\omega }_b\tilde{t},\hfill \end{array}$$(A.16)

$$\begin{array}{cc}& \frac{a^2}{r^2}cos\theta =\sum _{n=1}^{\infty }2n{J}_n^{\prime }(ne)cosn{\omega }_b\tilde{t},\hfill \end{array}$$(A.17)

$$\begin{array}{cc}& \frac{a^2}{r^2}sin\theta =\frac{\sqrt{1-e^2}}{e}\sum _{n=1}^{\infty }2n{J}_n(ne)sinn{\omega }_b\tilde{t},\hfill \end{array}$$(A.18)

where $\stackrel{\sim }{t}=t-T_0$, J_n(z) is a Bessel function of order n and J_n′(z) is its derivative with respect to z

Appendix B: Variances

B.1. Elliptic orbits

The time residual for the BT model is given by Eq. (38) in terms of six independent parameters: $\{\delta K,\delta a,\delta e,\delta {\alpha }_b,\delta {\eta }_b,\delta {\mathrm{\Theta }}^{\prime }\}$. The dependence on δa only comes from the dependence on ω_b in the last term of Eq. (7), which is a subdominant contribution with respect to the other terms. Here, for simplicity, we follow Blandford & Teukolsky (1976) and we neglect such contribution to estimate the variances of the remaining five parameters: $\{\delta K,\delta e,\delta {\alpha }_b,\delta {\eta }_b,\delta {\mathrm{\Theta }}^{\prime }\}$. The dominant part of the time residual can be written as

$$\begin{array}{cc}\hfill R^{\mathrm{BT}}& =\delta {\eta }_{b}sin{E}^{\prime }-\delta e[{\alpha }_b+\frac{sin{E}^{\prime }({\alpha }_{b}sin{E}^{\prime }-{\eta }_{b}cos{E}^{\prime })}{1-ecos{E}^{\prime }}]+\delta {\alpha }_b(cos{E}^{\prime }-e)\hfill \\ \hfill & -\frac{\delta {\mathrm{\Theta }}^{\prime }({\alpha }_{b}sin{E}^{\prime }-{\eta }_{b}cos{E}^{\prime })}{1-ecos{E}^{\prime }},\hfill \end{array}$$(B.1)

and the variances for the parameters for generic values of e, α_b and η_b are as follows

$$\begin{array}{cc}\hfill \mathrm{var}(\delta K)& =-\mathcal{Q}\{{\eta }_b^4[e^4(26-9\epsilon )+4e^2(8\epsilon -7)+8(\epsilon -1)]\hfill \\ \hfill & +2{\alpha }_b^2{\eta }_b^2[e^6+e^4(18-5\epsilon )+4e^2(6\epsilon -5)+8(\epsilon -1)]\hfill \\ \hfill & +{\alpha }_b^4[8e^8+e^6(32\epsilon -62)+e^4(74-65\epsilon )+4e^2(4\epsilon -3)+8(\epsilon -1)]\},\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta {\alpha }_b)& =\mathcal{Q}\{4{\eta }_b^4[e^2(\epsilon -3)-4\epsilon +4]+2{\alpha }_b^2{\eta }_b^2[7e^4-22e^2-2(e^4-7e^2+8)\epsilon +16]\hfill \\ \hfill & +{\alpha }_b^4[e^6(\epsilon -6)+e^4(30-16\epsilon )+8e^2(4\epsilon -5)-16(\epsilon -1)]\},\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta {\eta }_b)& =\mathcal{Q}\{4{\alpha }_b^4{(e^2-1)}^2[e^2(\epsilon -3)-4\epsilon +4]+2{\eta }_b^4[e^4+4e^2(\epsilon -2)-8\epsilon +8]\hfill \\ \hfill & +{\alpha }_b^2{\eta }_b^2[e^6(\epsilon -6)+e^4(42-20\epsilon )+e^2(52\epsilon -68)-32(\epsilon -1)]\},\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta e)& =\mathcal{Q}{e}^4\{{\alpha }_b^2(e^2-2)[e^2(\epsilon -2)-2\epsilon +2]-2{\eta }_b^2(e^2+2\epsilon -2)\},\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta {\mathrm{\Theta }}^{\prime })& =\mathcal{Q}\epsilon e^2\{2{\alpha }_b^2(e^2-1)[e^2(\epsilon -2)-2\epsilon +2]+{\eta }_b^2(e^2-2)(e^2+2\epsilon -2)\},\hfill \end{array}$$(B.6)

where

$$\begin{array}{c}\hfill \mathcal{Q}\equiv \frac{{ϵ}^2}{d{n}_c[e^4+4e^2(\epsilon -2)-8\epsilon +8]{[{\eta }_b^2-{\alpha }_b^2(e^2-1)]}^2},\end{array}$$(B.7)

and $\epsilon \equiv \sqrt{1-e^2}$

B.2. Near-circular orbits

In the ELL1 timing model it is not difficult to obtain an analytic estimate for the variance of the six orbital parameters. These expressions are however very long and not illuminating, and for simplicity here we only quote the results for η ≃ κ ≃ 0

$$\begin{array}{cc}\hfill \mathrm{var}(\delta K)=& \frac{{ϵ}^2(34x^4{\omega }_b^4+61x^2{\omega }_b^2+76)}{d{n}_c(34x^4{\omega }_b^4-11x^2{\omega }_b^2+4)},\hfill \end{array}$$(B.8)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta s)=& \frac{4s^2{ϵ}^2(5x^4{\omega }_b^4+80x^2{\omega }_b^2+32)}{81d{n}_c{x}^6{\omega }_b^4},\hfill \end{array}$$(B.9)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta x)=& \frac{20{ϵ}^2}{9d{n}_c},\hfill \end{array}$$(B.10)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta \eta )=& \frac{32{ϵ}^2(x^2{\omega }_b^2+1)}{d{n}_c{x}^2(34x^4{\omega }_b^4-11x^2{\omega }_b^2+4)},\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta \kappa )=& \frac{32{ϵ}^2}{9d{n}_c{x}^4{\omega }_b^2},\hfill \end{array}$$(B.12)

$$\begin{array}{cc}\hfill \mathrm{var}(\delta {\mathrm{\Psi }}^{\prime })=& \frac{4{ϵ}^2(17x^2{\omega }_b^2+2)}{d{n}_c{x}^2(34x^4{\omega }_b^4-11x^2{\omega }_b^2+4)}.\hfill \end{array}$$(B.13)

These expressions can be further simplified by noticing that the quantity

$$\begin{array}{c}\hfill {\omega }_b^2{x}^2=\frac{G{M}_2}{a}\frac{{sin}^2\iota }{(1+\frac{M_1}{M_2})}\end{array}$$(B.14)

is very small, so that at leading order we have:

$$\begin{array}{cc}& \mathrm{var}(\delta K)=\frac{19{ϵ}^2}{d{n}_c},\mathrm{var}(\delta s)=\frac{128s^2{ϵ}^2}{81d{n}_c{x}^6{\omega }_b^4},\mathrm{var}(\delta x)=\frac{20{ϵ}^2}{9d{n}_c},\hfill \\ \hfill & \mathrm{var}(\delta \eta )=\frac{8{ϵ}^2}{d{n}_c{x}^2},\mathrm{var}(\delta \kappa )=\frac{32{ϵ}^2}{9d{n}_c{x}^4{\omega }_b^2},\mathrm{var}(\delta {\mathrm{\Psi }}^{\prime })=\frac{2{ϵ}^2}{d{n}_c{x}^2}.\hfill \end{array}$$(B.15)

We note that the variances of δs and δκ are relatively large since powers of the small quantity xω_b appear in the denominator. This is expected because the dependence on s only enters in the factor ω_b in the last term of Eq. (14), which (as in the BT model) is a subdominant contribution. By ignoring such contribution the calculation of the variances of the remaining five independent parameters yields:

$$\begin{array}{cc}& \mathrm{var}(\delta K)=\frac{(19+18{\kappa }^2){ϵ}^2}{d{n}_c},\mathrm{var}(\delta x)=\frac{2{ϵ}^2}{d{n}_c},\hfill \\ \hfill & \mathrm{var}(\delta \eta )=\frac{2{ϵ}^2(4+{\eta }^2+4{\kappa }^2)}{d{n}_c{x}^2},\mathrm{var}(\delta \kappa )=\frac{2{ϵ}^2(4+4{\eta }^2+{\kappa }^2)}{d{n}_c{x}^2},\mathrm{var}(\delta {\mathrm{\Psi }}^{\prime })=\frac{2{ϵ}^2}{d{n}_c{x}^2}.\hfill \end{array}$$(B.16)

Therefore, except for δκ, this relatively simpler calculation gives similar results. The reason why the variance of δκ becomes smaller by neglecting the contribution of δs is that at leading (non-vanishing) order in xω_b the term with δs, which is suppressed by xω_b, becomes completely degenerate with the one with δκ (they have the same time dependence given by sin2Ψ′), which does not vanish for xω_b → 0. It is only after including the $\mathcal{O}(1)$ correction in xω_b for δκ that those parameters can be distinguished

Appendix C: Bayesian approach

C.1. Perturbed BT model

The model for the variations of the orbital parameters for eccentric orbits including the ULDM contribution is given by:

$$\begin{array}{cc}\hfill \delta {\mathrm{\Theta }}^{\prime }(t)& =\delta {\mathrm{\Theta }}_0^{\prime }+\frac{5}{2}\alpha {\omega }_b\frac{{\mathrm{\Phi }}_0\varrho }{m}[sin(mt+\mathrm{{\rm Y}})-sin(m{T}_0+\mathrm{{\rm Y}})]+\alpha {\mathrm{\Phi }}_0\varrho \sum _{n=1}^{+\infty }{B}_{n}sin(n{\omega }_b\tilde{t})cos(mt+\mathrm{{\rm Y}})\hfill \\ \hfill & +\alpha \frac{{\mathrm{\Phi }}_0\varrho }{2}\sum _{n=1}^{+\infty }(A_n-n{\omega }_b{B}_n)\{\frac{sin(mt+\mathrm{{\rm Y}}+n{\omega }_b\tilde{t})-sin(m{T}_0+\mathrm{{\rm Y}})}{m+n{\omega }_b}\hfill \\ \hfill & +\frac{sin(mt+\mathrm{{\rm Y}}-n{\omega }_b\tilde{t})-sin(m{T}_0+\mathrm{{\rm Y}})}{m-n{\omega }_b}\}-\frac{3}{2}{\omega }_b{\int }_{T_0}^t\frac{\delta a(t^{\prime })}{a}d{t}^{\prime },\hfill \end{array}$$(C.1a)

$$\begin{array}{cc}\hfill \frac{\delta a(t)}{a}& =\frac{\delta a_0}{a}+\frac{{\dot{a}}_0\tilde{t}}{a}-2\alpha {\mathrm{\Phi }}_0\varrho [cos(mt+\mathrm{{\rm Y}})-cos(m{T}_0+\mathrm{{\rm Y}})]\hfill \\ \hfill & -8\alpha {\mathrm{\Phi }}_0\varrho \sum _{n=1}^{+\infty }{J}_n(ne)[cos(mt+\mathrm{{\rm Y}})cos(n{\omega }_b\tilde{t})-cos(m{T}_0+\mathrm{{\rm Y}})]\hfill \\ \hfill & +6\alpha {\mathrm{\Phi }}_0\varrho {\omega }_b\sum _{n=1}^{+\infty }n{J}_n(ne)\{\frac{cos(mt+\mathrm{{\rm Y}}+n{\omega }_b\tilde{t})-cos(m{T}_0+\mathrm{{\rm Y}})}{m+n{\omega }_b}\hfill \\ \hfill & -\frac{cos(mt+\mathrm{{\rm Y}}-n{\omega }_b\tilde{t})-cos(m{T}_0+\mathrm{{\rm Y}})}{m-n{\omega }_b}\},\hfill \end{array}$$(C.1b)

$$\begin{array}{cc}\hfill \delta e(t)& =\delta e_0+{\dot{e}}_0\tilde{t}+\alpha {\mathrm{\Phi }}_0\varrho \frac{1-e^2}{e}\sum _{n=1}^{+\infty }{J}_n(ne)\{\frac{n{\omega }_b-2m}{n{\omega }_b+m}[cos(mt+\mathrm{{\rm Y}}+n{\omega }_b\tilde{t})-cos(m{T}_0+\mathrm{{\rm Y}})]\hfill \\ \hfill & +\frac{2m+n{\omega }_b}{n{\omega }_b-m}[cos(mt+\mathrm{{\rm Y}}-n{\omega }_b\tilde{t})-cos(m{T}_0+\mathrm{{\rm Y}})]\},\hfill \end{array}$$(C.1c)

$$\begin{array}{cc}\hfill \delta \omega (t)& =\delta {\omega }_0+{\dot{\omega }}_0\tilde{t}+\alpha {\mathrm{\Phi }}_0\varrho \frac{\sqrt{1-e^2}}{e}\sum _{n=1}^{+\infty }{J}_n^{\prime }(ne)\{\frac{n{\omega }_b-2m}{n{\omega }_b+m}[sin(mt+\mathrm{{\rm Y}}+n{\omega }_b\tilde{t})-sin(m{T}_0+\mathrm{{\rm Y}})]\hfill \\ \hfill & -\frac{2m+n{\omega }_b}{n{\omega }_b-m}[sin(mt+\mathrm{{\rm Y}}-n{\omega }_b\tilde{t})-sin(m{T}_0+\mathrm{{\rm Y}})]\},\hfill \end{array}$$(C.1d)

where $\stackrel{\sim }{t}=t-T_0$, A_n and B_n are defined in Eq. (31), $\delta a_0$, $\delta e_0$ and $\delta {\omega }_0$ account for a constant shift due to the error in the determination of the parameters evaluated at T₀, and the derivatives of the orbital parameters on the right-hand side account for the possibility of having another perturbation to the system (on top of the ULDM effect) which produces a secular variation. The variations $\delta \mathrm{\Omega }$ and $\delta s$ do not have contributions from ULDM. We note that in the two-step approach presented in the main text we marginalise over those errors for each variable; we can therefore discard ULDM-dependent terms that can be absorbed in the errors of the orbital parameters, since the ULDM contribution cannot be determined separately

For the sake of completeness we write below the ULDM contribution decomposed into its X and Y components:

$$\begin{array}{cc}\hfill \delta S^{\mathrm{DM}}(t)& =h^S(t)+{\ell }^S(t),\hfill \end{array}$$(C.2)

$$\begin{array}{ccc}\hfill h^S(t)& \equiv \hfill & \hfill \alpha [A_X^S(t)X+A_Y^S(t)Y],\end{array}$$(C.3)

$$\begin{array}{cc}\hfill A_X^S(t)& =H_X^S(t)-H_X^S(T_{\mathrm{ref}}),\hfill \end{array}$$(C.4)

$$\begin{array}{cc}\hfill A_Y^S(t)& =H_Y^S(t)-H_Y^S(T_{\mathrm{ref}}),\hfill \end{array}$$(C.5)

$$\begin{array}{ccc}\hfill {\ell }^S(t)& \equiv \hfill & \hfill \alpha [L_X^S(t)X+L_Y^S(t)Y],\end{array}$$(C.6)

where h^S(t) is the characteristic signals caused by ULDM on the variables δS that enter the likelihood in Eq. (54), for the BT model T_ref = T₀,

$$\begin{array}{cc}\hfill H_X^{{\mathrm{\Theta }}^{\prime }}(t)& =\sqrt{2}{\omega }_b\frac{{\mathrm{\Phi }}_0}{m}sinmt+\frac{{\mathrm{\Phi }}_0}{\sqrt{2}}\sum _{n=1}^{+\infty }{B}_{n}sin(n{\omega }_b\tilde{t})cosmt+\frac{{\mathrm{\Phi }}_0}{2^{3/2}}\sum _{n=1}^{+\infty }(A_n-n{\omega }_b{B}_n)\hfill \\ \hfill & \times \{\frac{sin(mt+n{\omega }_b\tilde{t})}{m+n{\omega }_b}+\frac{sin(mt-n{\omega }_b\tilde{t})}{m-n{\omega }_b}\}+H_{P_b,X}^{{\mathrm{\Theta }}^{\prime }}(t),\hfill \end{array}$$(C.7)

$$\begin{array}{cc}\hfill H_Y^{{\mathrm{\Theta }}^{\prime }}(t)& =\sqrt{2}{\omega }_b\frac{{\mathrm{\Phi }}_0}{m}cosmt-\frac{{\mathrm{\Phi }}_0}{\sqrt{2}}\sum _{n=1}^{+\infty }{B}_{n}sin(n{\omega }_b\tilde{t})sinmt+\frac{{\mathrm{\Phi }}_0}{2^{3/2}}\sum _{n=1}^{+\infty }(A_n-n{\omega }_b{B}_n)\hfill \\ \hfill & \times \{\frac{cos(mt+n{\omega }_b\tilde{t})}{m+n{\omega }_b}+\frac{cos(mt-n{\omega }_b\tilde{t})}{m-n{\omega }_b}\}+H_{P_b,Y}^{{\mathrm{\Theta }}^{\prime }}(t),\hfill \end{array}$$(C.8)

$$\begin{array}{cc}\hfill H_{P_b,X}^{{\mathrm{\Theta }}^{\prime }}(t)& =-\frac{{\mathrm{\Phi }}_0{\omega }_b}{2^{3/2}}\{\frac{-7sinmt}{m}+\sum _{n=1}^{+\infty }\frac{6J_n(ne)}{(n{\omega }_b-m)(m+n{\omega }_b)}\hfill \\ \hfill & \times [\frac{(m-n{\omega }_b)(2m-n{\omega }_b)sin(mt+tn{\omega }_b-T_{0}n{\omega }_b)}{m+n{\omega }_b}\hfill \\ \hfill & +\frac{(m+n{\omega }_b)(2m+n{\omega }_b)sin(mt-tn{\omega }_b+T_{0}n{\omega }_b)}{m-n{\omega }_b}]\},\hfill \end{array}$$(C.8)

$$\begin{array}{cc}\hfill H_{P_b,Y}^{{\mathrm{\Theta }}^{\prime }}(t)& =-\frac{{\mathrm{\Phi }}_0{\omega }_b}{2^{3/2}}\{\frac{-7cosmt}{m}+\sum _{n=1}^{+\infty }\frac{6J_n(ne)}{(n{\omega }_b-m)(m+n{\omega }_b)}\hfill \\ \hfill & \times [\frac{(m-n{\omega }_b)(2m-n{\omega }_b)cos(mt+tn{\omega }_b-T_{0}n{\omega }_b)}{m+n{\omega }_b}\hfill \\ \hfill & +\frac{(m+n{\omega }_b)(2m+n{\omega }_b)cos(mt-tn{\omega }_b+T_{0}n{\omega }_b)}{m-n{\omega }_b}]\},\hfill \end{array}$$(C.10)

$$\begin{array}{cc}\hfill H_X^a(t)& =-\sqrt{2}{\mathrm{\Phi }}_{0}acosmt-4\sqrt{2}{\mathrm{\Phi }}_{0}a\sum _{n=1}^{+\infty }{J}_n(ne)cosmtcos(n{\omega }_b\tilde{t})\hfill \\ \hfill & +3\sqrt{2}{\mathrm{\Phi }}_{0}a{\omega }_b\sum _{n=1}^{+\infty }n{J}_n(ne)\{\frac{cos(mt+n{\omega }_b\tilde{t})}{m+n{\omega }_b}-\frac{cos(mt-n{\omega }_b\tilde{t})}{m-n{\omega }_b}\},\hfill \end{array}$$(C.11)

$$\begin{array}{cc}\hfill H_Y^a(t)& =\sqrt{2}{\mathrm{\Phi }}_{0}asinmt+4\sqrt{2}{\mathrm{\Phi }}_{0}a\sum _{n=1}^{+\infty }{J}_n(ne)sinmtcos(n{\omega }_b\tilde{t})\hfill \\ \hfill & -3\sqrt{2}{\mathrm{\Phi }}_{0}a{\omega }_b\sum _{n=1}^{+\infty }n{J}_n(ne)\{\frac{sin(mt+n{\omega }_b\tilde{t})}{m+n{\omega }_b}-\frac{sin(mt-n{\omega }_b\tilde{t})}{m-n{\omega }_b}\},\hfill \end{array}$$(C.12)

$$\begin{array}{cc}\hfill H_X^e(t)& ={\mathrm{\Phi }}_0\frac{1-e^2}{\sqrt{2}e}\sum _{n=1}^{+\infty }{J}_n(ne)\{\frac{n{\omega }_b-2m}{n{\omega }_b+m}cos(mt+n{\omega }_b\tilde{t})+\frac{2m+n{\omega }_b}{n{\omega }_b-m}cos(mt-n{\omega }_b\tilde{t})\},\hfill \end{array}$$(C.13)

$$\begin{array}{cc}\hfill H_Y^e(t)& =-{\mathrm{\Phi }}_0\frac{1-e^2}{\sqrt{2}e}\sum _{n=1}^{+\infty }{J}_n(ne)\{\frac{n{\omega }_b-2m}{n{\omega }_b+m}sin(mt+n{\omega }_b\tilde{t})+\frac{2m+n{\omega }_b}{n{\omega }_b-m}sin(mt-n{\omega }_b\tilde{t})\},\hfill \end{array}$$(C.14)

$$\begin{array}{cc}\hfill H_X^{\omega }(t)& ={\mathrm{\Phi }}_0\frac{\sqrt{1-e^2}}{\sqrt{2}e}\sum _{n=1}^{+\infty }{J}_n^{\prime }(ne)\{\frac{n{\omega }_b-2m}{n{\omega }_b+m}sin(mt+n{\omega }_b\tilde{t})-\frac{2m+n{\omega }_b}{n{\omega }_b-m}sin(mt-n{\omega }_b\tilde{t})\},\hfill \end{array}$$(C.15)

$$\begin{array}{cc}\hfill H_Y^{\omega }(t)& ={\mathrm{\Phi }}_0\frac{\sqrt{1-e^2}}{\sqrt{2}e}\sum _{n=1}^{+\infty }{J}_n^{\prime }(ne)\{\frac{n{\omega }_b-2m}{n{\omega }_b+m}cos(mt+n{\omega }_b\tilde{t})-\frac{2m+n{\omega }_b}{n{\omega }_b-m}cos(mt-n{\omega }_b\tilde{t})\};\hfill \end{array}$$(C.16)

for a, e and ω, the ℓ^S(t) contributions vanish while ℓ^Θ′(t) involves a linearly growing term

$$\begin{array}{c}\hfill L_X^{{\mathrm{\Theta }}^{\prime }}(t)=\frac{3{\omega }_b}{2a}{H}_X^a(T_0)\tilde{t},L_Y^{{\mathrm{\Theta }}^{\prime }}(t)=\frac{3{\omega }_b}{2a}{H}_Y^a(T_0)\tilde{t}.\end{array}$$(C.17)

As already mentioned, in the two-step approach, in the definition of h^S(t) a constant and a term that is linearly growing (in time) are irrelevant since they can be absorbed into a redefinition of the nuisance parameters. The reason we keep the constant terms in h^S(t) (given by H_X^S(T_ref) and H_X^S(T_ref)) is because some variations present resonances and the constant parts are needed to obtain finite results at the resonance points

C.2. Perturbed ELL1 model

The model for the variations of the orbital parameters for near-circular orbits including the ULDM contribution is given by:

$$\begin{array}{cc}\hfill \delta {\mathrm{\Psi }}^{\prime }(t)& =\delta {\mathrm{\Psi }}_{\mathrm{asc}}^{\prime }-\frac{3{\omega }_b}{2}(\frac{\delta a_{\mathrm{asc}}}{a}+2\alpha {\mathrm{\Phi }}_0\varrho cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}}))\overline{t}-\frac{3{\omega }_b}{4}\frac{{\dot{a}}_{\mathrm{asc}}}{a}{\overline{t}}^2\hfill \\ \hfill & +\frac{11{\omega }_b}{2m}\alpha {\mathrm{\Phi }}_0\varrho [sin(mt+\mathrm{{\rm Y}})-sin(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})],\hfill \end{array}$$(C.18)

$$\begin{array}{cc}\hfill \delta x(t)& =\delta x_{\mathrm{asc}}+{\dot{x}}_{\mathrm{asc}}\overline{t}-2\alpha x{\mathrm{\Phi }}_0\varrho [cos(mt+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})],\hfill \end{array}$$(C.19)

$$\begin{array}{cc}\hfill \delta \kappa (t)& =\delta {\kappa }_{\mathrm{asc}}+{\dot{\kappa }}_{\mathrm{asc}}\overline{t}+\alpha {\mathrm{\Phi }}_0\varrho \frac{\kappa }{{\kappa }^2+{\eta }^2}\{\frac{\kappa ({\omega }_b+2m)}{2({\omega }_b-m)}[cos(mt-{\omega }_b\overline{t}+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]\hfill \\ \hfill & +\frac{\kappa ({\omega }_b-2m)}{2({\omega }_b+m)}[cos(mt+{\omega }_b\overline{t}+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]-\frac{\eta ({\omega }_b+2m)}{2({\omega }_b-m)}[sin(mt-{\omega }_b\overline{t}+\mathrm{{\rm Y}})\hfill \\ \hfill & -sin(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]+\frac{\eta ({\omega }_b-2m)}{2({\omega }_b+m)}[sin(mt+{\omega }_b\overline{t}+\mathrm{{\rm Y}})-sin(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]\},\hfill \end{array}$$(C.20)

$$\begin{array}{cc}\hfill \delta \eta (t)& =\delta {\eta }_{\mathrm{asc}}+{\dot{\eta }}_{\mathrm{asc}}\overline{t}+\alpha {\mathrm{\Phi }}_0\varrho \frac{\eta }{{\kappa }^2+{\eta }^2}\{\frac{\kappa ({\omega }_b+2m)}{2({\omega }_b-m)}[cos(mt-{\omega }_b\overline{t}+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]\hfill \\ \hfill & +\frac{\kappa ({\omega }_b-2m)}{2({\omega }_b+m)}[cos(mt+{\omega }_b\overline{t}+\mathrm{{\rm Y}})-cos(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]-\frac{\eta ({\omega }_b+2m)}{2({\omega }_b-m)}[sin(mt-{\omega }_b\overline{t}+\mathrm{{\rm Y}})\hfill \\ \hfill & -sin(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]+\frac{\eta ({\omega }_b-2m)}{2({\omega }_b+m)}[sin(mt+{\omega }_b\overline{t}+\mathrm{{\rm Y}})-sin(m{T}_{\mathrm{asc}}+\mathrm{{\rm Y}})]\},\hfill \end{array}$$(C.21)

where $\overline{t}=t-T_{\mathrm{asc}}$, the quantities $\delta a_{\mathrm{asc}}$, $\delta x_{\mathrm{asc}}$, $\delta {\kappa }_{\mathrm{asc}}$ and $\delta {\eta }_{\mathrm{asc}}$ account for possible errors evaluated at $T_{\mathrm{asc}}$ and, as for the BT model, possible secular variations due to other effects are also included in the model. The corresponding X and Y components defined in Eq. (C.2), where T_ref = T_asc in the ELL1 model, are given by

$$\begin{array}{cc}\hfill H_X^{{\mathrm{\Psi }}^{\prime }}(t)& =\frac{11{\omega }_b}{2^{3/2}m}{\mathrm{\Phi }}_{0}sin(mt),H_Y^{{\mathrm{\Psi }}^{\prime }}(t)=\frac{11{\omega }_b}{2^{3/2}m}{\mathrm{\Phi }}_{0}cos(mt),\hfill \end{array}$$(C.22)

$$\begin{array}{cc}\hfill H_X^x(t)& =-\sqrt{2}x{\mathrm{\Phi }}_{0}cos(mt),H_Y^x(t)=\sqrt{2}x{\mathrm{\Phi }}_{0}sin(mt),\hfill \end{array}$$(C.23)

$$\begin{array}{cc}\hfill H_X^{\kappa }(t)& =\frac{{\mathrm{\Phi }}_0\kappa }{2^{3/2}({\kappa }^2+{\eta }^2)}\{\frac{({\omega }_b+2m)}{({\omega }_b-m)}[\kappa cos(mt-{\omega }_b\overline{t})-\eta sin(mt-{\omega }_b\overline{t})]\hfill \\ \hfill & +\frac{({\omega }_b-2m)}{({\omega }_b+m)}[\kappa cos(mt+{\omega }_b\overline{t})+\eta sin(mt+{\omega }_b\overline{t})]\},\hfill \end{array}$$(C.24)

$$\begin{array}{cc}\hfill H_Y^{\kappa }(t)& =\frac{{\mathrm{\Phi }}_0\kappa }{2^{3/2}({\kappa }^2+{\eta }^2)}\{\frac{-({\omega }_b+2m)}{({\omega }_b-m)}[\eta cos(mt-{\omega }_b\overline{t})+\kappa sin(mt-{\omega }_b\overline{t})]\hfill \\ \hfill & +\frac{({\omega }_b-2m)}{({\omega }_b+m)}[\eta cos(mt+{\omega }_b\overline{t})-\kappa sin(mt+{\omega }_b\overline{t})]\},\hfill \end{array}$$(C.25)

$$\begin{array}{cc}\hfill H_X^{\eta }(t)& =\frac{\eta }{\kappa }{H}_X^{\kappa }(t),H_Y^{\eta }(t)=\frac{\eta }{\kappa }{H}_Y^{\kappa }(t).\hfill \end{array}$$(C.26)

The linear terms for x,  κ and η, are zero, ℓ^S = 0, while for Ψ′ we have

$$\begin{array}{c}\hfill L_X^{{\mathrm{\Psi }}^{\prime }}(t)=-\frac{3{\omega }_b}{2^{1/2}}{\mathrm{\Phi }}_{0}cos(m{T}_{\mathrm{asc}})\overline{t},L_Y^{{\mathrm{\Psi }}^{\prime }}(t)=\frac{3{\omega }_b}{2^{1/2}}{\mathrm{\Phi }}_{0}sin(m{T}_{\mathrm{asc}})\overline{t}.\end{array}$$(C.27)

We note that the signals for the parameters η and κ exhibit a resonant term, which comes from the first resonance $m={\omega }_b$ and is the only one that survives in the e → 0 limit

C.3. Integrals for marginalisation

The integral expression for the quantities needed to evaluate the Bayes’ factor ℬ are:

$$\begin{array}{cc}\hfill u& ={\int }_0^{T_{\mathit{obs}}}\frac{{(G{h}^S)}^2(t)}{2{\sigma }_S^{2}d}\mathrm{d}t,\hfill \end{array}$$(C.28)

$$\begin{array}{cc}\hfill u_X& ={\int }_0^{T_{\mathit{obs}}}\frac{(G{A}_X^S)(t)(G{h}^S)(t)}{2{\sigma }_S^{2}d}\mathrm{d}t,\hfill \end{array}$$(C.29)

$$\begin{array}{cc}\hfill u_Y& ={\int }_0^{T_{\mathit{obs}}}\frac{(G{A}_Y^S)(t)(G{h}^S)(t)}{2{\sigma }_S^{2}d}\mathrm{d}t,\hfill \end{array}$$(C.30)

where all quantities are fixed to their fiducial values and for a given function of time f(t),

$$\begin{array}{cc}\hfill (Gf)(t)& =f(t)-\sum _{i=m_S}^2{f}_i(t)\frac{{\int }_0^{T_{\mathit{obs}}}{f}_i(\tau )f(\tau )d\tau }{{\int }_0^{T_{\mathit{obs}}}{f}_i^2(\tau )d\tau },\hfill \\ \hfill f_0(\tau )& =\frac{{\tau }^2}{T_{\mathit{obs}}}-\tau +\frac{T_{\mathit{obs}}}{6},f_1(\tau )=\tau -\frac{T_{\mathit{obs}}}{2}{f}_2(\tau )=T_{\mathrm{obs}},\hfill \end{array}$$(C.31)

where m_S = 0 if Ξ₂ ≠ 0, m_S = 1 if Ξ₂ = 0 and Ξ₁ ≠ 0, and m_S = 2 if Ξ₂ = 0 = Ξ₁ and only Ξ₀ ≠ 0

Appendix D: Resonances: going beyond P_b in the BT model

In this appendix we discuss on the importance of taking into account the ULDM effect on orbital parameters other than the orbital period P_b (which is the most studied one) for eccentric systems. We focus here on system J1903+0327, whose parameters are given in Table E.2, and compute the sensitivity curves using the effect of ULDM on the orbital parameter Θ′. The variance of δΘ′ is given in Eq. (B.6) while the model can be computed from Eq. (C.1a) where A_n and B_n are defined in Eq. (31)

Upon using Eq. (C.1b) into Eq. (C.1a) we can write δΘ′ as

$$\begin{array}{c}\hfill \delta {\mathrm{\Theta }}^{\prime }={\mathrm{\Xi }}_0+{\mathrm{\Xi }}_{1}t-2\pi \frac{{\dot{P_b}}_0}{P_b^2}{t}^2+h^{{\mathrm{\Theta }}^{\prime }}(t),\end{array}$$(D.1)

where we have used the fact that

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}=\frac{2}{3}\frac{\dot{P_b}}{P_b}+\frac{\alpha }{3}\dot{\mathrm{\Phi }},\end{array}$$(D.2)

and we have included the parameter ${\dot{P_b}}_0$ in order to describe all possible secular variations that are not caused by ULDM (e.g. the Shklovskii effect Shklovskii (1970) or the galactic acceleration Damour & Taylor (1991)). Here we collect in Ξ₀ all the constant contributions to δΘ′, regardless of whether these come from ULDM or some other effects; all the contributions growing linearly with t are included in the Ξ₁t term. The relevant ULDM signal is characterised by the functions h^Θ′(t) and h_Pb^Θ′(t) given in Appendix C as

$$\begin{array}{cc}\hfill h^{{\mathrm{\Theta }}^{\prime }}(t)& =\alpha [(H_X^{{\mathrm{\Theta }}^{\prime }}(t)-H_X^{{\mathrm{\Theta }}^{\prime }}(T_0))X+(H_Y^{{\mathrm{\Theta }}^{\prime }}(t)-H_Y^{{\mathrm{\Theta }}^{\prime }}(T_0))Y],\hfill \end{array}$$(D.3)

$$\begin{array}{cc}\hfill h_{P_b}^{{\mathrm{\Theta }}^{\prime }}(t)& =\alpha [(H_{P_b,X}^{{\mathrm{\Theta }}^{\prime }}(t)-H_{P_b,X}^{{\mathrm{\Theta }}^{\prime }}(T_0))X+(H_{P_b,Y}^{{\mathrm{\Theta }}^{\prime }}(t)-H_{P_b,Y}^{{\mathrm{\Theta }}^{\prime }}(T_0))Y],\hfill \end{array}$$(D.4)

where the functions on the right-hand-side are defined in Eq. (C.7). Here h^Θ′(t) describes the effect of ULDM on Θ′ given in Eq. (C.1a) and h_Pb^Θ′(t) takes into account only the ULDM effect on P_b obtained from the last term in Eq. (C.1a) and Eq. (D.2)

We note that the constant parameters Ξ₀ includes the error in the determination of Θ′₀, and the coefficient Ξ₁ of the linear term depends on the error in the determination of a₀. In order to obtain the sensitivity curve for α, we use the variance of δΘ′ given in Eq. (B.6) and marginalise over the two nuisance parameters Ξ₀ and Ξ₁ in the model of Θ′ using uniform priors for both of them. We assume here that any possible secular effect on P_b characterised by $\dot{P_{b0}}$ in Eq. (D.1) can be independently determined with significant precision so that we can subtract them and then set $\dot{P_{b0}}=0$ (more in Appendix E where we show that this simplification is only relevant for masses that are very close to the resonance; this is the case because ULDM produces a secular effect at m = nω_b, but not for other values of the mass)

In order to illustrate the difference in sensitivity that could be obtained by performing an analysis with h^Θ′(t) in comparison to h_Pb^Θ′(t), we restrict ourselves to the study of the signal at low masses, where the contribution of the first harmonic with n = 1 in Eq. (C.1a) dominates the signal and accounts for a resonant effect for m = ω_b. Within the scope of this work it is enough to focus on the contribution of the first n = 1 harmonic, so we neglect higher ones. The result is shown in Fig. (D.1), where we plotted the sensitivity curves for five different random realisations of the values of 𝜚 and Υ, for which we used different colours. The dashed lines were obtained using h^Θ′(t) while the solid lines correspond to h_Pb^Θ′(t). We see the effect on P_b dominates, so the estimated sensitivity obtained from both h_Pb^Θ′(t) and h^Θ′(t) are about the same. The use of h_Pb^Θ′(t) instead of h^Θ′(t) underestimate (overestimate) the sensitivity at the lower (higher) masses away from the resonant band but only by an order one factor

Appendix E: Marginalisation and the resonant peaks in the two-step approach

In this appendix we analyse the sensitivity curves near the resonances. To this end, we consider the effect on δΘ′ at low masses around the first resonant peak m = ω_b and neglect higher harmonics. For the sake of generality, we model the signal for the second step as

$$\begin{array}{c}\hfill \delta {\mathrm{\Theta }}^{\prime }(t)={\mathrm{\Xi }}_0+{\mathrm{\Xi }}_{1}t+{\mathrm{\Xi }}_2{t}^2+h^{{\mathrm{\Theta }}^{\prime }}(t)+{\ell }^{{\mathrm{\Theta }}^{\prime }}(t),\end{array}$$(E.1)

where h^Θ′(t) and ℓ^Θ′(t) are defined as in Eq. (C.2) and are given in Eqs. (D.3) and (C.17), respectively. Here the parameters Ξ_i (i = 0, 1, 2) do not include the linear function in time coming from the ULDM contribution, which is given by ℓ^Θ′(t), nor the secular quadratic term that shows up at the resonant mass

Fig. E.1 shows the sensitivity curves as a function of the mass for the system J1903+0327 around the resonant mass m = ω_b for five realisations of the parameters 𝜚 and Υ. Firstly, we see that the shape of the curve around the resonant mass depends strongly on the values of the parameters 𝜚 and Υ. Secondly, since we do not know the nuisance parameters Ξ_i (i = 0, 1, 2), an appropriate procedure to obtain a conservative sensitivity curve for α is to marginalise over such irrelevant parameters. However, in Fig. E.1 it can be seen that the differences between marginalising over the nuisance parameters, or neglecting them, are significant only in a tiny mass range near resonance.¹⁰ This is as expected, because only close to the resonance the ULDM produces a secular drift of the orbital period which can be swept in the marginalisation of the quadratic term Ξ₂

Appendix F: A one-step approach to estimate the sensitivity (cross-verification)

In this section we fit time residuals directly to the timing model (one-step approach) in order to construct the sensitivity curves. This approach is based on less robust assumptions, in particular we assume here that time residuals are due solely to the effect of ULDM on the binary’s orbital parameters and are subject to a Gaussian (white) noise distribution. We focus our analysis on the 22 ELL1(H) pulsars included in the NANOGrav 2023 dataset (we refer to Table E.1) and then extend this analysis to the three binary pulsars with higher eccentricity of Table E.2

Disadvantages. In the two-step approach we have introduced nuisance parameters to model possible errors or constant shifts in the determination of the orbital parameters and secular drifts from other sources, over which we have marginalised. The one-step method we describe here does not take into account any effect other than ULDM and noise, i.e. no such errors and secular drifts are included. Consequently, the sensitivity curves produced by this approach tend to be more optimistic (or less conservative) compared to those derived from the two-step method

Advantages. This method allows us to jointly consider variations of all orbital parameters, enhancing sensitivity—a capability lacking in the previous method. The marginalisation over 𝜚 and Υ using a Gaussian prior can be done fully analytically, giving insight into how the sensitivity depends on each variable such as observation time, noise magnitude, or cadence. Importantly, plotting the sensitivity curves is not time consuming, compared to the two-step method. Lastly, since not all binary pulsars are of the same quality in terms of their ability to detect ULDM, this method can be used to quickly determine which pulsars are the most suitable to study ULDM, allowing us to restrict ourselves to a selected subset of pulsars. This knowledge can be further exploited when using the two-step method

Timing model

In line with our assumptions, we decompose time residual R into two components

$$\begin{array}{c}\hfill R=R^{\mathrm{DM}}+R^{\mathrm{PN}}.\end{array}$$(F.1)

Here, R^DM represents the contribution of ULDM in the ELL1 model, while R^PN accounts for pulsar noise. Within the ELL1 model given in Eq. (39), at leading order in the small quantities η, κ, xω_b and $t\dot{\nu }/\nu$, the variation of δx, δΨ′,δη, δκ leads to the following expression for the time residual R^ELL1

$$\begin{array}{c}\hfill R^{\mathrm{ELL}1}=sin{\mathrm{\Psi }}^{\prime }\delta x+xcos{\mathrm{\Psi }}^{\prime }\delta {\mathrm{\Psi }}^{\prime }-x\frac{cos(2{\mathrm{\Psi }}^{\prime })+3}{2}\delta \eta +x\frac{sin(2{\mathrm{\Psi }}^{\prime })}{2}\delta \kappa ,\end{array}$$(F.2)

while R^ELL1 transforms into R^DM when δx, δΨ′,δη, δκ arise from ULDM, given by Eq. (C.18). Since we are assuming here that the variations only come from the ULDM, in such perturbation we set to zero both the error and the derivative of the orbital parameters at T_asc, i.e. δx = h^x etc

We note that the function δΨ′(t) given in Eq. (C.18) includes a constant term and a linear term (in time) that depend on the ULDM field, while the ULDM signal considered in the two-step method, namely h^Ψ′(t) defined in Eq. (C.2), does not incorporate those contributions, since they have been absorbed into the nuisance parameters. The presence of the linear term significantly enhances the sensitivity, but the result should be taken with a grain of salt because it assumes that other linear contributions are either subdominant or can be independently measured and subtracted away, whereas we expect that it is indistinguishable from an error in the semi-major axis $a_1$. Nonetheless, in order to understand the physical effect of the linear term, we present the sensitivity plots using the one-step method both with the presence of the linear term and without it. The sensitivity curves obtained in this section (without the linear term) become very close to the ones obtained from the two-step method

Noise model

Since we have already described dark matter component (we refer, for instance, to Appendix C), let us turn our attention to the noise model. We consider uncorrelated white Gaussian noise with correlation matrix

$$\begin{array}{c}\hfill C_{ai,bj}={ϵ}_a^2{\delta }_{\mathit{ab}}{\delta }_{\mathit{ij}},\end{array}$$(F.3)

where $a,b$ are pulsar indices and $i,j$ are TOAs indices, ϵ_a is the noise amplitude of pulsar a and the δs are Kronecker deltas

Concentrating on a single pulsar a, the likelihood function P_a of time residuals, denoted as {R_i}_a, where curly brackets denote a set in which i = 1, …, n_a, with parameters p_a ≡ (α,𝜚_a,Υ_a) is

$$\begin{array}{c}\hfill P_a\left({\{R_i\}}_a|p_a\right)=\frac{1}{{(2\pi )}^{\frac{n_a}{2}}{ϵ}_a^{n_a}}exp[-\frac{1}{2}\sum _{i=1}^{n_a}\frac{{(R_{\mathit{ai}}-R_{\mathit{ai}}^{\mathrm{DM}})}^2}{{ϵ}_a^2}].\end{array}$$(F.4)

The likelihood function combined for all pulsars a = 1, …, N_p is as follows

$$\begin{array}{c}\hfill P(\{R_{\mathit{ai}}\}|p)\equiv {\mathrm{\Pi }}_{a=1}^{N_p}{P}_a\left({\{R_i\}}_a|p_a\right).\end{array}$$(F.5)

Sensitivity limit (δ-function priors)

We apply Bayes’ theorem to Eq. (F.5)

$$\begin{array}{c}\hfill P(p|\{R_{\mathit{ai}}\})=P(\{R_{\mathit{ai}}\}|p)\frac{P_0\left(p\right)}{P_0(\{R_{\mathit{ai}}\})},\end{array}$$(F.6)

where P₀(p) denotes the prior distribution for the unknown parameters p ≡ (α, {𝜚_a},{Υ_a}), and P₀({R_ai}) is the evidence. In what follows, we consider δ-function priors

$$\begin{array}{c}\hfill P_0\left(p\right)={\mathrm{\Pi }}_{a=1}^{N_p}{P}_0\left(p_a\right)={\mathrm{\Pi }}_{a=1}^{N_p}\delta ({\varrho }_a-{\varrho }_a^{\prime })\delta ({\mathrm{{\rm Y}}}_a-{\mathrm{{\rm Y}}}_a^{\prime }),\end{array}$$(F.7)

which leads to the probability distribution function of α

$$\begin{array}{c}\hfill P(\alpha |\{R_{\mathit{ai}}\})=\frac{1}{\sqrt{2\pi }{\sigma }_{\alpha }}exp\left[-\frac{{(\alpha -{\alpha }_{\mathrm{peak}})}^2}{2{\sigma }_{\alpha }^2}\right],\end{array}$$(F.8)

where the peak and standard deviation are

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{peak}}\equiv {\sigma }_{\alpha }^2\sum _{a=1}^{N_p}\frac{{\alpha }_{\mathrm{peak},a}}{{\sigma }_{\alpha ,a}^2},\frac{1}{{\sigma }_{\alpha }^2}\equiv \sum _{a=1}^{N_p}\frac{1}{{\sigma }_{\alpha ,a}^2},\end{array}$$(F.9)

where

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{peak},a}\equiv \frac{{\mathit{R}}_a.{\mathit{Q}}_a}{{\mathit{Q}}_a^2},{\sigma }_{\alpha ,a}\equiv \frac{{ϵ}_a}{\sqrt{{\mathit{Q}}_a^2}},\end{array}$$(F.10)

where Q_a ≡ R_a^DM/α is a function describing the ULDM signal, but its magnitude is independent on the coupling constant (but still depends on 𝜚_a and Υ_a). With this in hand, we are in a position to establish a simple sensitivity limit for the detection. Following the well-known 3σ-rule for the Gaussian distribution, to ensure the exclusion of the null case the requirement

$$\begin{array}{c}\hfill |{\alpha }_f|>3{\sigma }_{\alpha },\end{array}$$(F.11)

where α_f represents the genuine value in nature, should be met. Indeed, if we use the Bayes’ factor to determine the detection limit, i.e

$$\begin{array}{c}\hfill \mathcal{B}=\frac{\int P(\{R_{\mathit{ai}}\}|{\alpha }_f,\{{\varrho }_a\},\{{\mathrm{{\rm Y}}}_a\})P_0(\{{\varrho }_a\},\{{\mathrm{{\rm Y}}}_a\}){\mathrm{\Pi }}_{a=1}^{N_p}d{\varrho }_{a}d{\mathrm{{\rm Y}}}_a}{P(\{R_{\mathit{ai}}\}|{\alpha }_f=0)},\end{array}$$(F.12)

where P({R_ai}|α_f = 0) is the likelihood function of pure white Gaussian noise (that is, the hypothesis that R_a contains α_f = 0), and we apply R_a = R_a^PN + α_fQ_a , we find

$$\begin{array}{c}\hfill \mathcal{B}\equiv {\mathrm{\Pi }}_{a=1}^{N_p}{\mathcal{B}}_a\equiv {\mathrm{\Pi }}_{a=1}^{N_p}exp[{\alpha }_f^2\frac{{\mathit{Q}}_a^2}{2{ϵ}_a^2}+\frac{{\mathit{R}}_a^{\mathrm{PN}}·{\mathit{Q}}_a}{{ϵ}_a^2}{\alpha }_f].\end{array}$$(F.13)

We note that this is a noise-dependent quantity. To get rid of this dependency, we average over multiple noise realisations. This is not an unambiguous procedure; one way to average is as follows

$$\begin{array}{c}\hfill \overline{\mathcal{B}}\equiv \int \mathcal{B}P(\{R_{\mathit{ai}}^{\mathrm{PN}}\}){\mathrm{\Pi }}_{a=1}^{N_p}d{R}_{a1}^{\mathrm{PN}}\cdots d{R}_{a{n}_a}^{\mathrm{PN}}=exp\left[\sum _{a=1}^{N_p}\frac{{\mathit{Q}}_a^2}{{ϵ}_a^2}{\alpha }_f^2\right],\end{array}$$(F.14)

or we can average e.g. logarithms as

$$\begin{array}{c}\hfill ln\mathbb{B}\equiv \int ln\mathcal{B}P(\{R_{\mathit{ai}}^{\mathrm{PN}}\}){\mathrm{\Pi }}_{a=1}^{N_p}d{R}_{a1}^{\mathrm{PN}}\cdots d{R}_{a{n}_a}^{\mathrm{PN}}=\sum _{a=1}^{N_p}\frac{{\mathit{Q}}_a^2}{2{ϵ}_a^2}{\alpha }_f^2,\end{array}$$(F.15)

yielding, in that order

$$\begin{array}{cc}\hfill |{\alpha }_f|& =\sqrt{ln\overline{\mathcal{B}}}{\sigma }_{\alpha }\doteq 2.63{\sigma }_{\alpha },\hfill \end{array}$$(F.16)

$$\begin{array}{cc}\hfill |{\alpha }_f|& =\sqrt{2ln\mathbb{B}}{\sigma }_{\alpha }\doteq 3.72{\sigma }_{\alpha }.\hfill \end{array}$$(F.17)

with $\overline{\mathcal{B}}=1000$ and 𝔹 = 1000. The take-home message is that, when we consider the δ-function priors, the sensitivity limit is given by σ_α, with a numerical factor close to 3. Throughout the following discussion, we adhere to the factor of 3

We show a few sensitivity lines for two selected pulsars in Figs. F.1 and F.2, and the combined result of all 22 ELL1(H) pulsars in Fig. F.3. Indeed, we can observe that the so-called linear term has a significant effect on the estimated sensitivity derived from δΨ′, both qualitatively (the shape of the curve) and quantitatively (the curve’s normalisation). Fitting only the periodic component gives a less optimistic result and more similar to the results obtained by the two-step approach

The parameters δη and δκ are of particular interest due to the resonance occurring when m = ω_b. The depth of the signal relative to the line is approximately determined by the ω_bT factor, consistently around three to four orders of magnitude across all studied pulsars. This also suggests that the depth scales with the length of the observation

Large eccentricities

For binary pulsars described by the BT model, the dominant part of the time residual (ignoring the contribution of δK) comes from Eq. (B.1), whereas the variations due to ULDM can be obtained from (C.1) using the definition of α_b and η_b in (3) and recalling that η_b ≈ β_b. Explicitly,

$$\begin{array}{cc}\hfill R^{\mathrm{BT}}& =\frac{\delta x}{x}[{\eta }_{b}sin{E}^{\prime }+{\alpha }_b(cos{E}^{\prime }-e)]-\delta e[{\alpha }_b+\frac{e{\eta }_{b}sin{E}^{\prime }}{1-e^2}+\frac{sin{E}^{\prime }({\alpha }_{b}sin{E}^{\prime }-{\eta }_{b}cos{E}^{\prime })}{1-ecos{E}^{\prime }}]\hfill \\ \hfill & +\delta \omega x[-sin{E}^{\prime }sin\omega {(1-e^2)}^{1/2}+cos\omega (cos{E}^{\prime }-e)]-\frac{\delta {\mathrm{\Theta }}^{\prime }({\alpha }_{b}sin{E}^{\prime }-{\eta }_{b}cos{E}^{\prime })}{1-ecos{E}^{\prime }},\hfill \end{array}$$(F.18)

where the variations of δx/x = δa/a, δω, δe, and δΘ′ are given by Eq. (C.1) and, as in the ELL1 case, in this appendix we ignore possible errors in the orbital parameters and secular drifts from other sources, and include the linearly growing term in δΘ′ that depends on the ULDM field

Five sensitivity curves for δx, δω, δe, δΘ′ are displayed in Fig. F.4. The shape of the sensitivity curves derived from δΘ′ is non-trivial and can vary depending on the value of Υ. The sensitivity obtained from the combination of all parameters, primarily influenced by δΘ′, for all three high-eccentricity pulsars is depicted in Fig. F.5, focusing solely on the first three harmonics. In such figure we also show the sensitivity curves for δΘ′_Pb defined in Eq. (65). In agreement to what was shown in Appendix D, the curves for δΘ′_Pb only differ from the ones for δΘ′ by order-one factors. The peculiarities related to the linear terms that we mentioned above for δΨ′ are also present for δΘ′ in the BT model and only the ones with no linear terms resemble those found in the two-step method. The combination of all parameters, primarily influenced by δΘ′, for all three high-eccentricity pulsars is depicted in Fig. F.6, focussing solely on the first three harmonics

Sensitivity limit (Gaussian priors)

We aim to compute $\overline{\mathcal{B}}$ in the case of Rayleigh (for 𝜚_a) and uniform (for Υ_a) priors:

$$\begin{array}{c}\hfill P_0\left(p\right)\equiv {\mathrm{\Pi }}_{a=1}^{N_p}{P}_0\left(p_a\right)={\mathrm{\Pi }}_{a=1}^{N_p}\frac{1}{\pi }{e}^{-{\varrho }_a^2}{\varrho }_a.\end{array}$$(F.19)

In order to compute $\overline{\mathcal{B}}$, we need to compute the integrals over 𝜚_a, Υ_a and $R_{\mathit{ai}}^{\mathrm{PN}}$. Practically, it is advantageous to first perform the integration over the noise and after that over the variables $X_a\equiv \sqrt{2}{\varrho }_{a}cos{\mathrm{{\rm Y}}}_a$ and $Y_a\equiv \sqrt{2}{\varrho }_{a}sin{\mathrm{{\rm Y}}}_a$, leading to two separate Gaussian integrals. For this we decompose R^DM into its X and Y components as

$$\begin{array}{c}\hfill R^{\mathrm{DM}}=\alpha (Q_{X}X+Q_{Y}Y),\end{array}$$(F.20)

where for the ELL1 case we have

$$\begin{array}{c}\hfill Q_{X(Y)}^{\mathrm{ELL}1}=sin{\mathrm{\Psi }}^{\prime }{A}_{X(Y)}^x+xcos{\mathrm{\Psi }}^{\prime }[A_{X(Y)}^{{\mathrm{\Psi }}^{\prime }}+L_{X(Y)}^{{\mathrm{\Psi }}^{\prime }}]-x\frac{cos(2{\mathrm{\Psi }}^{\prime })+3}{2}{A}_{X(Y)}^{\eta }+x\frac{sin(2{\mathrm{\Psi }}^{\prime })}{2}{A}_{X(Y)}^{\kappa },\end{array}$$(F.21)

and for the BT model

$$\begin{array}{cc}\hfill Q_{X(Y)}^{\mathrm{BT}}& =\frac{A_{X(Y)}^a}{a}[{\eta }_{b}sin{E}^{\prime }+{\alpha }_b(cos{E}^{\prime }-e)]\hfill \\ \hfill & -A_{X(Y)}^e[{\alpha }_b+\frac{e{\eta }_{b}sin{E}^{\prime }}{1-e^2}+\frac{sin{E}^{\prime }({\alpha }_{b}sin{E}^{\prime }-{\eta }_{b}cos{E}^{\prime })}{1-ecos{E}^{\prime }}]\hfill \\ \hfill & +A_{X(Y)}^{\omega }x[-sin{E}^{\prime }sin\omega {(1-e^2)}^{1/2}+cos\omega (cos{E}^{\prime }-e)]\hfill \\ \hfill & -\frac{[A_{X(Y)}^{{\mathrm{\Theta }}^{\prime }}+L_{X(Y)}^{{\mathrm{\Theta }}^{\prime }}]({\alpha }_{b}sin{E}^{\prime }-{\eta }_{b}cos{E}^{\prime })}{1-ecos{E}^{\prime }}.\hfill \end{array}$$(F.22)

After integration over noise and X and Y we obtain

$$\begin{array}{c}\hfill \overline{\mathcal{B}}=exp\left(\frac{1}{4}\frac{{\alpha }^4}{{\mathrm{\Sigma }}^4}\right),\end{array}$$(F.23)

while Σ is defined as

$$\begin{array}{c}\hfill \frac{1}{{\mathrm{\Sigma }}^4}\equiv \sum _{a=1}^{N_p}\frac{1}{{\mathrm{\Sigma }}_a^4},\end{array}$$(F.24)

with

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_a\equiv \frac{1}{\sqrt[4]{2}}\frac{1}{\sqrt{{\varrho }_{f,a}}}\frac{{ϵ}_a}{\sqrt[4]{{({\mathit{Q}}_{X,a}·{\mathit{S}}_a)}^2+{({\mathit{Q}}_{Y,a}·{\mathit{S}}_a)}^2}},\end{array}$$(F.25)

where the vector S is defined as

$$\begin{array}{c}\hfill \mathit{S}=\sqrt{2}({\mathit{Q}}_{X}cos{\mathrm{{\rm Y}}}_f+{\mathit{Q}}_{Y}sin{\mathrm{{\rm Y}}}_f),\end{array}$$(F.26)

while in the last equation we have ignored the index a, which indicates a particular pulsar. The vectors are indicated by an arrow above the symbols and a dot (⋅) stands for the scalar product

Let us compare the sensitivity limits obtained for both priors (δ-function and Gaussian) when considering a single pulsar

$$\begin{array}{cc}\hfill |{\alpha }_f{|}^{\delta }& =\sqrt{ln\overline{\mathcal{B}}}\frac{1}{{\varrho }_f}\frac{ϵ}{\sqrt{{\mathit{S}}^2}},\hfill \end{array}$$(F.27)

$$\begin{array}{cc}\hfill |{\alpha }_f{|}^G& =\sqrt[4]{4ln\overline{\mathcal{B}}}\frac{1}{\sqrt[4]{2}}\frac{1}{\sqrt{{\varrho }_f}}\frac{ϵ}{\sqrt[4]{{({\mathit{Q}}_X·\mathit{S})}^2+{({\mathit{Q}}_Y·\mathit{S})}^2}}.\hfill \end{array}$$(F.28)

For 𝜚_f large (small) the δ-sensitivity limit outperforms (does not outperform) the Gaussian prior. This is an intuitive result, for example when 𝜚_f is large, the signal should be strong, so the δ-prior should produce a low-lying sensitivity line, but a Gaussian prior tends to underestimate this signal by leaning more towards the average

We plot sensitivities obtained from both priors for random realisations of 𝜚_f (Rayleigh distribution) and Υ_f (flat distribution). The result is depicted on Fig. F.7. We can notice that most of the δ-sensitivities are above the sensitivities derived from the Gaussian distribution. This result comes out opposite to the one we arrived at in the main body of the paper. We explain this as follows—for a given value of Υ_f and m, we find the value of 𝜚_f^EQ so that both sensitivities come out the same

$$\begin{array}{c}\hfill {\varrho }_f^{\mathrm{EQ}}=\sqrt{\frac{ln\overline{\mathcal{B}}}{2}}F(m,\mathrm{{\rm Y}}),\end{array}$$(F.29)

with the function F(m, Υ) given as

$$\begin{array}{c}\hfill F(m,\mathrm{{\rm Y}})\equiv \frac{\sqrt{{({\mathit{Q}}_X·\mathit{S})}^2+{({\mathit{Q}}_Y·\mathit{S})}^2}}{\sqrt{{\mathit{S}}^2}}.\end{array}$$(F.30)

For the method described in the main text, the function F would be different, but the numerical factor $\sqrt{(ln\overline{\mathcal{B}})/2}$ is the same. The numerical values may therefore differ slightly, but this is sufficient to change the statistics of the sensitivity lines. For $\overline{\mathcal{B}}=1000$, Eq. (F.30) and 30 values of Υ (on the interval [0, 2π) with a constant step), we find that 𝜚_f^EQ ≥ 1.3. However, this implies P(𝜚 ≥ 1.3)≐0.18. Let us remind that the median of the distribution is $\varrho =\sqrt{ln2}\doteq 0.833$. Thus, the probability that the sensitivity curve for the δ-prior is lower than that of the Gaussian priors is less than 18%. If the factor F was just slightly smaller, it could lead to a probability larger than 50%. We expect the latter be the case for the method in the main text

All Tables

All Figures