Axion resonances in binary pulsar systems

TL;DR

This paper studies how ultra-light axion backgrounds can cause observable instantaneous and secular effects on binary pulsar systems, potentially allowing detection of axions through pulsar timing measurements.

Contribution

It introduces a detailed analysis of instantaneous axion-induced effects on binary pulsars using N-body simulations, extending beyond previous focus on secular resonance effects.

Findings

Instantaneous effects are significant over a wide range of configurations.

Signal-to-noise ratio increases with binary eccentricity.

Detection prospects depend on pulsar proximity and timing accuracy.

Abstract

We investigate the extent to which resonances between an oscillating background of ultra-light axion and a binary Keplerian system can affect the motion of the latter. These resonances lead to perturbations in the instantaneous time-of-arrivals, and to secular variations in the period of the binary. While the secular changes at exact resonance have recently been explored, the instantaneous effects have been overlooked. In this paper, we examine the latter using N-body simulations including the external oscillatory forcing induced by the axion background. While the secular effects are restricted to a narrow width near the resonance, the instantaneous changes, albeit strongest close to resonances, are apparent for wide range of configurations. We compute the signal-to-noise ratio (SNR) as a function of semi-major axis for a detection of axion oscillations through the R\{o} mer delay. The latter can be extracted from the time-of-arrivals of binary pulsars. The SNR broadly increases with increasing binary eccentricity in agreement with the secular expectation. However, we find that it differs significantly from the scaling a^{5/2} around the lowest orders of resonance. Future observations could probe these effects away from resonances and, therefore, constrain a much broader range of axion masses provided that binary pulsar systems are found near the central region of our Galaxy, and that the time-or-arrival measurement accuracy reaches < 10 ns

Equations

$$\text{mbox}resonancecondition:{\Omega }_{c}k=2m_a,k\in \mathbb{N}$$

$$\phi(\mathbf{x},t)=\psi(\mathbf{x})\cos\big{(}m_{a}t+\varphi({\bf x})\big{)}$$

$$P_{\mathrm{DM}}(\mathbf{x},t)\supset -\frac{1}{2}{m}_a^2{\psi }^2(\mathbf{x})\cos (2m_{a}t+2\phi ),$$

$$d{s}^2=(1+2\Psi )d{t}^2+(1-2\Phi ){\delta }_{ij}d{x}^{i}d{x}^j,$$

$$\Psi (\mathbf{x},t)\supset {\Psi }_{2c}(\mathbf{x})\cos ({\omega }_{a}t+2\phi )$$

$${\Psi }_{\text{2c}}(\mathbf{x})=\frac{\pi G}{m_a^2}{\rho }_{\text{DM}}(\mathbf{x}),{\rho }_{\text{DM}}(\mathbf{x})=\frac{1}{2}{m}_a^2{\psi }^2(\mathbf{x})$$

$$\mathbf{F}=-4\pi G{\rho }_{\mathrm{DM}}\mu \cos ({\omega }_{a}t+\phi )\mathbf{r}$$

$${\mathcal{H}}_{\text{pert}}=2\pi G{\rho }_{\mathrm{DM}}\mu \cos ({\omega }_{a}t+\phi )r^2.$$

$$ϵ\equiv \frac{{\mathcal{H}}_{\mathrm{pert}}^0}{{\mathcal{H}}_K^0}=\frac{2\pi G{\rho }_{\mathrm{DM}}\mu a_{\ast }^2}{GM\mu /a_{\ast }}=\frac{2\pi {\rho }_{\mathrm{DM}}{a}_{\ast }^3}{M},$$

$$a_{\ast }={\left(\frac{GM}{{\omega }_a^2}\right)}^{1/3}{k}^{2/3},k\in \mathbb{N}$$

$$\dot{T}=-6\frac{J_k(k{e}_{\ast })}{k}{T}^{2}G{\rho }_{\mathrm{DM}}\sin ({\gamma }_k+\phi ),$$

$$\mathcal{H}({\gamma }_k,{\hat{J}}_{\text{res}})$$

$$-\frac{3k^2{G}^2{M}^2{\mu }^3}{J_c^{\ast 4}}\frac{{\hat{J}}_{\text{res}}^2}{2}-\frac{2ϵ}{k^2}{J}_k(k{e}_{\ast })\cos ({\gamma }_k+\phi )\frac{GM\mu }{a_{\ast }}.$$

$$\frac{d{\gamma }_k}{d\tau }$$

$$\frac{d{\hat{J}}_{\text{res}}}{d\tau }$$

$$\frac{dt}{d\tau }=+\frac{\partial \mathcal{H}}{\partial (-E)}=1,$$

$$k\delta ({\hat{J}}_{\text{res}})\sim {\left(\frac{8ϵJ_k(k{e}_{\ast })}{3k^2}\right)}^{1/2}{J}_c^{\ast }.$$

$$\frac{\delta T}{T}=3\sqrt{\frac{8ϵJ_k(k{e}_{\ast })}{3k^2}},$$

$$\delta {\hat{J}}_{\text{res}}\cdot {\hat{J}}_{\text{res}}(t=0)\sim ϵ,$$

$$\text{mbox}i.e.\delta {\hat{J}}_{\text{res}}\sim {ϵ}^{1-p}.$$

$$\dot{T}\approx \frac{T_{\max }-T_0}{t_{\mathrm{res}}}=10^{-9}\frac{3T_0}{2t_{\mathrm{res}}}\approx 2\cdot 10^{-13}.$$

$$\left(\frac{S}{N}\right)^{2}=\frac{1}{\sigma_{\Delta}^{2}}\sum_{i=1}^{N}\big{[}\Delta t_{\text{\tiny TOA}}(t_{i})\big{]}^{2},$$

