Secular effects of ultralight dark matter on binary pulsars

Dark matter (DM) can consist of very light bosons behaving as a classical scalar field that experiences coherent oscillations. The presence of this DM field would perturb the dynamics of celestial bodies, either because the (oscillating) DM stress tensor modifies the gravitational potentials of the galaxy or if DM is directly coupled to the constituents of the body. We study secular variations of the orbital parameters of binary systems induced by such perturbations. Two classes of effects are identified. Effects of the first class appear if the frequency of DM oscillations is in resonance with the orbital motion; these exist for general DM couplings including the case of purely gravitational interaction. Effects of the second class arise if DM is coupled quadratically to the masses of the binary system members and do not require any resonant condition. The exquisite precision of binary pulsar timing can be used to constrain these effects. Current observations are not sensitive to oscillations in the galactic gravitational field, though a discovery of pulsars in regions of high DM density may improve the situation. For DM with direct coupling to ordinary matter, the current timing data are already competitive with other existing constraints in the range of DM masses ∼ 10 − 22 − 10 − 18 eV. Future observations are expected to increase the sensitivity and probe new regions of parameters. equations for quadratic coupling are obtained by replacing


I. INTRODUCTION
Dark matter (DM), a key ingredient in the success of the standard cosmological model, remains a mysterious component of our Universe. Its properties are still largely unknown, despite impressive progress in characterizing and testing the signatures of different DM candidates.
In this paper we consider the possibility that DM is composed of ultralight scalar particles. Natural candidates of this kind are axionlike and dilatonlike particles, whose theoretical and phenomenological motivation can be found in [1][2][3]. Due to huge phase space occupation numbers, DM in these models is described as a classical scalar field Φ [4][5][6].
Several observables have been identified to probe different mass ranges of ultralight dark matter (ULDM). Current constraints from observations of the cosmic microwave background and large scale structure on linear scales robustly exclude m Φ ≲ 10 −24 eV [7] (with a factor of 10 improvement in m Φ expected in the future [8]). The analysis of the Lyman-α forest [9][10][11] and galactic dynamics [12][13][14] further excludes m Φ ≲ 10 −22 eV and shows a tension between the data and predictions of ULDM models with m Φ ∼ 10 −22 −10 −21 eV. Complementary constraints in this mass range can come from pulsar-timing arrays [15,16]. Yet higher masses can be probed by the dynamics of star clusters [17] or the spectral shape of the 21-cm absorption feature [18,19] (see [20] for the discussion of the robustness of 21 cm bounds). Finally, ULDM masses up to m Φ ∼ 10 −18 eV might be explored using 21 cm intensity mapping [21]. Independently of their role as DM, the existence of ultralight bosons has potentially observable signatures due to the superradiance instability of rotating black holes [22][23][24][25][26]. These tests may probe (or rule out) the existence of scalar degrees of freedom with masses in the range m Φ ∼ 10 −21 −10 −10 eV.
The previous bounds are based only on the gravitational interaction of ULDM and hence apply to any ULDM model. But DM candidates may also be directly coupled to standard model (SM) fields through nongravitational interactions. An interesting possibility is a universal dilatonlike coupling between ULDM field Φ and ordinary matter which respects the weak equivalence principle (WEP). Still, this interaction violates the strong equivalence principle (SEP) with nonuniversal effective couplings to Φ being generated for objects with large gravitational self-energy [27]. Powerful constraints on the SEP violation and other consequences of universally coupled ultralight scalar fields (non-necessarily DM) come from the tests within the Solar System [28][29][30][31][32][33]. In addition, for the ULDM case, one can infer constraints on the direct coupling from the bounds on the stochastic gravity wave background [16,34] (see the discussion in [35]).
In our previous work [35] we showed that competitive constraints on ULDM in the mass range 10 −22 eV ≲ m Φ ≲ 10 −18 eV are obtained from the timing of millisecond binary pulsars. That work focused on ULDM universally coupled to the binary system members and studied the secular change in the orbital period induced by this coupling. The study was generalized in [36,37] to the case of ultralight vector and tensor DM. Recently Ref. [38] considered instantaneous variations of the binary orbit due to ULDM. These are, however, strongly suppressed compared to the secular effects.
In this paper we extend the analysis of [35] by including the possibility of nonuniversal couplings and studying the secular effect of ULDM on all orbital parameters appearing in the pulsar timing model [39][40][41]. It is worth stressing that the couplings of ULDM to the two components of the binary pulsar system are naturally expected to be different even if at the fundamental level ULDM couples to the SM fields universally. This is a consequence of the SEP violation and large gravitational binding energy of the pulsar constituting a few tens of percent of its mass. Irrespectively of whether the companion is a normal star, a white dwarf, or another neutron star (NS) with a different mass, its binding energy will generically be different.
We will see that nonuniversality of ULDM-pulsar coupling leads to a rich phenomenology that can be tested by current and future observations. In particular, in contrast to the universal coupling case, we find a nonzero drift of the orbital period for binaries with circular orbits. Given that the orbits of most of the binary pulsars are close to circular 1 [43,44], this greatly increases the number of systems that can be used in the analysis.
The paper is organized as follows. In Sec. II we summarize the relevant properties of the ULDM field in the Milky Way halo. The effective Lagrangian describing interaction of a binary system with ULDM is presented in Sec. III. We first deal with the case of pure gravitational interaction (Sec. III A) and then include the direct coupling (Sec. III B). Section IV contains the derivation of the key equations governing the secular evolution of the binary system. The case of only gravitational interaction is considered in Sec. IVA. We allow for direct coupling in Sec. IV B, studying first the ULDM-induced motion of the binary barycenter (Sec. IV B 1), next turning to the evolution of the orbital elements for linearly (Sec. IV B 2) and quadratically (Sec. IV B 3) coupled ULDM. In Sec. V we provide numerical estimates for the size of possible effects and compare them to the current sensitivity. We first consider in Sec. VA the effects that appear when the frequency of ULDM oscillations is in resonance with the orbital motion of the binary. We separately discuss the cases of pure gravitational interaction (Sec. VA 1), universal direct coupling (Sec. VA 2), and nonuniversal coupling (Sec. VA 3). In Sec. V B we study the timing constraints when the ULDM mass deviates from the resonant value. In Sec. V C we explore nonresonant secular effects that appear for quadratically coupled ULDM and lead to constraints in a wide range of ULDM masses. Section VI is devoted to conclusions and outlook. Two appendixes are added to make the paper self-contained. Appendix A summarizes the Keplerian dynamics and the formalism of osculating orbits. Appendix B contains the osculating orbit equations for a binary interacting with ULDM prior to extraction of secular contributions.

II. ULTRALIGHT DARK MATTER IN THE HALO
We assume that DM is described by a scalar field Φ minimally coupled to gravity. For gravity, we use the standard Einstein-Hilbert action, while for DM the action reads 2 We are interested in the mass range m Φ ∼ 10 −22 − 10 −18 eV. The occupation numbers for modes within the halo of a typical galaxy are sufficiently high that the classical field theory description is applicable. In the nonrelativistic limit, the field oscillates in time with frequency set by the mass, whereas its spatial gradients are suppressed. The rapid oscillations are conventionally factored out and the field is written as One can compute the stress energy tensor T μν of the scalar field and take the average hÁ Á Ái over the "fast" oscillations of frequency m Φ . At leading order in gradients of Φ 0 and Ξ one finds This leads to the identification of as the DM density and velocity, respectively. The spatial stress tensor reads where and dots stand for terms suppressed by gradients. We see that the ULDM field is characterized by a large oscillating pressure. This, however, vanishes upon averaging, so that at large time and length scales the field behaves as a pressureless fluid, similar to the standard cold DM. On the other hand, the gradient terms neglected in (6) survive the averaging and lead to deviations from cold DM behavior at scales comparable to the de Broglie wavelength of the field λ dB ≡ 2π=ðm Φ VÞ. The description of self-gravitating ULDM in the nonrelativistic regime which is valid at all scales is provided by the Schrödinger-Poisson system of equations for the complex amplitude Ψ and the Newtonian potential ϕ; see e.g., [6] for review. In a virialized halo one expects the DM field to consist of a superposition of waves with random phases. The amplitudes of the modes are expected to follow approximately the Maxwell distribution with the dispersion V 0 determined by the typical velocity at a given location in the halo. Such superposition produces in space an interference pattern of wave packets with characteristic size λ dB =2. Inside an individual wave packet the field performs coherent oscillations with the period t osc ¼ 2π=m Φ and the slowly changing amplitude and phase. Each coherence patch is characterized by a local DM density and velocity, which exhibit large (order-one) fluctuations between the patches. These expectations are born out by numerical simulations [45][46][47][48] that show a granular structure of DM distribution in halos characteristic of wave interference. Due to the motion of wave packets, the field at a fixed point in space loses coherence after time t coh ∼ λ dB =ð2VÞ. The local DM density averaged over timescales much bigger than t coh reproduces the smooth halo profile, whereas the local velocity averages to zero.
Let us compare the scales introduced above to the typical characteristics of binary pulsar systems. The ULDM oscillation period is which for m Φ ∼ 10 −22 −10 −18 eV is comparable to pulsar orbital periods. This allows for a resonant interaction between ULDM and a binary. The de Broglie wavelength of the ULDM field is set by the velocity dispersion V 0 of the halo which can be estimated from the virial velocity.
For the Milky Way we assume V 0 ∼ 10 −3 [49], so that This is typically much larger than the size of the orbit of the binary systems we consider, with semimajor axes always satisfying a ≲ 10 8 km. Also the coherence time, is larger than or comparable to the observation time.
Therefore we can assume that the binary system is located entirely inside a single ULDM coherence patch and the field oscillations remain coherent on the relevant timescales. It will be convenient to work in the binary rest frame. We will need the field together with its spatial gradients in the vicinity of the binary and will neglect higher order spatial derivatives. Thus we write where ⃗ V now denotes the relative velocity of DM with respect to the binary barycenter (BB) and we have introduced the vector to characterize the spatial variation of the field amplitude. The granular structure of the ULDM halo described above implies that the absolute value of ⃗ S is of order 2λ −1 dB . We will consider ρ DM , ⃗ S, ⃗ V, and ϒ in (11) as constant over the binary system size.
Most of the millisecond binary pulsars observed to date are located in the ∼2 kpc neighborhood of the Solar System. The Milky Way halo models give DM density in this regionρ DM ∼ 0.3-0.6 GeV=cm 3 [50][51][52]. In our numerical estimates we will use a conservative proxȳ ρ DM ¼ 0.3 GeV=cm 3 . It is worth mentioning, however, that the actual ULDM density in the vicinity of the binary may differ by a factor of a few fromρ DM due to large fluctuations between the coherence patches.
Throughout our analysis we will neglect any effect that the binary system may produce on the ULDM distribution (that is, the backreaction effect). To estimate when this approximation is valid, we follow Ref. [6] and treat the ULDM coherence patches as quasiparticles of mass The interaction between the binary and ULDM may be viewed as the scattering of quasiparticles that lasts for the coherence time t coh . It leads to the change of quasiparticle momentum and energy, are the time derivatives of the binary momentum and energy. The backreaction is negligible as long as Δp and ΔE are small compared to the typical quasiparticle momentum p 0 ∼ M qp V 0 and its total energy (including the rest mass) E 0 ∼ M qp . In this way we arrive at the conditions These will be verified a posteriori upon estimating _ p b and _ E b in Sec. V. Note that they are easier to satisfy for lighter ULDM, which is one of the reasons to restrict our study to the range m Φ ≲ 10 −18 eV.