$$\dot{\mathbf{\theta }}=\frac{\partial \mathcal{H}}{\partial \mathbf{J}},\dot{\mathbf{J}}=-\frac{\partial \mathcal{H}}{\partial \mathbf{\theta }},$$

$$\frac{r^2}{a_{\ast }^2}-1=\frac{1}{2}{\alpha }_0(e)+\sum _{n=1}^{\infty }{\alpha }_n(e)\cos (n{\theta }_c).$$

$${\alpha }_n(e)=\frac{1}{\pi }{\int }_0^{2\pi }\left(\frac{r^2}{a_{\ast }^2}-1\right)\cos (n{\theta }_c)d{\theta }_c.$$

$${\alpha }_0(e)$$

$${\alpha }_n(e)$$

$$\frac{r^2}{a_{\ast }^2}=1+\frac{3}{2}{e}^2-\sum _{n=1}^{\infty }\frac{4}{n^2}{J}_n(ne)\cos (n{\theta }_c),$$

$$\cos ({\omega }_{a}t$$

$${\mathcal{H}}_0=-\frac{G^2{M}^2{\mu }^3}{2J_c^2}-E.$$

$${\mathcal{H}}_{\text{pert}}\supset {\mathcal{H}}_{\text{res}}=-\frac{2ϵ}{k^2}{J}_k(k{e}_{\ast })\cos ({\omega }_{a}t-k{\theta }_c+\phi )\frac{GM\mu }{a_{\ast }},$$

$$\mathcal{H}({\theta }_c,J_c,t,-E)$$

Full text

1 Introduction

The nature of dark matter is one of the greatest puzzles in physics. Since its existence was conjectured by Zwicky [1] from virial equilibrium considerations in galaxy clusters, a plethora of models have been proposed, ranging from weakly interacting massive particles (WIMPs) [2, 3] to modified gravity theories such as MOND [4], but many of them are now severely constrained by astrophysical/cosmological or particle physics experiments. Another possibility is that dark matter is a Bose-Einstein condensate of light bosons such as axions [see 5, 6, 7, for early work]. Ultra-light axions with a mass $m_{a}\gtrsim 10^{-22}\ {\rm eV}$ much lighter than that of QCD axions could help resolve some of the small-scale problems of standard cold dark matter (CDM) such as the core-cusp issue etc. [see, e.g. 8, 9, 10, for recent comprehensive discussions]. However, large Compton wavelengths leading to a significant suppression of the small-scale power spectrum are incompatible with Lyman-$\alpha$ forest measurements, which set a lower bound of $m_{a}\gtrsim 10^{-21}\rm\ eV$ (95% C.L.) [11, 12].

Pulsar timing offers another avenue to probe the existence of coherent oscillations induced by ultra-light scalar fields [13, 14, 15]. Millisecond radio pulsars (MSPs, for review see [16]) are known as the most stable clocks in the universe on timescales longer than a few spin periods of neutron star. It is, therefore, not surprising that double neutron star systems and neutron star - white dwarfs systems are actively used to study gravity (for review see [17]). In particular, the first indirect detection of gravitational waves emission [18, 19] was based on detail measurements of the orbital period decrease in a binary pulsar. Another example are tests of the equivalence principle using a triple containing a MSP and two white dwarfs [20]. The central parsec of the Milky Way could contain up to 10000 MSPs and a part of them will be discovered with the next generation radio telescopes such as SKA [21].

When a binary pulsar or MSP with a white dwarf companion orbits in a background of ultra-light axions, the dynamics of the binary is affected by the axion oscillations especially if the axion mass $m_{a}$ is close to integer multiples of half the orbital frequency $\Omega_{c}$. When there is such a resonance, axion oscillations induce a secular perturbation to the orbital period which might be observable with pulsar timing experiments [14, 22].

In this paper, we investigate the instantaneous perturbation rather than the secular effect using direct N-body simulations, in which the axion background is treated as a small perturbation to the binary’s Keplerian solution. Our analysis is, therefore, also valid away from resonances, which are restricted to a very narrow range of orbital frequencies (or semi-major axes).

The paper is organized as follows. In section §2, we provide a general overview of ultra-light axion dark matter and the resonance phenomenon. In 2.1, we discuss the perturbations of a Newtonian binary, while 2.2 and 2.3 focus on the secular change of the orbital period, and the instantaneous effect of resonances, respectively. In section §3, we introduce the N-body simulations and present the expected signals. In section §4 we discuss the prospects for a detection of axion oscillations through the measurement of the Rømer delay in pulsar timing experiments. We summarize our conclusions in section §5.

2 Theoretical background

Ultra-light axion fields with mass $m_{a}$ in the range $10^{-22}-10^{-20}\ \rm eV$ oscillate on characteristic timescales of a few hours to a few days, which are of the same order as the orbital periods $T$ of typical binary pulsar systems. More precisely, when the orbital frequency $\Omega_{c}=2\pi/T$ multiplied by an integer equals twice the axion mass $m_{a}$,

[TABLE]

there is a resonance between the motion of the binary and the axion background. Let us first revisit the analysis of [13, 14, 15, 22].

Note that, unless stated otherwise, we use natural units $c=\hbar=k_{B}=1$ throughout.

2.1 Oscillatory Perturbations of a Newtonian Binary

The evolution of the ultra-light axion field is described by the Klein-Gordon equation, which solution is of the form [e.g. 8, for a review]

[TABLE]

where $\psi$ is a slowly varying complex amplitude (we thus ignore its time dependence since it is much longer than the timescales involved in our problem), $m_{a}$ is the axion mass and $\varphi({\bf x})$ is a position-dependent phase.

The energy density and pressure of the axion field can be extracted from the stress-energy tensor. They can be separated into two components: a time-independent piece, and an oscillating term arising from the coherent oscillations of the axion field. While the dominant contribution to the energy density slowly varies with time and scales like $\propto\psi^{2}({\bf x})\propto a^{-3}$, the dominant contribution to the pressure oscillates on a timescale $\propto m_{a}^{-1}$ according to

[TABLE]

that is, a characteristic frequency $\omega_{a}=2m_{a}$ twice the axion mass.

In order to investigate the effects of this perturbation, we consider a perturbed metric in the Newtonian gauge

[TABLE]

in which $\Psi$ is the Newtonian gravitational potential sourced by the dark matter distribution. The spatial part of Einstein’s equations shows that the time-independent contributions to the scalar potentials satisfy $\Psi=\Phi$ [the time-independent anisotropic stress induced by axion field is negligible, see e.g. 23], while the time-dependent pressure generates a time-dependent gravitational potential [13]

[TABLE]

with an amplitude

[TABLE]

The effect is proportional to the axion density and, therefore, is largest at the center of the axion solitonic cores where the density reaches $\rho_{\rm DM}({\bf x})\approx 10^{2}(m_{a}/10^{-22}\ \rm eV)^{2}\ \rm M_{\odot}\ pc^{-3}$ [24] for present-day Milky-Way size halos.

2.2 Secular change of the orbital period

Consider the motion of a binary with separation $a$, total mass $M$ and reduced mass $\mu$ in a background of ultra-light axions. In the local Fermi frame attached to the binary center-of-mass, the perturbative force induced by the oscillatory background of axions is [14]

[TABLE]

where $\omega_{a}=2m_{a}$. This adds a time-dependent piece to the Hamiltonian given by

[TABLE]

The strength of the perturbation is quantified by the ratio of the typical strength of the perturbative and Keplerian Hamiltonians $\mathscr{H}_{\rm pert}^{0}$ and $\mathscr{H}_{K}^{0}$:

[TABLE]

where $a_{*}$ is the semi-major axis (SMA) at the $k$-th order of resonance. From Kepler’s law, it is given by

[TABLE]

Unsuprisingly, we have $\epsilon\propto a_{*}^{3}$ since the relative amplitude of the axion oscillations depends on dark matter mass enclosed by the orbit.

We show in Appendix §A that, upon isolating the resonant part of the perturbed Hamiltonian, the change in the orbital period derived by [14] through a secular averaging is recovered:

[TABLE]

Here, $J_{k}$ is the $k$-th Bessel function, whereas $e_{*}$ and $\gamma_{k}$ are the corresponding eccentricity and angle evaluated for the resonance condition Eq. (2.1).

The impact of a resonance reveals itself in the long-term evolution of the system, rather than the instantaneous accelerations. Its increased influence on the motion of the system arises due to the extended coherence of the forces. However, we emphasize that the system also accumulates instantaneous changes which, as we shall see below, gather to somewhat different effects as one moves away from resonances.

2.3 Motion close to resonance

To reveal the effect of resonances (see e.g. [25]), we adopt the following strategy: we work in the extended phase-space, so that the collection of points that satisfy the resonance condition ($\omega_{a}$ being an integer multiple of $\Omega_{c}$) is an hypersurface of the extended phase space – the resonant surface. Since we are chiefly interested in points near this surface, we transform out the angle-action variables so that the resonance condition no longer is a condition on a linear combination of the angles, diagonalising it in effect. This results in a new set of angle-action variables, which includes a slowly evolving resonant angle and its conjugate momentum.

More precisely, the extended phase-space is constructed by introducing the canonical (angle-action) variables $(t,-E)$, where the time coordinate $t$ is now an additional angle-like variable, and $-E$ is the corresponding generalized momentum. A new variable $\tau$ is introduced to parametrize the phase space trajectories 111In analogy with Special Relativity, $\tau$ plays the role of proper time while $t$ is the time coordinate. Note that, in the analytical mechanics literature, both symbols are usually interchanged.. Furthermore, the Hamiltonian transforms as $\mathscr{H}\to\mathscr{H}_{\textrm{new}}=\mathscr{H}-E$. We shall henceforth drop the subscript ‘new’ on the extended Hamiltonian to avoid notational encumbrance. The resonant coordinates (say, for the $k$-th resonance) are defined as follows: the resonant angle is $\gamma_{k}=\omega_{a}t-k\theta_{c}$, where $\theta_{c}$ is the true anomaly, whereas the corresponding resonant action is $J_{\rm res}$. We refer the reader to Appendix A for more details.

The resonant coordinates $\gamma_{k}$ and $J_{\text{res}}$ do not change when the system is on the resonant surface, while they change slowly near resonance. Hence, we move on to expand the Hamiltonian about the resonant surface in powers of $\epsilon$, as well as average over the non-resonant angles (which evolve fast). All the non-resonant actions are integrals of motion up to first order in $\epsilon$, and assume their values at the resonant surface. Furthermore, the non-resonant part of the Hamiltonian oscillates rapidly in the non-resonant – or fast – angles, and average out to zero. Therefore, we conveniently define $J_{\textrm{res}}^{*}$ as the value of the resonant action $J_{\textrm{res}}$ at the resonant surface. If its deviation from this surface $\hat{J}_{\textrm{res}}=J_{\textrm{res}}-J_{\textrm{res}}^{*}$ is small (i.e. of order $\mathcal{O}(\sqrt{\epsilon})$) as the system evolves around the resonant surface, we can expand the Hamiltonian in powers of $\hat{J}_{\text{res}}$. Retaining terms up to quadratic order, that is, of order $\epsilon$, we arrive at (cf. Appendix A)

[TABLE]

On the right-hand side, the first line is an irrelevant constant while the second line is the Hamiltonian of a simple pendulum, in which $\hat{J}_{\textrm{res}}$ plays the role of the momentum and $\gamma_{k}$ the angle subtended relative to the vertical axis. The equations of motion are

[TABLE]

These equations are supplemented by

[TABLE]

which shows that $\tau\equiv t$ without loss of generality. There are two fixed points $(\gamma_{k},\hat{J}_{\text{res}})=(-\varphi\pm\pi,0)$ and $(-\varphi,0)$, respectively stable and unstable. They correspond to the resonant (phase space) orbits, which arise when i) the resonance condition is satisfied and ii) the oscillatory force induced by the axion condensate is maximally outwards (stable fixed point) or inwards (unstable fixed point) at periapsis. Overall, $\mathscr{H}$ generates two distinct classes of orbits: libration (small oscillations about the stable equilibrium) very near the resonance; and circulation or rotation (the pendulum completes whole revolutions) as one moves away from it. A phase space portrait is shown in Fig.1.