III. INTERACTION OF DARK MATTER WITH COMPACT BODIES
To model the binary system, we treat the pulsar and its companion as point particles. These will always interact with DM gravitationally. In addition, there may be an interaction between the bodies and Φ from a direct coupling. In all cases, in the limit of point particles the effective action for the bodies in the system will be where τ A is the proper time and M A ðΦÞ is a function of the DM field [27]. In order to characterize the main potentially observable effects that the DM field produces on a binary system, it will be enough to work at leading (Newtonian) order in the velocities of the members of the system, and at first order in the corrections due to the presence of Φ (in an analogous way as customarily done for gravitational waves). In the cases where the description of the orbital motion requires the inclusion of relativistic corrections, our results can be generalized in a systematic way by means of the post-Newtonian formalism [53].
A. Pure gravitational interaction: Oscillations in the spacetime curvature Let us start with the case when there is no direct coupling between DM and the bodies, so that their masses are independent of Φ, M A ðΦÞ ≡ M A . Following [15,35] we begin by computing the perturbation of the metric produced by the oscillating wave (11) in the Newton gauge, We are interested in resonant effects due to the oscillating pressure (7). Thus we can neglect the gradient contributions into the DM energy-momentum tensor which are suppressed by the DM velocity. The linearized Einstein equations read where G is Newton's gravitational constant. Their solution splits into the time-independent part and an oscillatory homogeneous contribution, 3 ψ ¼ψð⃗ xÞ þψðtÞ; ð17Þ and similarly for ϕ. The time-independent part locally has the form where the constant vector A i and tensor B ij correspond to the acceleration and tidal forces produced by the galactic gravitational field. These are standard effects in general relativity (GR) [54] and are taken into account in the analysis of pulsar timing. 4 In what follows we focus on the effects due to the oscillating nature of the ULDM field, and thus omit the time-independent part of the gravitational potentials. From (16b) we find that the oscillating part of the potential ψ satisfiesψ As will be clear in a moment, this expression is sufficient to study the dynamics of the binary system within our 3 More precisely, the oscillatory contribution has a weak dependence on the spatial coordinates suppressed by the DM density gradients, which we neglect in our derivation. 4 We note that in the case of ULDM, large density fluctuations between coherence patches can lead to an additional acceleration of the binary as a whole, on top of the acceleration in the largescale galactic field. This effect, however, appears degenerate with the gravitational pull of other celestial bodies.
approximations. In particular, we will not need the Newtonian potentialφ.
The equations of motion for the binary are conveniently formulated in the Fermi normal coordinates ft F ; ⃗ ξg associated with the BB rest frame. This corresponds to the geodesic motion of a free particle in the metric (15) with the initial conditions ⃗ Including the Newtonian potential ϕ N generated by the orbiting bodies, the metric in the new coordinates reads which is independent ofφðtÞ. The dynamics of the binary system is obtained from the geodesic deviation equations, where one adds the modifications to the curvature δR μ νρσ produced by the DM field [55]. Defining ⃗ r ≡ ⃗ ξ 1 − ⃗ ξ 2 , one gets an extra acceleration δr i ¼ −δR i 0j0 r j . Since the perturbation of the Riemann tensor is gauge independent at the linearized level, we can evaluate it in the Newton gauge (15). In this way we arrive aẗ where the dot stands for the derivative with respect to t F and M T ≡ M 1 þ M 2 . Alternatively, from the nonrelativistic limit of Eq. (14) one can derive the Lagrangian whose variation reproduces Eq. (22). Notice that the orbital dynamics decouples from the center-of-mass motion. The BB coordinate obeys the equation̈⃗R ¼ −̈ψ ⃗ R, which admits the trivial solution ⃗ R ¼ 0. In other words, BB indeed follows the geodesic which we have chosen as the origin of our Fermi frame. We are going to see that the situation is different if ULDM couples nonuniversally to the bodies in the binary.
B. Direct coupling of ULDM to matter: Violation of the equivalence principle We now include the possibility of direct coupling encapsulated by the dependence of the masses of the bodies on the ULDM field. We assume that the change of the masses due to this coupling is small, so that it can be Taylor expanded in the amplitude of the field. We analyze the cases of linear and quadratic couplings, whereM A is the mass of the body in the absence of the field. Part of the sensitivities α A and β A arises from the coupling of the field Φ with the SM fundamental fields. This coupling may be universal (that is, given by the same interaction with the same coupling constant for all fields) or not. For instance, any direct coupling of Φ with the Higgs field will produce a nonuniversal coupling to the masses of electrons, protons, and neutrons (the main sources of mass in the stars); see e.g., [56].
The microscopic couplings are typically inversely proportional to some high-energy scale Λ, which sets the cutoff of the effective low-energy theory. Generically, one expects on dimensional grounds α A ∼ Λ −1 , β A ∼ Λ −2 . In this case, the linear coupling will dominate the low-energy phenomenology and the quadratic interaction can be neglected. However, the linear coupling may be forbidden if the underlying theory possesses an appropriate symmetry, for example, is invariant under Φ ↦ −Φ. Then the dominant term in the interaction between Φ and SM will be quadratic. We study the two types of couplings in (25) separately, as representatives of these two classes of theories.
For compact objects, such as NSs, the mass also has a contribution coming from the self-gravity of the body. This contribution couples to Φ differently from matter and thus generally violates the SEP even if the WEP is satisfied. To be more specific, consider the case when Φ couples to the SM fields universally through a Jordan-Fierz metric g μν ¼ A 2 ðΦÞg μν and define If the self-gravity of the body is sufficiently weak, the linear sensitivity can be expanded 5 as [27] 5 The calculation in [27] assumes a constant background field far from the body, rather than an oscillating wave. Since in our case the timescale of the oscillations is large compared to the scales relevant for the inner structure of the body and the corresponding near field physics, the evolution of the field can be considered as adiabatic.
where c 1 < 0, c 2 , etc., are order-one coefficients and s A ≃ GM A =r A is the ratio of the gravitational binding energy to the rest mass of the body (here r A is the size of the body). In case α ¼ 0 a similar expression holds for the quadratic sensitivity β A upon replacing α by β. The binding energy is negligible for ordinary stars and white dwarfs, whereas for neutron stars the typical values are s A ≃ 0.1-0.2 depending on the NS mass [27]. This implies that the two bodies comprising a binary pulsar system will typically have a 10% difference in their sensitivities, and similarly for the case of quadratic coupling. 6 The preceding logic does not take into account nonperturbative effects that can further boost the difference in the sensitivities of compact bodies. A well-known example of such an effect is spontaneous scalarization [57,58] which typically occurs for sufficiently large negative β. In this case a NS develops a nonzero scalar charge α A ≠ 0 even if the fundamental linear coupling of the SM fields to Φ vanishes, α ¼ 0. The physical origin of this effect is easy to understand. The coupling of the scalar field to nuclear matter inside the star gives a contribution to the effective scalar mass, where T m μ μ ¼ −ρ m þ 3p m is the trace of the matter energymomentum tensor and ρ m and p m are the density and pressure in the NS interior. For ρ m > 3p m and large negative β this mass becomes tachyonic, implying an instability of the trivial field configuration Φ ¼ 0 inside the NS. A necessary and sufficient condition for the development of the instability is that the corresponding "Compton wavelength" should be smaller than the size of the body, 1= ffiffiffiffiffiffiffiffiffiffi ffi jm 2 eff j p ≲ r A . In this case the NS generates a nontrivial scalar profile which corresponds to an effective scalar charge. 7 As it is clear from the above reasoning, spontaneous scalarization can happen also for β > 0 if the nuclear equation of state is such that the pressure in the NS interior is large, p m > ρ m =3 [61][62][63][64]. If instead p m is less than ρ m =3, positive values of β will increase the mass of the scalar field inside the NS. This expels the scalar field from the star interior and reduces the sensitivities [57]. If m 2 eff ≫ 1=r 2 A , only the outer layer of the NS of thickness ∼m −1 eff will interact with the scalar field. The upshot of this discussion is that nonperturbative phenomena can dramatically modify the interaction of the scalar field Φ with NSs, compared to its bare coupling to the fundamental SM fields. This strengthens the role of pulsar timing as a probe of scalar-tensor theories complementary to the laboratory tests or tests performed in a weak gravitational field [58,65,66]. A caveat is that the nonperturbative effects, and hence the bounds on the fundamental parameters α and β inferred from them, strongly depend on the assumed equation of state of nuclear matter [66,67]. For this reason, we will adopt a model-independent perspective and search for constraints on the effective parameters α A and β A rather than on fundamental model parameters (cf. [68]). Finally, if nonperturbative effects are absent, one can assume the typical values given by Eqs. (27) and (28).
The effective Lagrangian describing the nonrelativistic motion of a binary system interacting with ULDM in Fermi coordinates has the form (23) with the masses promoted to the field-dependent expressions (25). In deriving the equations of motion we keep only terms that are first order in the sensitivities, leading order in the nonrelativistic approximation, and up to linear in the spatial gradient of Φ. Combining the equations for ⃗ ξ 1 and ⃗ ξ 2 and using expression (11) for the ULDM field we obtain in the case of linear coupling where the BB position is defined using the unperturbed masses, and Φ and _ Φ are evaluated at ⃗ ξ ¼ 0. We have introduced the following notations: 6 The sensitivities vanish in the extreme case of a black hole due to the no-hair theorems [27]. Thus, the asymmetry (28) is expected to be of order one in the case of a NS-black hole binary. 7 It has been further shown that in a binary NS-NS system a nontrivial scalar profile due to scalarization of one of the binary members can lead to induced scalarization of the other with the transition to a phase with nonzero scalar charge occurring dynamically [59,60].
Two comments are in order. First, we observe that the orbital dynamics and the BB motion are coupled whenever the interaction with the scalar field is nonuniversal, Δα ≠ 0. Second, ⃗ R ¼ 0 is no longer a solution of the equations of motion due to an extra (nonuniversal) force exerted by the ULDM field on the system as a whole. In other words, the BB does not move along a geodesic and cannot be chosen as the origin of the Fermi normal coordinates. The Fermi system constructed in the previous subsection and used to derive (23) should be understood as associated with the unperturbed motion of the BB obtained if one neglects all terms proportional to α A . The values of ⃗ R and _ ⃗ R induced in this frame are OðαÞ. Thus, the terms containing these quantities on the right-hand side (RHS) of (30) are Oðα 2 Þ and will be neglected in what follows. We will discuss below to what extent the modification of the BB trajectory at order OðαÞ implied by (30a) may be relevant for observations.
For the case of quadratic coupling the equations read where we have already neglected the terms proportional to ⃗ R, _ ⃗ R. The coefficients β T , β μ , and Δβ are defined similarly to (32). Note that these equations can be obtained from (30) by the substitution To simplify the notations, we will omit the overbar on the values of the masses at Φ ¼ 0 in the rest of the paper.
Equations (30) and (33) describe the perturbation of the Keplerian motion due to the interaction with the ULDM field. The orbital motion will also be affected by the GR effects, such as the periastron precession and emission of gravitational waves. As we mentioned above, as long as both GR and ULDM corrections are small, they sum up linearly and can be treated independently. Thus, we are not going to consider the GR corrections explicitly in our analysis, focusing on the new effects due to ULDM. When comparing our results to the data, we will use the residuals derived after subtracting the known GR contributions.

IV. THE PERTURBED KEPLERIAN PROBLEM
We now derive the modifications to Keplerian orbits from the terms in Eqs. (30b) or (33b) depending on Φ. We treat the new terms as perturbations and work in the framework of osculating Keplerian orbits [69,70]. Our calculation is analogous to the case of a gravitational wave perturbation of Ref. [71].
We start with unperturbed Keplerian motion and introduce Cartesian coordinate system ðx; y; zÞ centered on the body with mass M 2 with the axes x and y lying in the orbital plane and the unit vectorx pointing toward the pericenter. We also use cylindrical coordinates ðr; θ; zÞ, so that the position of the body M 1 is given by ⃗ r ¼ r cos θx þ r sin θŷ. The system ðx; y; zÞ is rotated with respect to the observer reference frame ðX; Y; ZÞ by the Euler angles Ω (longitude of the ascending node), ι (inclination), and ω (longitude of the pericenter); see Fig. 1. The motion within the orbital plane is characterized by the orbit's semimajor axis a, its eccentricity e, and the time of the pericenter passage t 0 ; see Appendix A for more details.
A perturbing force ⃗ F leads to deviations from the Kepler motion. The idea of the osculating orbits method is to find at every given instance a Keplerian orbit that touches the actual trajectory. The parameters of the latter orbit-its orbital elements-evolve with time according to the equations [71] (see Appendix A for an outline of the derivation), FIG. 1. Orientation of the binary system orbit in the fundamental reference frame ðX; Y; ZÞ. The orbital frame ðx; y; zÞ is centered at the body with mass M 2 . It has x and y axes in the plane of the orbit, with the x-axis pointing toward the pericenter. The angle between the axis X and the unit vectorn pointing toward the ascending node is denoted by Ω, the inclinations between the planes xy and XY are denoted by ι, and the angle betweenn and the x-axis is ω. θ is the angle between the direction toward the body M 1 and the x-axis. We also show the vectorb that will be used in Sec. IV B. The observer looks at the system from Z → −∞. where is the binary orbital frequency and we have decomposed the perturbing force in the cylindric coordinates, In Eqs. (34e) and (34f) we have introduced the argument of the pericenter, and the parameter related to the pericenter passage time t 0 . Working at leading order in ⃗ F we can evaluate the RHS of Eqs. (34) on the unperturbed Keplerian orbit. We will be interested in secular effects accumulating over many orbits.
Oscillations of ULDM will also produce nonsecular effects that can, in principle, be used to further probe the existence of Φ [38]. However, one expects nonsecular perturbations to be suppressed compared to the secular ones.
Equations (34) are to be supplemented by Eqs. (30a) or (33a) describing acceleration of the BB. We now consider separately the cases of purely gravitational interaction

A. Interaction only through gravity
We consider first the case of no direct coupling α A ¼ β A ¼ 0. As discussed in Sec. III A, in this case the BB follows a geodesic in the galactic gravitational field (including the Newtonian potential from DM).
For the orbital dynamics, we read out the components of the perturbing force from Eq. (22): The vanishing of F z implies that the orbital elements Ω and ι describing the orientation of the orbital plane with respect to the reference coordinate system will remain constant. This can be traced back to the absence of a preferred direction. The nontrivial orbital equations reduce to Note that Eqs. (39b) and (39a) imply the conservation of the combination að1 − e 2 Þ which is connected to the angular momentum L of the osculating orbit; see Eq. (A5a).
Using Eq. (19) and the Fourier decomposition of the variables that describe the unperturbed Keplerian motion [see Eqs. (A6) in Appendix A], we see that there is a secular effect when a resonant condition m Φ ≈ Nω b =2, N ∈ N is satisfied. More precisely, let us define which we will assume to be small, δω 2 ≪ 2m Φ . One computes the secular contributions into the temporal derivatives of the orbital parameters by averaging over time intervals Δt satisfying P b =N ≪ Δt ≪ 2π=δω 2 , where we have introduced the orbital period The result is h_ where J N ðxÞ is the Bessel function, J 0 N ðxÞ is its derivative with respect to the argument, and is the combination describing the phase difference between the ULDM-induced metric oscillations and the binary orbital motion. Expanding J N ðNeÞ and J 0 N ðNeÞ for small values of e, one observes that the effect on the orbital period increases with e and vanishes as e → 0. Therefore, we expect to obtain the strongest constraints from measurements of h _ P b i for highly eccentric systems. The pericenter advance h _ ϖi is, on the contrary, enhanced for orbits with low eccentricity, notably for N ¼ 1. However, measurements of h _ ϖi for nearly circular orbits suffer from large uncertainties. Therefore, mildly eccentric orbits are potentially the best candidates to constrain this effect.
The numerical estimates of the effects (42) and their comparison with the precision of pulsar timing measurements will be discussed in Sec. VA 1.

B. Direct coupling
We now turn on the direct coupling between ULDM and the binary. We will assume that it generates stronger effects than the purely gravitational interaction of the previous subsection, which will be ignored.

Motion of the binary barycenter
We start by analyzing the effects on the motion of the BB. Consider first the linear coupling case. At leading order in α A one can replace the dynamical variables on the RHS of (30a) by their unperturbed values. Then, using the Keplerian expression for relative acceleration of the binary, the BB equation can be cast in the form We are interested in the secular contribution that survives after averaging over the orbital motion and the oscillations of the field Φ. Only the last term produces such a contribution provided the system is close to the resonance. The latter corresponds to the condition, Using the expression (11) for the ULDM field and Eqs. (A6) for the Fourier decomposition of the Keplerian variables, one averages over time intervals Δt, such that P b =N ≪ Δt ≪ 2π=δω 1 . This yields where We observe that a nonuniversal coupling leads to a secular acceleration of the BB along a direction in the orbital plane.
If N ¼ 1, the effect survives for nearly circular orbits, We have introduced here the time-dependent unit vector (see Fig. 1) 8b ðtÞ ¼x cos γ 1 ðtÞ −ŷ sin γ 1 ðtÞ: ð51Þ Note that the secular BB acceleration is suppressed by the small frequency difference δω 1 . In particular, it vanishes exactly at resonance. As will be discussed in Sec. V, this makes the BB motion less sensitive to the ULDM effects than the orbital dynamics. For quadratic coupling one uses Eq. (33a). The square of the field amplitude at the binary system location entering the RHS can be written as The oscillating term leads to the same phenomenology as in the case of linear coupling. The resonance condition now coincides with that for pure gravitational interaction [see Eq. (40)], and the corresponding secular BB acceleration is obtained from (47) by the replacement 8 This can be equivalently written in the form, bðtÞ ¼n cosðγ 1 ðtÞ − ωÞ − ðẑ ×nÞ sinðγ 1 ðtÞ − ωÞ; ð50Þ which shows that for N ¼ 1 resonancesbðtÞ is regular in the limit of strictly circular orbits. In this limit the pericenter passage time t 0 entering into the definition of γ 1 ðtÞ, as well as the pericenter longitude ω, are not separately well-defined. On the other hand, the combination ω − ω b t 0 is regular and equals the angle at t ¼ 0 between the radial vector ⃗ r and the direction toward the ascending node of the orbit. Thus, the quantities that remain regular for circular orbits can depend only on the difference γ 1 ðtÞ − Nω. This is indeed the case of Eq. (50) An important difference is introduced by the constant term in (52), which gives rise to a secular BB acceleration irrespective of whether the system is in resonance or not, In particular, the acceleration will be present even for the ULDM mass m Φ well above the orbital frequency. In fact, it will be present even for solitary pulsars and stars. 9 The effect has a clear physical interpretation. Recall that the masses of the binary constituents depend on the local DM density, whereas ⃗ S points toward decreasing DM density. For β T > 0 (β T < 0) the system tends to minimize its mass by moving toward a region of lower (higher) ρ DM .

Evolution of the orbital elements: Linear coupling
We first derive the orbital equations in the linear coupling case and then discuss the modifications for quadratic coupling. From Eq. (30b) we obtain the components of the perturbing force, These are to be inserted into Eqs. (34). The result is quite lengthy and not particularly illuminating; for completeness we give it in Appendix B. What we are interested in are the secular contributions which appear when the system satisfies the resonant condition (46). Averaging over multiple orbits we obtain the following equations: Here we have introduced the time-dependent vectors, In deriving (56) we used the Fourier decomposition of the Keplerian motion; see Eqs. (A6). Note that for j ⃗ Sj ∼ m Φ j ⃗ Vj, which is typical of the ULDM halo, the magnitude of ⃗ U, ⃗ U 0 is given by the value of j ⃗ Vj in units of the orbital velocity, As discussed in Sec. II, we expect j ⃗ Vj ∼ 10 −3 , and hence the terms multiplying Δα in (56) can be as large as Oð10Þ for binaries with long orbital periods.
Equations (56) describe the new effects due to the ULDM coupling. To obtain the full evolution of the orbital parameters, one has to add the standard contributions due to the deviations from Keplerian dynamics, notably, the general relativistic corrections. Assuming that all corrections are small, they can be combined linearly.
In what follows we focus on two representative cases: Case 1: Δα ¼ 0 or negligible ULDM gradients ∇Φ ¼ 0, eccentric binaries. The ULDM interaction in this case does not introduce any preferred direction, so the orbital elements Ω and ι describing the orientation of the orbital plane remain constant. The effect on the orbital period P b ∝ a 3=2 in this case was already discussed in [35]. Similarly to what happens for pure gravitational interaction, it vanishes for circular orbits, e ¼ 0. The use of the osculating method allows us to derive how all the orbital parameters are modified in the presence of Φ. We postpone the numerical estimates of the new effects and the discussion of resulting bounds on the ULDM couplings till Sec. VA 2. Case 2: Δα ≠ 0 and ∇Φ ≠ 0, nearly circular orbits. This case is important because many binary pulsars are in systems with low eccentricities [43,44]. By inspection of Eqs. (56) we see that for e ≪ 1 nontrivial effects survive only for the resonances on the principal (N ¼ 1) and the second (N ¼ 2) harmonics. We consider them in turn. For N ¼ 1, e → 0, the secular equations simplify to wherebðtÞ is the unit vector introduced in (51) and dots in (59e) stand for terms of order Oð1Þ at e → 0. Note that, unlike the case of universal coupling, we now have a nontrivial effect on the semimajor axis and hence on the orbital period P b . Besides, the orbital plane will be rotating whenever the ULDM velocity and/or density gradient have nonvanishing projections on its normal. A comment is in order. As pointed out in footnote 8, the quantity γ 1 is ill-defined for circular orbits and should appear in the equations for physical observables only in the combination γ 1 − Nω. Recalling the expression (50) forb we see that this requirement is satisfied for all Eqs. (59), except those for the eccentricity and the argument of the pericenter, Eqs. (59b) and (59e). The resolution of this apparent inconsistency lies in the observation that the perturbation theory developed in terms of e and ϖ breaks down when e → 0. Instead, a pair of good variables is provided by [41] κ ¼ e cos ω; η which are the components of the Runge-Lenz vector along the directionsn and ðẑ ×nÞ, respectively (see Appendix A).
To obtain the secular evolution of the new observables we can use Eqs. (59b) and (59e) at small, but finite e and then take the limit e → 0. This yields the equations which depend on γ 1 − ω, as they should. which replacesb for N ¼ 2. We observe that if ⃗ V and/or ⃗ S have nonvanishing projections on the orbital plane, the orbit gets polarized. This is reminiscent of the Nordtvedt effect [72] in SEP violating theories that was studied in the context of binary pulsars by Damour and Schäfer [73]. In the latter case one considers an approximately constant extra force ⃗ F DS ¼ Δ⃗ g generated by the differential acceleration of the two binary members in the galactic gravitational field. For weakly circular orbits the resulting secular equations are with no change in the other orbital elements. One may wonder why the effects of the constant (in the case of Damour-Schäfer) and oscillating (in our case) forces turn out to be similar. To understand this one recalls that what matters for the osculating orbit equations (34) are the components of the force in the frame ðr;θÞ attached to the orbiting bodies. A constant force will rotate in this frame clockwise (fromθ tor) with frequency ω b . This resonates with cos θ and sin θ factors in Eqs. (34b) and (34e) to produce the secular drift of the eccentricity parameters (64). Consider now the N ¼ 2 ULDM resonance. The perturbing force oscillates in the fixed reference frame with nearly twice the orbital frequency m Φ ≈ 2ω b . Its components in the rotating frame ðr;θÞ contain harmonics with frequencies m Φ þ ω b ≈ 3ω b and m Φ − ω b ≈ ω b . The latter harmonic corresponds to a contribution that rotates in the plane ðr;θÞ counterclockwise and resonates with the θ-dependent factors in Eqs. (34b) and (34e). Apart from the flip of the angular velocity of the force in the ðr;θÞ plane, which only affects the overall coefficient in the secular equations, the similarity with the Damour-Schäfer case is clear.