During such periodic motions, which take place on a timescale $\propto 1/\sqrt{\epsilon}$, one can derive the typical variation in $k\hat{J}_{\textrm{res}}=J_{c}^{*}-J_{c}$ by considering e.g. the average kinetic energy of a pendulum initially at rest at the unstable equilibrium point (i.e. at the top). The result is

[TABLE]

The scaling $k\delta(\hat{J}_{\text{res}})\propto\sqrt{\epsilon}J_{c}^{*}$ is actually valid within the entire region of libration, and inside the region of rotation so long as corrections of order $\mathcal{O}(\hat{J}_{\text{res}}^{3})$ to (2.12) are negligible.

Recall that $J_{c}$ is associated with the semi-major axis of the binary, and thus determines its period according to Kepler’s third law. Using this relation between $J_{c}$ and $T$, we observe that the period changes by

[TABLE]

over one revolution of the pendulum. However, a single such revolution corresponds to many orbits of the binary and, therefore, represents a cumulative effect of order $\sqrt{\epsilon}$. It is noteworthy that our particular $\mathscr{H}_{\textrm{pert}}$, being of an oscillatory nature, does not induce a net drift in phase space – if the system starts out in the vicinity of the resonant surface, it will remain there but undergo oscillations about it with a magnitude $\sim\sqrt{\epsilon}$. The pendulum approximation developed in this section is valid at all times, albeit for systems near a resonance solely. For $k=1$, $e_{*}=1/2$ and $\epsilon\approx 10^{-13}$, we obtain $\delta T\approx 4\times 10^{-7}~{}T$. This change is accumulated over about $10^{6.5}$ periods. This is considerably greater than the $\mathcal{O}(\epsilon)$ change accumulated over a few orbital times (i.e. a time interval much smaller than the libration timescale), but is only visible when observing the binary for a fairly long amount of time (depending on the precision of the observations, of course).