Evolution of the orbital elements: Quadratic coupling
In the case of quadratic coupling we encounter two types of effects: resonant and nonresonant ones. The former are due to the oscillating contribution in Φ 2 =2 [see Eq. (52)] and appear if m Φ is close to Nω b =2. They are the same as for linear coupling. The corresponding secular equations are obtained from (56) by the substitution (53) supplemented with Then, the rest of the discussion of the previous subsection applies without change.
The nonresonant effects come from the first term in (52). From Eq. (33b) we see that it leads to the perturbing force in the orbital equation, where the superscript "nr" stands for nonresonant. The first contribution has the same form as the Newtonian force and can be absorbed into the redefinition of the total mass. We will drop it in what follows. The second term is present only for nonuniversal coupling and represents an approximately constant 10 force pointing along the gradient of ULDM density. This is exactly the situation considered by Damour and Schäfer [73], with the only difference that, unlike their case, the differential acceleration need not point toward the galactic center. The resulting phenomenology is wellknown. Being conservative, the new force does not affect the orbital period. However, it induces an evolution of the eccentricity and of the orientation of the orbit. The full set of secular equations reads In the limit of circular orbits only h_ ei nr and h _ ϖi nr survive and combine into equations of the form (64) for the parameters κ and η with the replacement ⃗ F DS ↦ ðρ DM =m 2 Φ ÞΔβ ⃗ S. Due to nonresonant effects, timing even of a single binary pulsar allows us to probe a wide range of ULDM masses m Φ ≪ T −1 obs V −2 0 , where T obs is the observation time 10 It varies on timescales of the ULDM coherence time which we assume to be much longer than the orbital period. and V 0 is the DM virial velocity. The upper limit comes from the requirement that T obs is less than the ULDM coherence time t coh [see Eq. (10)], so that the ULDM density gradient does not change during the observational campaign. This requirement is implied by the standard pulsar timing procedure which assumes that the time derivatives of the orbital elements can be treated as constant. In principle, it can be relaxed if the analysis is appropriately modified to allow for time variation of ⃗ S (see more on this below).

V. NUMERICAL ESTIMATES
In this section we estimate the size of the effects derived above and compare it with the precision achieved in pulsar timing measurements. We do not aim at a systematic exploration of the bounds set by observations on the ULDM properties which is beyond the scope of this paper. Our goal is to illustrate the potential reach of pulsar timing constraints and identify the most promising directions for further analysis.
The standard timing models for binary pulsars parametrize the time of arrival of the ith pulse as a function of i, the mass of the pulsar M 1 , its companion M 2 , and the set of parameters fXg which include the pulsar spin frequency, orbital elements, BB velocity, and acceleration with respect to the Solar System, etc. [39][40][41]. Secular variations are introduced phenomenologically by allowing a given parameter X to be a linear function of time, XðtÞ ¼ X 0 þ _ Xt. In our context, this approximation is valid if the system is sufficiently close to resonance (for resonant effects) or if the observation time does not exceed the ULDM coherence time (for nonresonant effects). Violation of these conditions leads to a loss of sensitivity that we will qualitatively discuss later in this section. A systematic study should incorporate the modifications we derived, including modulations due to deviations from resonance or loss of coherence of the ULDM field, in the timing package [41]. We leave this study for the future.
Within the assumption of a constant secular drift of the orbital parameters Ref. [39] estimated the statistical uncertainty of the measured drift in a given campaign. 11 We quote the results for the time derivatives of the orbital period and eccentricity which are typically the two quantities measured with the highest precision, 12 where the coefficients C P and C e depend on the eccentricity and the orientation of the orbit and are typically of order one, _ n dat is the rate at which the timing data are acquired with the precision δt, and T obs is the total time of observation. Here "d" stands for days and "yr" for years. Equations (68) set the benchmark for the magnitude of the effects that can potentially be detected. It is worth noting that the numbers we are using are rather conservative. The future square kilometer array (SKA) is expected to reach the timing precision of δt ∼ 80-230 ns with 10 min of integration of SKA [74], which suggests that the data acquisition rate _ n dat can exceed 1=d. We observe from Eqs. (68) that the statistical uncertainty in determination of _ P b decreases with observation time faster than that of _ e. However, the measured value _ P meas b differs from the intrinsic one _ P int b due to a number of systematic effects. The most important among them are related to the motion of the BB with respect to the Solar System barycenter: an apparent change of the orbital period due to the galactic acceleration [75] and the Shklovskii effect [76]. The corresponding contributions have to be subtracted to find _ P int b . This introduces additional uncertainty in the comparison of the theoretical predictions made in terms of _ P int b with the data. On the other hand, the measurements of _ e are not contaminated by any known external systematic effect [77]. This gives them an advantage over measurements of _ P b .

A. Resonant effects
We start with the resonant effects and assume that the deviations from the resonance is small enough so that the period of secular modulations is significantly longer than the duration of the observational campaign. More precisely, in the numerical estimates we will assume jδω 1;2 j ≲ 1=T obs . Then the modulation phase γ 1;2 in the secular equations can be treated as constant which corresponds to a constant drift of the orbital parameters. 11 Reference [39] considers the binary pulsar PSR1913 + 16 and makes some simplifications specific to this system in the derivation of the covariance matrix for the orbital parameters. Still, the final results can be applied to other systems as orderof-magnitude estimates. 12 Another well-measured quantity is the precession of the pericenter _ ω. This, however, receives large GR contributions [54] and is used to extract the total mass of the binary M T . Thus, it does not provide independent constraints on deviations from GR.

Pure gravitational coupling
This case was considered in Sec. IVA. Replacing the semimajor axis a in favor of the orbital period P b using expression (41) and substituting realistic values into Eqs. (42) we obtain the estimates As already pointed out in [35], detection of these effects presents a strong challenge. It would require a determination of h _ P b i with the accuracy of 10 −17 for systems with orbital periods of ∼100 days. This accuracy may be achieved in the future for the double pulsar [78]. The latter, however, has P b ∼ 0.1 d, which makes the effects described by Eqs. (70) completely negligible. Equation (68a) shows how the precision deteriorates for the nonrelativistic systems where the effect of ULDM is more sizable.
The consideration of the rest of the orbital parameters shown in (70) does not change the previous conclusion.
Hence, unless a binary system is found in a region of large DM density, the prospects to detect the ULDM-induced oscillations of the galactic gravitational field with binary pulsars are slim.

Universal direct coupling
We now show that timing of binary pulsars can give relevant bounds on direct ULDM-SM couplings. Let us first neglect the difference in the effective couplings of the two binary members. Setting For the case of quadratic coupling substitution (65) yields Comparison of the above expressions with the estimated precision (68) suggests that pulsar timing can probe the values of the linear (quadratic) coupling down to 13 α ∼ 10 −23 GeV −1 (β ∼ 10 −32 GeV −2 ). The range of ULDM masses accessible at the maximal sensitivity is restricted by the resonant condition, which reads obs ) for the Nth resonance. In systems with large eccentricity several first resonances have comparable constraining power.
As an example, let us consider the binary pulsar J1903 þ 0327. The system consists of a millisecond pulsar with mass M 1 ≃ 1.67 M ⊙ and a main-sequence companion star with M 2 ∼ M ⊙ . The orbital period and eccentricity are P b ¼ 95 d, e ¼ 0.44. For this system both _ P b and _ e have been measured using the data collected over T obs ∼ 3 yr [79]: Taking conservatively the measured value of _ P b as the upper limit on the possible effects of ULDM we obtain the constraints shown in Table I. The constraints from the drift of the eccentricity are a factor of 20 weaker.
The numbers in the third column of the table can be compared with the bounds on linearly coupled fields following from the measurements of light deflection by the Sun [33] and planetary dynamics [32]. These are at the level α ≲ 10 −21 GeV −1 in the considered mass range. We see that the pulsar bounds are competitive. Moreover, as emphasized in Sec. III, they probe effective scalar couplings in the regime of strong gravity, which are in general different from the weak-field coupling tested by the Solar System observations. For the quadratically coupled ULDM the existing bounds in this mass range are very mild (see the discussion in [35]). This makes binary pulsar timing particularly promising for testing this type of models. As highlighted in [77], the determination of the orbital parameters of J1903 þ 0327 is expected to improve in the future. Apart from the sheer increase of statistics due to a longer campaign, the improvement will come from infrared interferometric observations of the main sequence star companion, allowing one to better determine the location of the system in the Galaxy and the orientation of its orbit.
For systems with nearly circular orbits a universal interaction with ULDM does not lead to any change in the orbital period. Then the strongest effect is the drift of eccentricity occurring at the N ¼ 1 resonance. To illustrate the ensuing constraints, we consider the system J1713 þ 0747 that has been timed for more than two decades, T obs ∼ 20 y [80][81][82]. It consists of a millisecond pulsar with mass M 1 ≃ 1.3 M ⊙ in a wide nearly circular orbit (P b ≃ 67.8 d, e ≃ 7 × 10 −5 ) with a white dwarf companion of mass M 2 ≃ 0.29 M ⊙ . The residual drifts of its orbital period and eccentricity parameters (after subtracting the known standard contributions) have been constrained at the level [82]: Comparing the two latter bounds with the predicted effect of ULDM, Eqs. (61), and analogous equations for the quadratic coupling, we obtain the constraints presented in Table II. These bounds are comparable to those in Table I, but apply to different resonant ULDM masses, as a consequence of different orbital periods of J1713 þ 0747 and J1903 þ 0327. We stress that they come entirely from the measurement of the eccentricity parameters. Note that, unlike the bounds in Table I, these constraints are independent of the relative phase γ 1;2 between the orbital motion and the ULDM oscillations. Indeed, the expressions (61) for the drift of the eccentricity parameters κ and η involve both sin γ 1;2 and cos γ 1;2 . Hence, by taking the combination ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h_ κi 2 þ h_ ηi 2 p the dependence on γ 1;2 can be eliminated.  This under the assumption that the phase γ 1;2 takes a random value between 0 and 2π, so that sin γ 1;2 is a number of order one.

Nonuniversal coupling
All orbital elements of a binary system will in general drift under the influence of ULDM, once the nonuniversality of the coupling is taken into account. In our preliminary analysis we focus on the orbital period and eccentricity. We also restrict ourselves to systems with nearly circular orbits, which constitute the majority of observed pulsar binaries [43,44]. As derived in Sec. IV B 2, only resonances on the first two harmonics survive in this case. The precise magnitude of the effects depends on the directions of the ULDM velocity and density gradient with respect to the orbital plane. It is convenient to introduce the following notations: and V is the absolute value of the relative velocity between the ULDM and the binary. Generically, one expects the quantities q 1 ,q 1;κ ,q 1;η , q to be Oð1Þ numbers. Note that q depends only on the projections of ⃗ V and ⃗ S on the orbital plane and is independent of the phase difference between the ULDM oscillations and the orbital motion. Substituting the numerical values into Eqs. (59a) and (62) and using the definitions of the vectorsb andb from (50) and (63) we obtain The expression for h_ ηi is the same as (75b) with the replacementq 1;κ ↦q 1;η . Using the substitution (53) one obtains the equations for the case of quadratic coupling, h_ κi ≃ 11.9 × 10 −20 s −1 Δβ and similarly for h_ ηi withq 2;κ changed toq 2;η . The coefficients q 2 ,q 2;κ ,q 2;η have the same form as Eqs. (73) with γ 1 replaced by γ 2 .
To illustrate the attainable sensitivity, let us again consider the system J1713 þ 0747. The constraints on the drift of its orbital period and eccentricity yield the limits on the ULDM couplings shown in Table III. The limits set by _ P b are the strongest current constraints in the relevant mass ranges, for both the quadratic and the linear couplings. In particular, they are by a few orders of magnitude stronger than the bounds in Table II obtained for the same system assuming universal couplings. On the other hand, the constraints in Table III do not restrict the values of the ULDM couplings directly, but only their difference for the two members of the binary. As discussed TABLE III. Limits on the difference of the ULDM couplings to the members of the binary pulsar system J1713 þ 0747. For masses corresponding to the N ¼ 1 resonance the limits come from the measurement of _ P b , whereas the N ¼ 2 resonances are constrained by _ κ and _ η. We assume the ULDM density at the system location ρ DM ¼ 0.3 GeV=cm 3 and the relative velocity between the binary and ULDM V ¼ 10 −3 .  in Sec. III B, in general one expects Δα ∼ 0.1α and similarly for β, but the precise relation may depend on the model. Besides, the bounds in Table III depend on unknown order-one coefficients q 1;2 , so inferring constraints on ULDM couplings from them requires an additional assumption that these coefficients do not vanish accidentally. In this respect, the bounds in Table II are more robust. The bounds from _ κ and _ η listed in Table III are significantly weaker compared to those from _ P b , but they probe different masses. In general, constraining ULDM couplings from eccentricity measurements can have certain advantages. First, as already mentioned, such measurements are less subject to systematic uncertainties than those of _ P b . Second, the ULDM contribution into _ κ, _ η has one less power of ð P b 100 d Þ which makes it more important than _ P b for not so slow systems. Yet more powerful constraints can potentially be derived from a dedicated analysis of timing data using the complete set of equations for the orbital parameters.
The analogy with the Nordtvedt effect mentioned at the end of Sec. IV B 2 may suggest to consider the measurements of the absolute value of eccentricity, rather that its time derivative, as a possible probe of ULDM along the lines of the Damour-Schäfer test of SEP violation [73]. Notice, however, that the standard Damour-Schäfer test assumes the rate of periastron advance _ ω PN , given by the post-Newtonian GR corrections, to be sufficiently large so that the projection of the galactic acceleration ⃗ g onto the orbit can be considered constant over timescales 2π= _ ω PN . The corresponding assumption in our case would be that the RHS of Eqs. (62) can be regarded as approximately constant over the same timescales implying the condition _ γ 1;2 ðtÞ ¼ δω 1;2 ≪ _ ω PN . In other words, the system must be extremely close to resonance, which is a highly improbable situation. Thus, the Damour-Schäfer-type analysis does not appear promising for the study of resonant ULDM effects. Its relevance for probing nonresonant effects will be discussed in Sec. V C.
Finally, let us discuss the resonant BB acceleration which appears in the case of nonuniversal coupling and which we have disregarded so far. Using Eq. (49) for circular orbits we obtain for the N ¼ 1 resonance, Notice that we have normalized the frequency of secular modulations δω 1 to the value 0.1 yr −1 which ensures that the acceleration stays constant over an observational campaign of T obs ∼ 10 yr. To estimate the experimental precision, with whicḧR can be measured, we use the following reasoning. An acceleration along the line of sight contributes through the time dependent Doppler effect into the apparent change of the pulsar spin frequency _ ν meas =ν ∼̈R k . The latter can be measured in a given campaign with the uncertainty [39] δ_ ν ν This provides us with a proxy for̈R statistical uncertainty. We see that for the interesting values of Δα ∼ 10 −23 GeV −1 it is almost 4 orders of magnitude bigger than the expected effect. A similar result holds for the case of quadratic coupling. Moreover, using (78) as an estimate of δ̈R is certainly overly optimistic. The ULDM contribution must be disentangled from other sources of _ ν meas , such as the intrinsic spin-down, the galactic acceleration, and the Shklovskii effect. In principle, this might be possible by observing the secular modulations. However, for a fixed campaign duration T obs , observing modulations leads to a further increase of the statistical error (see below). We conclude that the BB acceleration has less constraining power than the study of the orbital elements and can be safely neglected.