In addition to the narrow libration region of width $\sim\sqrt{\epsilon}$, fast circulation far away from resonances affect the instantaneous position of the binary. This effect is not strictly taken into account in the resonant formalism developed in [25] (but cf. appendix B), but we can treat it approximately as follows. Suppose the system is given an initial condition $\mathcal{O}(\epsilon^{p})$ from the resonant surface, where $0<p<1/2$. This is still near the resonance. However, (2.12) will hold only if the ”error” term $\mathcal{O}(\hat{J}_{\textrm{res}}^{3})$ is negligible relative to the ”kinetic” term $\propto\hat{J}_{\text{res}}^{2}$ and the potential term $\propto\epsilon$. The other harmonics in $\mathscr{H}_{\textrm{pert}}$ involving both $\beta$ and $\gamma_{k}$, still average out to zero on a dynamical timescale and, thus, are irrelevant for the long-time evolution of the system. To assess whether $\mathcal{O}(\hat{J}_{\textrm{res}}^{3})$ can be neglected, notice that, if $\hat{J}_{\textrm{res}}=c_{1}\epsilon^{p}$ at $t=0$, then the ”kinetic” term dominates over both the ”potential” and ”error” term since it is of order $\mathcal{O}(\epsilon^{2p})\gg\mathcal{O}(\epsilon,\epsilon^{3p})$ at $t=0$. If we restrict ourselves to the range $1/3<p<1/2$, then the “potential” term is the second dominant term whereas the ”error” term is still negligible, so that (2.12) is recovered.

Let us pursue the discussion with this additional restriction for simplicity. In this regime, the variation of $\hat{J}_{\textrm{res}}^{2}$ is $\mathcal{O}(\epsilon)$, which in turn implies

[TABLE]

One is led to the conclusion that, if one plots the variation in $J_{c}$ – which is linear in $\hat{J}_{\textrm{res}}$ – as a function of the initial distance from $J_{c}^{*}$, this variation (owing to the proximity to the resonant surface) behaves like $\epsilon/\left|J_{c}(t=0)-J_{c}^{*}\right|$ (cf. figure 3).

Initial conditions further away than $\mathcal{O}(\epsilon^{1/3})$ from the resonant surface are not captured by this approximation. In reality, $\epsilon$ is extremely small, whence the observed binaries will be most likely off resonance. We, therefore, require the aid of numerical simulations to study the behaviour far from resonance. This is the subject of sec. 3 below.

3 N-body simulations

In this Section, we numerically examine the dynamics of binary pulsar systems under the influence of scalar field oscillations. We consider axion masses in the range $10^{-22}-10^{-20}\ \rm eV$ bracketing the limits set thus far by the Lyman-$\alpha$ forest measurements [12, 11]. Furthermore, we purposely pick up a large dark matter density of $\rho_{\rm DM}=5\cdot 10^{3}\ \rm M_{\odot}\ pc^{-3}$ to ensure that the simulation results are not significantly affected by numerical noise. For an axion mass $m_{a}=10^{-21}\ \rm eV$, such a density is admittedly achieved only within the hypothetical, $\sim 1\ \rm pc$ - radius solitonic core residing in the central region of the Milky-Way. Notwithstanding, we emphasize that our choice of $\rho_{\rm DM}$ mainly serves the purpose of illustrating the response of the binary system near and away from resonances.