Estimate of backreaction on ULDM distribution
In all previous calculations we have assumed that the ULDM configuration is not modified by the interaction with the binary system. The conditions for this approximation are formulated in (13). Let us verify them. The average rate of the momentum transfer from ULDM to the binary is given by the acceleration of its barycenter, _ p b ¼ M T ḧRi, and the first condition in (13) becomes On the other hand, the energy transfer affects primarily the energy of the Keplerian motion. The latter is related to the binary semimajor axis by the standard expression, Connecting the latter to the orbital period, we obtain that the second condition in (13) amounts to We see that both conditions are safely fulfilled in the examples considered above.

B. Detuning from resonance
We now consider the case when the modulation period T mod defined in Eq. (69) is comparable to or shorter than the observation time T obs , but still longer than the binary period. More precisely, we are interested in the behavior of the system at timescales t, such that Therefore, secular effects can still be isolated by averaging over a time interval Δt encompassing several orbital periods around t, but shorter than T mod , The new feature is that averaged quantities obtained in this way will now depend on time due to the detuning from resonance. We want to understand how the timing sensitivity changes in this situation. For simplicity, we focus on the variation of the orbital period in the case of universal linearly coupling. Keeping in the expression for h _ P b i only the contribution of the resonant mode N and integrating this expression from 0 to time t, we obtain the orbital period as a function of time, where ΔP b ðtÞ ≡ P b ðtÞ − P b ð0Þ. The behavior of ΔP b during an observational campaign with T obs ∼ 20 y is illustrated in Fig. 2 for three choices of δω 1 . We observe that the amplitude of the modulations decreases inversely proportional to δω 1 . Thus, although the modulations provide a characteristic signature of ULDM, they become harder to detect for a given precision of the P b measurements. The extraction of the sensitivity to this signal from observations requires the reanalysis of the data with a modification of the data analysis package (e.g., TEMPO2 [41]). We leave this for future work. Here we provide a rough estimate using the following strategy. We chop the oscillating signal into intervals of almost resonant behavior, bound the secular drift of the orbital parameters using data from each interval, and combine the measurements from different intervals assuming that their errors are uncorrelated. In other words, we consider the signal as a triangle wave, starting at an initial time T 0 and with _ P b constant over typical intervals of duration T I ¼ T mod =4. For simplicity, we assume that T 0 coincides with the first point of the measurement. This can be achieved by scanning the data for different T 0 . The error on _ P b obtained by fitting the data in an interval T I is read off from (68a), where ðδ _ P b Þ T obs is the would-be uncertainty for the case when _ P b stays constant for the whole observational campaign. The total sensitivity is then estimated by dividing ðδ _ P b Þ I by the square root of the number of intervals, where in the final step we used (69). This formula applies when δω 1 ≥ π=2T obs ; otherwise, one should just use ðδ _ P b Þ T obs . We see that the uncertainty in the determination of _ P b grows as the system deviates more from resonance, which degrades the sensitivity to ULDM couplings.
To illustrate this point, we consider the binary pulsar J1903 þ 0327 which we already discussed in Sec. VA 2. Current measurements over T obs ∼ 3 yr constrain the drift of its orbital period at the level ðδ _ PÞ T obs ≲ 10 −10 . Using this value we obtained the estimates of the constraints on the FIG. 2. The relative change of the orbital period due to the ULDM coupling as a function of time. The three lines correspond to cases differing by the detuning from resonance: exact resonance, δω 1 ¼ 0 (black solid line); mild detuning, δω 1 ¼ π=2T obs (red dashed line); strong detuning, δω 1 ¼ 8π=T obs (blue dotted line). We consider the N ¼ 1 resonance and take the initial resonant phase γ 1 ð0Þ ¼ π=2. Other parameters are e ¼ 0.44, ULDM coupling at resonances listed in Table I. In Fig. 3 we show with black lines how the constraints degrade when the ULDM mass m Φ is away from an exact resonance. The curves are obtained by taking Eq. (86) for the upper limit on the amplitude of h _ P b i, which is admittedly a rather crude estimate. Thus, the bounds in Fig. 3 should be interpreted as indicative. A robust exclusion requires a dedicated analysis of J1903 þ 0327 timing data, which is outside the scope of this paper. We have restricted our analysis to the vicinity of resonances where a single harmonic dominates. Thus, the bounds apply within resonant bands of width δω 1 ≪ ω b ; they are shown by solid lines. In principle, one can extend the calculation beyond the resonant case by summing up the contributions from multiple harmonics or solving the orbital motion numerically [38]. We do not perform such a calculation in this paper. Instead, to get a rough idea of the expected result, we extrapolate the resonant expressions to δω 1 ∼ ω b ; this extrapolation is shown with dashed lines.
For comparison, we also show how the constraints can be improved assuming J1903 þ 0327 or a system with similar parameters is timed for T obs ∼ 20 yr and the drift of its orbital period is constrained at the level ðδ _ P b Þ T obs ≲ 10 −13 (blue curves). The sensitivity to the parameter α precisely at resonances (at the dips of the V-shaped curves in Fig. 3) scales with the observation time as T −5=2 obs . In contrast, it scales as T −5=2 obs · T 2 obs ¼ T −1=2 obs , outside of the resonance; cf. Eq. (86). Thus, the V-shaped exclusion regions become deeper and have sharper tips as the observation time increases. The plot suggests that future data even for a single system will be able to cover a substantial part of the ULDM parameter space. This, together with the possibility of detecting many systems with a variety of orbital parameters, motivates a more dedicated study.

C. Nonresonant effects for quadratic coupling
As we discussed previously, in the case of the quadratic coupling the ULDM force acting on the bodies has a nonoscillating contribution proportional to the DM density gradients. The latter are set by the de Broglie wavelength of ULDM in the halo. Thus, the vector where V 0 is the DM virial velocity, is expected to have order-one norm. The nonresonant force affects both the motion of the BB, Eq. (54), and the orbital elements, Eq. (67). Let us start with the first effect. Substituting the numerical values we obtain According to the discussion at the end of Sec. VA 3, detecting this small extra acceleration on top of the standard contributions will be extremely challenging. Still, it is worth stressing that this effect exists not only for binaries, independently of their orbital frequency, but even for solitary pulsars. Hence, one could think of increasing the sensitivity by a joint analysis of many systems. We do not pursue this possibility here. The study of the orbital parameters appears more promising. Here the leading effect is the drift of the eccentricity given by Eq. (67b), which upon inserting the numbers takes the form For nearly circular orbits this expression should be substituted by the equations for _ κ and _ η, which differ only by the replacement s y ↦ ⃗ s · ðẑ ×nÞ and s y ↦ −⃗ s ·n, respectively. Consider, for example, the system J1713 þ 0747. From the constraints on the drift of its eccentricity parameters we obtain the bound where we have assumed the typical values for the DM density and virial velocity, Importantly, this bound applies in a wide range of ULDM masses, where it is comparable to or stronger than the other existing bounds (discussed e.g., in Ref. [35]). The limitation comes from the standard assumption, made in inferring the bounds on the eccentricity drift, that the time derivatives of the orbital elements are constant over the observation time, which for J1713 þ 0747 is T obs ≃ 20 yr. This requires the constancy of the ULDM gradients on the same timescales.
In other words, the ULDM coherence time t coh should be longer than T obs . Recalling the expression (10) one obtains the condition (91).
In principle, one could extend the reach of pulsar constraints to higher masses by generalizing the timing analysis to nonconstant parameter drift along the lines of Sec. V B. However, there is one more restriction that prevents one from significantly extending the accessible parameter space, at least for the case of J1713 þ 0747. This is a non-negligible backreaction of the binary on the ULDM configuration at high m Φ . Indeed, substituting the expression for the barycenter acceleration (88) into the condition (79) we find that the backreaction is small provided that Assuming Δβ ∼ 0.1β T , we see that for the case of J1713 þ 0747 neglecting the backreaction in the derivation of the bound (90) is indeed justified for DM masses satisfying the condition (91). However, for higher masses one should include the backreaction, which is beyond the scope of this paper.
It is worth comparing the constraint (90) to the bounds on Δβ derived from the consideration of resonant effects in J1713 þ 0747 listed in Table III. We observe that the resonant drift of the orbital period leads to a stronger bound within the N ¼ 1 resonance band. However, for all other masses the broadband constraint (90) dominates. In particular, the N ¼ 2 resonant contribution to the eccentricity drift turns out to be smaller (by about a factor of 7) than the nonresonant one, as it is apparent from comparison of Eqs. (76b) and (89).
Observation of binary pulsars with nearly circular orbits yields yet another type of constraints. The interaction with ULDM leads to the growth of eccentricity, even if originally it was exactly zero. The growth continues until the ULDM force can no longer be considered as constant after the coherence time t coh , or until the orbit significantly changes the orientation due to the post-Newtonian periastron precession _ ω PN [73]. Thus, the typical increment of the eccentricity induced by ULDM is where h_ ei is given by Eq. (89), t coh is given by Eq. (10), and the precession time reads in GR, Imposing that Δe be smaller than the measured value of the eccentricity for a given system yields constraints on the ULDM coupling.
As an example we consider the system J0024-7201X [44,83] which has the orbital period P b ≃ 10.9 d and very low measured eccentricity e ¼ ð4.8 AE 1.5Þ × 10 −7 [84]. It consists of a pulsar with mass M 1 ≃ 1.4 M ⊙ and a companion with M 2 ≃ 0.4 M ⊙ ; the timing data were collected over a period T obs ≃ 15 yr. We find that for this system _ ω PN −1 is shorter than t coh for m Φ ≲ 6.2 × 10 −21 eV. Therefore, we impose h_ ei _ ω −1 PN < e for ULDM masses below this value and h_ eit coh < e for higher masses. This gives the bounds jΔβ⃗ sj < where we have assumed our usual reference values for ρ DM and V 0 . The accessible ULDM mass range is limited from above by the requirement of negligible backreaction. One observes that for m Φ ≲ 2 × 10 −20 eV, the bound (95) is somewhat stronger than the one obtained from a direct timing measurement, Eq. (90). It is, however, less robust (see e.g., [77]). Indeed, one cannot exclude a cancellation between the change in the eccentricity induced by ULDM and its initial value, so that the net eccentricity happens to be small during the observational epoch. This caveat can be removed by considering a population of low-eccentricity binary pulsars and performing a statistical analysis along the lines of [85][86][87]. We leave the exploration of this promising direction for the future.

VI. CONCLUSIONS AND OUTLOOK
In this paper we have made a detailed study of the influence of ULDM candidates on the evolution of binary systems and have identified potentially observable signatures in high-precision timing measurements of binary pulsars. We have focused on secular effects that accumulate over many orbital periods. We have shown that such effects appear due to the time evolution of the DM field Φ when the frequency of the Φ oscillations is in resonance with the orbital motion. Another class on (nonresonant) secular effects arises when Φ couples quadratically to the masses of the binary system members.
The resonant effects appear both due to the oscillations of the curvature of the spacetime generated by the DM stress energy tensor and from the possible direct coupling of the oscillating DM field to the SM fields. The first effect is generic for any ULDM candidate as it is based purely on the gravitational interaction of DM with the stars. However, from our estimate of the resulting drift of the orbital elements of binary pulsars [Eqs. (70)] we have concluded that it is too small to be constrained by current data. The prospects for reaching the required sensitivity in the near future are not very optimistic. Still, nature has been quite indulgent so far in providing convenient binary pulsar systems to test fundamental laws. Hopefully, the new generation of radio surveys, with SKA as the flagship [88], will also benefit from this generosity regarding DM.
We have performed a comprehensive analysis of the scenario where ULDM couples directly to SM through the dilatonic portal (linear or quadratic). We have argued that even if the interaction is universal (WEP preserving) at the fundamental level, the effective couplings to the binary system members will be nonuniversal due to violation of SEP. We have derived the complete set of equations for the orbital parameters in this scenario, both for the resonant [Eqs. (56)] and for the nonresonant [Eqs. (67)] cases. We have identified the two most prominent effects. The first one is the time evolution of the orbital period P b in the resonant case. In contrast to the previous studies that considered only universally coupled ULDM, we have found that this effect does not vanish for circular orbit due to the nonuniversal coupling. The second effect is the change of the orbital eccentricity e that appears in the resonant, as well as nonresonant, cases. It is similar to the polarization of orbits obtained in Refs. [36,37] for the scenarios with vector and tensor DM. We have discussed its relationship to the Damour-Schäfer effect in SEP violating theories [73].
We have estimated the sensitivity of the timing measurements to the above effects. For this purpose, we used the phenomenological pulsar timing model developed in [39][40][41]. In this model, the orbital parameters are allowed to change during the campaign at a constant rate; the time derivatives of different parameters are treated as independent variables. We have found that inside the resonant bands the strongest constraints come from the drift of the orbital period, provided the systematic error in the determination of _ P b is kept under control. One of the main foregrounds here is an apparent change of P b due to the acceleration of the binary with respect to the Solar System barycenter [75]. A better determination of the galactic acceleration in the vicinity of the binary (e.g., by tracking of nearby stars) would help to reduce this systematics. We have estimated the dependence of the constraints on the deviation of the ULDM mass m Φ from the resonance value (see Fig. 3). Despite a significant loss of sensitivity outside the resonant bands, we proposed that the combination of timing data from multiple binaries with different orbital periods may efficiently constrain wide regions in the parameter space of the ULDM candidates.
Eccentricity measurements are more robust against foreground contamination [77]. They provide complementary information that potentially can help to disentangle the ULDM effects from possible systematics. Moreover, they strongly constrain the nonresonant eccentricity drift present in the case of quadratically coupled Φ. These nonresonant constraints have an important advantage: the timing data even for a single binary pulsar allow one to put bounds on the quadratic coupling in a broad range of ULDM masses, m Φ ≲ 10 −18 eV. The accessible mass interval is restricted only by the requirement that the backreaction of the binary system on the field Φ should be negligible. All in all, we have concluded that pulsar timing can probe an interesting portion of the ULDM parameter space, yet unconstrained by other observations.
We have also pointed out that the direct coupling leads to an additional secular acceleration of the binary barycenter which, however, appears to be unobservable.
There are several future directions to explore. An immediate next step would be the incorporation of the effects derived in this work into the timing package TEMPO2 [41]. This will allow us to explore how the sensitivity can be improved by taking into account the correlated drift of all orbital elements and put rigorous constraints outside the resonant bands. Another promising direction is the statistical analysis of a population of loweccentricity binaries in order to constrain the nonresonant effects of quadratically coupled ULDM along the lines of the Damour-Schäfer test [73,[85][86][87]. It would also be interesting to explore whether ULDM may help to explain the gaps in the P b − e plane in the population of binary pulsars [89,90].
We have focused on binary pulsars with nonrelativistic orbits. The reason is not only simplicity, but the fact that most of the effects we discussed become weaker as the orbital period decreases. However, given that many fast binaries present excellent timing opportunities, it will be worth extending our calculations to include post-Newtonian corrections. Furthermore, the interaction of a binary system with Φ will lead to the emission of dipolar radiation which is tightly constrained [91]. It is important to understand the implications of this process in the present context. Finally, our analysis has been formulated in terms of effective ULDM couplings to the masses of the binary system members. It will be interesting to relate these effective couplings to the fundamental parameters in specific models, in particular when nonperturbative scalarization effects may show up.
Before closing, let us notice that all our estimates assume the value ρ DM ¼ 0.3 GeV=cm 3 corresponding to the expectations in the neighborhood of the Solar System. An enhancement of this number would generate stronger constraints. In particular, the discovery of binary pulsars in the center of the galaxy would boost all our bounds and estimates.
The eccentric anomaly E is implicitly defined as a function of time by the equation where ω b is the orbital frequency given by Eq. (35). The eccentric anomaly changes by 2π over one orbital period. The motion is completely characterized by six constant parameters (orbital elements): semimajor axis a, eccentricity e, the longitude of the ascending node Ω, the inclination angle ι, the argument of the pericenter ω (not to be confused with the orbital frequency), and the time of the pericenter passage t 0 .
The system possesses two vector integrals of motion: the angular momentum, and the Runge-Lenz vector, The orbital elements are related to them in the following way: sin ι cos ω ¼x ·Ẑ ¼ In the main text we need the Fourier decompositions of the Keplerian functions rðtÞ, θðtÞ, and several related quantities. These are conveniently derived by rewriting the relevant Fourier integrals as integrals over the eccentric anomaly E and using Eqs. (A1a), (A1c), and (A2). In this way one obtains [92] r a ¼ 1 þ e 2 2 − 2e X ∞ n¼1 J 0 n ðneÞ n cos nω bt ; ðA6aÞ r 2 a 2 ¼ 1 þ