For our fiducial model, we choose an axion mass $m_{a}=10^{-21}\ \rm eV$ and a binary mass $M=2M_{\odot}$, so that the resonant values of the semi-major axis satisfy $a_{*}=0.205k^{2/3}\rm AU$. Higher binary masses would physically make more sense since the typical mass of Neutron stars is expected to be larger than the Chandrasekhar mass of $1.4M_{\odot}$. Nevertheless, this approximation does not hinder our illustration purposes. Higher binary masses will simply shift the resonances as can be seen from Eq. (2.10).

To examine the dynamics of the binary pulsar perturbed by the oscillating axion field, we integrate the binary with the extra instantaneous force given in Section §2.2 using direct N-body simulations. For the N-body integration, we use the publicly available code REBOUND [26]. We use IAS15, a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits [27]. Our numerical simulations can capture the impact of the axion oscillations even near resonances, as the system is integrable to first order in the perturbation (see appendix B).

Fig. 2 shows the evolution of the SMA as a function of time for resonant and non-resonant orbits with various eccentricities. The resonant orbits with SMA $a_{*}=0.2052\ \rm AU$ accumulate instantaneous changes to achieve a fractional perturbation of order $|a(t)/a_{0}-1|\sim 10^{-9}$ after a few hundred periods. They exhibit a periodic pattern of characteristic timescale of $\sim 1000\ \rm yr$ consistent with the $\sin(\gamma_{k}+\varphi)$ dependence of the period derivative $\dot{T}$. By contrast, the non-resonant simulations (dashed) have no long time coherence (i.e. $a/a_{0}-1$ fluctuates on a very short timescale relative to the resonant orbits) and only lead to a fractional perturbation $|a(t)/a_{0}-1|$ of order $\sim 10^{-13}$, four orders of magnitude lower than in the resonant case, and still two orders of magnitude lower when compared to a resonant and nearly circular orbit ($e=0.01$).

Note that, in both cases, the circular orbits remain almost unperturbed, while larger eccentricities result in larger effects, as expected from [14].

In order to explore the structure of the resonances as a function of mass and SMA, we run a grid of initial conditions and plot the maximal change $\Delta a/a_{0}$, where $\Delta a\equiv(\max(a)-\min(a))/2$. We sample three axion masses $m_{a}=10^{-20},10^{-21},10^{-22}\ \rm eV$. For each mass, we initialize the orbit with $2000$ different values of the SMA, with a log-uniform sampling centered on the lowest resonant peaks. This leads to a total of $6000$ simulations. The initial eccentricity is $e=0.5$ and all other angles are set to zero. The total time of integration is $1000$ orbital periods, and is determined separately for each binary.

The results are presented in Fig. 3. The maximal relative variation of the SMA, $\Delta a/a_{0}$, exhibits a series of resonant peaks, whose width and amplitude hardly changes as the axion mass is varied, except for an overall frequency shift. This shows that the resonant pattern is driven by the ratio $\Omega/\omega_{a}$. The structure of those peaks is given by equation (2.16), when transforming $T$ to $a$ via Kepler’s third law. Note also that $\Delta a/a_{0}$ scales like $\propto a^{3/2}$ when $a\gg a_{0}$. This reflects the scaling $\dot{T}\propto T^{2}$ in Eq.(2.11), which implies $\dot{a}\propto a^{5/2}$ through Kepler’s 3rd law.

In order to compare our findings to the secular change derived in [14], we use the results obtained in Fig. 2 in the resonance regime. Namely, we note that the maximal change of $a$ is about $a_{\rm max}\sim a_{0}(1+10^{-9})$. Converting this change into a perturbation of the orbital period, we find $T_{\rm max}\approx T_{0}(1+1.5\times 10^{-9})$. As a result, the derivative of the (secular averaged) period at the fundamental resonance ($k=1$) can be estimated as

[TABLE]

where, in our case, the unperturbed period is $T_{0}\approx 24\ \rm days$ and the time to reach the maximum resonance is $t_{\rm res}\approx 500\ \rm yr$. Ignoring the slowly varying term $\sin(\gamma_{k}+\varphi)$ in Eq. (2.11), we estimate the period derivative derived in [14] as $\dot{T}=2.83\cdot 10^{-13}$. This is likely an overestimate since the sine modulation $\sin(\gamma_{k}+\varphi)$ must gradually decrease during the evolution until $t=t_{\rm res}$, so that it vanishes at $T_{\rm max}$. The period derivative reflects this behaviour and also decreases monotonically to zero over the same time interval.

To conclude, our numerical simulations reproduce the secular results of [14] at exact resonance with reasonable accuracy. Most importantly, they allow us to explore the instantaneous changes induced by the axion oscillations near and far away from resonances, which was not considered in [14]. Our next goal is to calculate the signal-to-noise ratio (SNR) for a detection of these instantaneous changes in the measurements of the time of arrival (TOA) signal of an hypothetical binary pulsar. We will focus on an axion mass $m_{a}=10^{-21}\ \rm eV$ for illustration.

4 Prospects for detection with pulsar timing experiments

A detailed introduction to the pulsar timing technique can be found in, e.g., [28, 29]. The timing method is based on repeated measurements of the TOA of high signal-to-noise average pulses (whose shape is fairly stable). This technique has been used to monitor the TOAs of particularly interesting pulsars over the past decades. The data is further processed using pulsar timing packages, the most popular one being TEMPO2 [30].

The best timing accuracy achieved thus far by the Parkes radio telescope is approximately tens of $\mu$s [31] for 100-1000 of individual observations. The Square Kilometer Array (SKA [32]) might achieve a timing accuracy as high as $5$ ns, while routinely decrease the uncertainty in the TOA down to values of order 10 – 100 ns for a large number of MSPs [33].

Various non-relativistic and relativistic effects contribute to fluctuations in the TOA. The Newtonian contribution includes the Rømer time delay, while the relativistic corrections include the Einstein and Shapiro time delay [29]. Since we have explored the impact of axion oscillations in the Newtonian regime, we will focus on the Rømer time delay. A similar analysis can be carried out for relativistic pulsars using Post-Newtonian corrections to the Newtonian dynamics. This analysis is out of the scope of the current paper and should be performed elsewhere.

4.1 Signal-to-noise for the Rømer delay

We analyze the typical pattern of the TOA of the average pulses of a typical binary in and off resonance and discuss the SNR and the possibility to constrain the mass and density of axion dark matter cores.

The square of the SNR for a detection of axion coherent oscillations in the TOA measurements is

[TABLE]

$\Delta t_{\text{\tiny TOA}}(t_{i})=|\delta\textbf{\em r}(t_{i})|/c$ (note that we reintroduced the speed of light $c$ here, for the convenience of the calculation) is the difference in TOA induced by the axion oscillations, $|\delta\boldsymbol{r}(t_{i})|$ is the difference at time $t_{i}$ between the perturbed and unperturbed orbits (staring with the same initial conditions at time $t=0$), $t_{j}=j\Delta$ are the times at which the observations are performed, $N=t_{\rm obs}/\Delta$ is the total number of measurements and $\sigma_{\Delta}$ is the error on the measurement. The latter can be as small as $\sigma_{\Delta}=10^{-6}$ s when the pulse shape is averaged over a time interval $\Delta=10$ s. In practice however, TOA measurements are performed only during a fraction of the total observational time. Therefore, we shall make the conservative choice $\sigma_{\Delta}=10^{-6}$ s and $\Delta=10^{4}$ s which corresponds to a total precision of $\sigma_{\Delta}/\sqrt{N}\approx 5$ ns in our calculation of the SNR. On the one hand, this precision is comparable to that expected to be reached by the SKA on a short timescale (forty hours) for average brightness MSPs. On the other hand, it is rather conservative, given the total observation time of 10 years. We choose this accuracy because MSPs discovered in the Galactic centre region are expected to be dimmer than the known MSP population, so it should take more time to achieve a similar timing precision.

Fig. 4 shows the obtained signal to noise ratio pattern of binary pulsar as a function of their initial semi-major axis end eccentricity. For each eccentricity, we sample 1200 different values of the semi-major axis from a log-uniform grid in the range $a_{0}/\rm AU\in[0.15,0.9]$. We integrate each orbit twice, with and without the perturbations and record the difference in the positions $\delta\textbf{\em r}(t_{i})$ between the orbits, and plot the SNR defined in equation (4.1). Furthermore, we assume a DM density of $\rho_{\rm DM}=5\cdot 10^{3}\ \rm M_{\odot}\ pc^{-3}$, which corresponds to the central density of an axion core in a Milky-Way size halo [24]. This large density ensures that our measurements are well above the numerical noise. The overall amplitude of the SNR scales with $\epsilon\propto\rho_{\text{DM}}$.

The SNR displayed in Fig. 4 show, as expected, that the more eccentric binary pulsars ($e\gtrsim 0.75$) can be detected more easily, since higher orders of resonance are present and the overall SNR level away from the resonant peaks is fairly high. Intuitively, we might expect this behaviour since the instantaneous orbital velocity of more eccentric binary ’samples’ a wider range and could be sensitive to several resonances even away from resonance. By contrast, resonances in systems with lower eccentricities $e\lesssim 0.3$ become harder to detect. Moreover, as discussed in [14] and as we can conclude from considerations of conservation of angular momentum, (nearly) circular orbits are not significantly affected by resonances. It is important to stress that binary pulsars would most likely be found away from the resonant peaks. Overall, the ratio between the SNR for a system at resonance, and the SNR for a system at the mid-point between two adjacent resonances, can vary considerably. It is around two orders of magnitude for the lowest resonance, and generally decreases for higher order resonances.

Finally, let us also mention that, for SMA $a\gg a_{0}$ corresponding to much higher orders of resonance (not shown in Fig.4), we checked that the SNR behaves like $\propto a^{5/2}$, in agreement with the secular estimate of [14]. However, the average SNR departs significantly from this scaling around the lowest order resonances shown in Fig. 4.

4.2 Signal-to-noise as a function of axion mass

In Fig. 5, the SNR is shown as a function of the axion mass for a given dark matter density $\rho_{\rm DM}=5\cdot 10^{3}\ \rm M_{\odot}\ pc^{-3}$, and a few binary pulsar systems with identical (unperturbed) SMA but different eccentricities. As a result, $\epsilon$ is constant in these simulations. Therefore, the plateau at low values of $m_{a}$ indicates that the SNR is proportional to $\epsilon$ times a function which asymptotes to a constant for $m_{a}\ll\Omega$. The situation is less clear on the high frequency side of the fundamental reasonance, also because we have not pushed our simulations to axion mass values $m_{a}\gg\Omega$.

For the largest eccentricity, our large dark matter density yields a SNR of order unity away from the resonances. In the central region of our galaxy, $\rho_{\rm DM}\sim 1\ \rm M_{\odot}\ pc^{-3}$ (assuming there is no axion core), which translates into a SNR of order $10^{-4}$. If binary pulsar systems are found within $\sim 0.1{\rm Kpc}$ of the galactic center, then the required sensitivity of SNR$\sim 1$ could be achieved with a twofold improvement on the measurement error $\sigma_{\Delta}$ and a very dense sampling of the TAO (i.e. a very short time interval $\Delta$ between measurements). Note that our analysis does not take into account correlations in the TOA of different binary pulsars [see, e.g. 34, for related work]

5 Summary and conclusions

In this paper, we analyzed the motion of binary pulsars perturbed by an ultra-light axion oscillating background. We used numerical integrations to track the dynamics of the binary and its features. While we recovered the secular result of [14] near resonance, we emphasized that the near resonance orbits undergo a libration around a (stable) fixed point which is, strictly speaking, one of the only two possible orbits at exact resonance. These arise when the resonance condition is satisfied, and when the oscillatory force induced by the axion condensate is maximal (in absolute value) at the periapsis. In addition, we explored in details the short-term, instantaneous features near and far away from resonance which were not considered in previous literature. Although we adopt the ultra-light axions model throughout, our results straightforwardly generalize to other (possibly interacting) models of Bose-Einstein condensate (BEC) dark matter.

While the short-term effects and subsequent time-of-arrival measurements are stronger near the narrow resonances, especially for eccentric binaries, they are also apparent away from the resonant peaks. To quantify their imprint in real data, we considered time-of-arrival (TOA) measurements and computed the signal-to-noise (SNR) ratio required to observe this effects by current pulsar timing techniques. A detection - or lack - of anomalies in the TOA of binary pulsars would constrain the axion mass. The SNR can vary considerably depending on whether one sits near a resonance, or away from them. For the largest value of the eccentricity considered here ($e=0.75$), we find that the broadband SNR level remains approximately constant (resonances excluded) across two orders of magnitude in axion mass. Though current observations could probe these effects away from resonances only at densities of $\rho_{\rm DM}\approx 10^{4}\ \rm M_{\odot}\ pc^{-3}$, future instruments could reduce the sensitivity by at least one or two orders of magnitude. This does not take into account the possibility of cross-correlating TOA measurements from different binary systems.

Acknowledgments

We would like to thank Hagai B. Perets for helpful discussions. E.G. and Y.B.G. acknowledge support from the Technion Jacobs scholarship; V.D. acknowledges support by the Israel Science Foundation (grant no. 1395/16).

Appendix A The resonant bit of the perturbation

We work in the extended phase space $(\boldsymbol{\theta},\mathbf{J})$, where $\boldsymbol{\theta}$ denotes the four angles $(\theta_{a},\theta_{b},\theta_{c},t)$ and $\mathbf{J}$ the four actions $(J_{a},J_{b},J_{c},-E)$. The angle-action variables $(\theta_{a},\theta_{b},\theta_{c},J_{a},J_{b},J_{c})$ are the Delaunay elements of the orbit. In our simplified planar setup, $(\theta_{a},J_{a})$ are trivially conserved. Hamilton equations read

[TABLE]

where a dot designate a derivative w.r.t. the affine variable $\tau$ introduced to parametrize the phase space trajectory of the system. Note that $\dot{t}=1$ since the extended Hamiltonian $\mathscr{H}$ is constructed such that it depends on $E$ only through an additive term $-E$. As a result, $\dot{\mathscr{H}}=0$ along any orbit of the extended phase space.

The $r^{2}$-dependence of the perturbation Hamiltonian can be expanded into a Fourier series. As a result, $\mathscr{H}_{\text{pert}}$ depends on the actions $J_{b}$ and $J_{c}$, but only on the $2\pi$-periodic angle $\theta_{c}$ ($\theta_{b}$ is cyclic, so that $J_{b}$ is a constant of motion). The smallness of $\epsilon$ and the time-periodic nature of $\mathscr{H}_{\text{pert}}$ guarantee that the actual motion remains close to the unperturbed Keplerian orbit.

More precisely, since $\mathscr{H}_{\text{pert}}$ is an even function of $r^{2}$, the coefficients of the sine terms all come out zero upon integrating over one (unperturbed) orbit, and we are left with the cosine series

[TABLE]

Here, $a_{*}$ is the semi-major axis at the $k$-th resonance, see Eq.(2.10). An inner product with $\cos(k\theta_{c})$ gives

[TABLE]

Substitution of the eccentric anomaly $\eta$ turns this expression into

[TABLE]

Therefore

[TABLE]

where $J_{n}(x)$ are Bessel functions. In order to identify the resonant piece of the Hamiltonian, we use the trigonometric identity

[TABLE]

Here, $\varphi$ is the phase of the axion field relative to the unperturbed orbit (whose periapsis is assumed to occur at $t=0$). We assume that both $\omega_{a}$ and $\Omega_{c}=\dot{\theta_{c}}$ are positive (the same works if one is negative, or both) in the unperturbed system described by the Hamiltonian

[TABLE]

The resonance condition (2.1) thus reads $\omega_{a}=k\Omega_{c}({\bf J}^{*})$. An asterisk will hereafter designate quantities evaluated at the resonant surface. This defines the resonant surface $J_{c}=J_{c}^{*}\equiv(kG^{2}M^{2}\mu/\omega_{a})^{1/3}$, which is a line in the plane $(J_{c},E)$. Picking up the zero energy level (which is arbitrary) further selects a point $(J_{c}^{*},E^{*})$ on this surface.

To understand how $J_{c}$ and $E$ evolves near resonance, it is enough to consider the resonant piece of the perturbation Hamiltonian,

[TABLE]

where $e\approx e_{*}$ is the eccentricity at resonance. The difference $\mathscr{H}_{\text{pert}}-\mathscr{H}_{\text{res}}$ oscillates on a short (dynamical) timescale and, therefore, can be neglected. Therefore, we write the perturbed Hamiltonian $\mathscr{H}=\mathscr{H}_{0}+\mathscr{H}_{\text{pert}}$ (in the extended phase space) near the $k$-th resonance as

[TABLE]

The resulting equations of motion for the actions are

[TABLE]

The physical interpretation of these equations is straightfoward: $\dot{E}$ represents the rate at which energy is injected into the system (averaged over one cycle as we ignore the fast oscillations), whereas $\dot{J}_{c}$ describes the variation of the usual energy constructed from the conservative part of the potential (i.e. the 2-body gravitational interaction). Conservation of energy ensures that they are proportional to each other, i.e. $\frac{\dot{E}}{\omega_{a}}=\frac{\dot{J}_{c}}{k}$. The relation

[TABLE]

along with the second Hamilton’s equation $\dot{J}_{c}=-\partial\mathscr{H}_{\textrm{res}}/\partial\theta_{c}$, leads to the equation (2.11) derived in [14].

Consider now an orbit initially on the resonant surface, i.e. $(J_{c}(0),E(0))=(J_{c}^{*},E^{*})$. The resonant condition implies that, at leading order in $\epsilon$, $|E|$ and $|J_{c}|$ grow linearly with time. As a result, the orbit will gradually leave the resonance surface. However, the argument of the $\sin$ factor also starts evolving, and this will bring the orbit back towards the resonant surface. Since the derivatives $\dot{E}$ and $\dot{J}_{c}$ are in phase, the actual motion is an oscillation around $(J_{c}^{*},E^{*})$ on a line at an angle $\arctan(k/\omega_{a})$ with the resonant surface, i.e. $\frac{E-E^{*}}{\omega_{a}}=\frac{J_{c}-J_{c}^{*}}{k}$.

In order to reveal the oscillations around the resonant surface, it is convenient to work with (non-)resonant canonical variables. We define the (non-)resonant angle $\gamma_{k}$ (resp. $\beta$) through the canonical transformation

[TABLE]

The corresponding actions are

[TABLE]

The new set of variables $(\gamma_{k},\beta,J_{\text{res}},J_{\text{nr}})$ captures the near resonance orbits defined by the Hamiltonian (A.8), which have $J_{\text{nr}}=J_{\text{nr}}^{*}=0$ and $J_{\text{res}}$ oscillating around $J_{\text{res}}^{*}$. Only the orbits with $\gamma_{k}+\varphi=0$ or $\pm\pi$ never leave the resonant surface. They have $J_{\text{res}}=J_{\text{res}}^{*}$, and correspond to a stable ($\gamma_{k}+\varphi=\pm\pi$) and an unstable ($\gamma_{k}+\varphi=0$) fixed point in the Poincaré section $(\beta,J_{\text{nr}})=(0,0)$ of the phase space.

To proceed further, we recast the resonant part of the perturbation Hamiltonian into the form

[TABLE]

and use (A.12) to express the actions $J_{c}$ and $-E$ around the resonant surface as

[TABLE]

The final approximations arise, again, from the fact that the non-resonant actions oscillate rapidly around their resonant value. On defining $\hat{J}_{\text{res}}=J_{\text{res}}-J_{\text{res}}^{*}$, the previous relations allow us to expand the unperturbed part (A.6) of the Hamiltonian around $J_{c}^{*}$ and $-E^{*}$. Since

[TABLE]

the resonant condition implies that the term linear in $\hat{J}_{\text{res}}$ vanishes. This leads to equation (2.12), which shows that the near resonance orbits are described by a simple pendulum Hamiltonian. The only trajectories that are at exact resonance - the resonant orbits - correspond to the fixed points mentioned above.

To conclude this Section, we note that the fast oscillations can be treated - in a first approximation - as stochastic fluctuations around the smooth phase space trajectories described by (2.12).

Appendix B An Integral of Motion

Three integrals of motion are trivial to find – those are the two independent components of the angular momentum, and the extended phase space Hamiltonian. To further investigate the motion of the binary under the perturbation in equation (2.7), let us try to find an integral of motion that takes into account all the resonances.

The method is described in [35]. We consider a function $I(\boldsymbol{\theta},\mathbf{J})=I_{0}(J_{c})+\epsilon I_{1}(\boldsymbol{\theta},\mathbf{J})$. For it to be an invariant, it needs to have vanishing Poisson bracket with the Hamiltonian, that is, $\{I,\mathscr{H}\}=0$. This implies that, in the cosine expansion of appendix A,

[TABLE]

where $A_{n}$ comes from the Fourier expansion of $\mathscr{H}_{\textrm{pert}}=\sum_{n=-\infty}^{\infty}\epsilon A_{n}\cos(n\theta_{c}+\omega_{a}\tau)$, related to $a_{n}(e)$ via equations (2.9) and (A.3). We set $\varphi=0$ without restricting the validity of our argument.

The important part is that the left-hand-side of equation (B.1) vanishes at the resonances, so we take $\frac{dI_{0}}{dJ_{c}}$ to be a function that vanishes there, too. We choose

[TABLE]

whence (up to an arbitrary constant)

[TABLE]

where $\textrm{Si}(a,z)=\int_{0}^{z}t^{a-1}\sin tdt$. $I$ is a constant of motion to first order in $\epsilon$; together with the angular momentum components $J_{a},J_{b}$ and the extended phase-space Hamiltonian, we have $4$ integrals of motion, and therefore the system is integrable to first order in $\epsilon$; the system is restricted to move along phase space trajectories with $I=\textrm{const},~{}\mathscr{H}=\textrm{const}$. The total variation in $J_{c}$ may therefore be estimated by taking the variation of equation (B.3), and using the implicit function theorem.