Neutron stars and Pulsar timing arrays as Axion giant gyroscopes

Yiming Liu 7520220161@bit.edu.cn School of Physics, Beijing Institute of Technology, Beijing, 100081, China    Jinneng Luo 3120221503@bit.edu.cn School of Physics, Beijing Institute of Technology, Beijing, 100081, China    Sichun Sun sichunssun@gmail.com School of Physics, Beijing Institute of Technology, Beijing, 100081, China
Abstract

We consider the three-dimensional rotating motions of neutron stars blown by the “axion wind”. Neutron star precession and spin can change from the magnetic moment coupling to the oscillating axion background field, in analogy to the gyroscope motions with a driving force and the laboratory Nuclear Magnetic Resonance(NMR) detections of the axion. This effect modulates the pulse arrival time of the pulsar timing arrays. It shows up as a signal on the timing residual and two-point correlation function on the recent data of Nanograv and PPTA. The current measurement of PTAs can thus cast constraints on the axion-nucleon coupling as gann1012GeV1similar-tosubscript𝑔annsuperscript1012superscriptGeV1g_{\text{ann}}\sim 10^{-12}\text{GeV}^{-1}italic_g start_POSTSUBSCRIPT ann end_POSTSUBSCRIPT ∼ 10 start_POSTSUPERSCRIPT - 12 end_POSTSUPERSCRIPT GeV start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT

preprint: APS/123-QED

I Introduction

Oscillating background fields across the cosmic space, which behave like matter, are interesting dark matter candidates Hu et al. (2000); Hui et al. (2017); Arias et al. (2012); Feng (2010); Gibney (2020). The most common ones are either axionlike pseudoscalars or scalar fields such as moduli Preskill et al. (1983); Abbott and Sikivie (1983); Dine and Fischler (1983); Svrcek and Witten (2006). Since axion has a large parameter space, various axion detection schemes are proposed across the spectrum111https://github.com/cajohare/AxionLimits, from the quantum laboratory and collider studies to astrophysical observations. Inspired by the nucleon spin precession detection of axion particles such as CASPEr Garcon et al. (2017a), here we extend the detection principle using the microscopic nucleon spins to the macroscopic neutron stars, such that the three-dimensional rotating motion of spinning neutron stars blown by the axion wind can yield signals for a large range of parameter space for the axion mass and decay constant. Here in this paper we mainly discuss the case that ultralight axion wind blows on the pulsars in PTA, which yields the common spectrum on the residual with fc=ω2π2.42 nHz(m1023eV)subscript𝑓𝑐𝜔2𝜋similar-to-or-equals2.42 nHz𝑚superscript1023eVf_{c}=\frac{\omega}{2\pi}\simeq 2.42\text{\,nHz}\left(\frac{m}{10^{-23}\text{% eV}}\right)italic_f start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = divide start_ARG italic_ω end_ARG start_ARG 2 italic_π end_ARG ≃ 2.42 nHz ( divide start_ARG italic_m end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 23 end_POSTSUPERSCRIPT eV end_ARG ) and two-point correlation function curve. We also briefly comment on the modulation of the pulsar period with different axion masses e.g.QCD axion 242 MHz(m106eV)242 MHz𝑚superscript106eV242\text{\,MHz}\left(\frac{m}{10^{-6}\text{eV}}\right)242 MHz ( divide start_ARG italic_m end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT eV end_ARG ), or the resonant frequency with pulsar period 2.42 kHz(m1011eV)2.42 kHz𝑚superscript1011eV2.42\text{\,kHz}\left(\frac{m}{10^{-11}\text{eV}}\right)2.42 kHz ( divide start_ARG italic_m end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 11 end_POSTSUPERSCRIPT eV end_ARG ).

Recently several Pulsar Timing Arrays (PTAs) observatories including NANOGrav(North American Nanohertz Observatory for Gravitational Waves) Agazie et al. (2023a), Parkes PTA (PPTA) Reardon et al. (2023), European PTA (EPTA) Antoniadis et al. (2023a, b) and Chinese PTA (CPTA) Xu et al. (2023), report strong evidence for gravitational wave (GW) background at the nano-Hertz waveband, which is consistent with previous claims Arzoumanian et al. (2020); Chen et al. (2021); Goncharov et al. (2021); Antoniadis et al. (2022). Extensive studies on potential sources of the observed GW background Cang et al. (2023); Franciolini et al. (2023a, b); Liu et al. (2023); Ellis et al. (2023); Wu et al. (2023a); Li et al. (2023a); Sun et al. (2022); Li et al. (2023b); Battista and De Falco (2021); De Falco and Battista (2023); Konoplya and Zhidenko (2023); Ben-Dayan et al. (2023); Balaji et al. (2023); Kohri and Terada (2021); Inomata et al. (2023a); Vagnozzi (2021); Benetti et al. (2022); Vagnozzi (2023); Guo et al. (2023); Oikonomou (2023a, b); Choudhury (2023); Choudhury et al. (2023); Bhattacharya et al. (2023); Madge et al. (2023); Sun and Zhang (2021), including supermassive black holes Afzal et al. (2023); Middleton et al. (2021); Agazie et al. (2023b); Antoniadis et al. (2023c), merging PBHs Depta et al. (2023); Gouttenoire et al. (2023), phase transitions Bian et al. (2021); Arzoumanian et al. (2021); Xue et al. (2021); Wang (2022a) and axion topological defects Wang (2022b); Ferreira et al. (2023); Inomata et al. (2023b), etc are conducted. Nongravitational wave signals, such as the ultralight dark matter can also be observed through the timing residuals Khmelnitsky and Rubakov (2014); Porayko and Postnov (2014); Porayko et al. (2018); Sun et al. (2022); Wu et al. (2023b); Omiya et al. (2023); Nomura et al. (2020).

Most of the pulsars we observed are at the distances of order kpc from the Earth. The de Broglie wavelength of the ultralight dark matter is around λdB=2πmv4kpc(1023eVm)(103v)subscript𝜆dB2𝜋𝑚𝑣similar-to-or-equals4kpcsuperscript1023eV𝑚superscript103𝑣\lambda_{\text{dB}}=\frac{2\pi}{{m}v}\simeq 4\text{kpc}\left(\frac{10^{-23}% \text{eV}}{m}\right)\left(\frac{10^{-3}}{v}\right)italic_λ start_POSTSUBSCRIPT dB end_POSTSUBSCRIPT = divide start_ARG 2 italic_π end_ARG start_ARG italic_m italic_v end_ARG ≃ 4 kpc ( divide start_ARG 10 start_POSTSUPERSCRIPT - 23 end_POSTSUPERSCRIPT eV end_ARG start_ARG italic_m end_ARG ) ( divide start_ARG 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT end_ARG start_ARG italic_v end_ARG ). Fuzzy dark matter oscillates with the frequency 2.42 nHz(m1023eV)2.42 nHz𝑚superscript1023eV2.42\text{\,nHz}\left(\frac{m}{10^{-23}\text{eV}}\right)2.42 nHz ( divide start_ARG italic_m end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 23 end_POSTSUPERSCRIPT eV end_ARG ) and the oscillation period Tc=fc113.1yr(1023eVm)subscript𝑇𝑐superscriptsubscript𝑓𝑐1similar-to-or-equals13.1yrsuperscript1023eV𝑚T_{c}={f}_{c}^{-1}\simeq 13.1\text{yr}\left(\frac{10^{-23}\text{eV}}{m}\right)italic_T start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = italic_f start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ≃ 13.1 yr ( divide start_ARG 10 start_POSTSUPERSCRIPT - 23 end_POSTSUPERSCRIPT eV end_ARG start_ARG italic_m end_ARG ). The field behaves like pressureless cold dark matter on the cosmic scale. Notably, the oscillation frequencies of fuzzy dark matter models fell into the sensitive region of PTA observations.

II Axion background coupling to the pulsar magnetic moment

Taking into account the huge occupation number of the dark matter particle in the universe the axion field in the Galaxy can be viewed as a classical field and be approximated byDuan et al. (2023); Khmelnitsky and Rubakov (2014)

a(t,x)a0cos(ωt+mavx+ϕ0)𝑎𝑡𝑥subscript𝑎0𝜔𝑡subscript𝑚𝑎𝑣𝑥subscriptitalic-ϕ0\displaystyle a(t,\vec{x})\approx a_{0}\cos\left(-\omega t+m_{a}\vec{v}\cdot% \vec{x}+\phi_{0}\right)italic_a ( italic_t , over→ start_ARG italic_x end_ARG ) ≈ italic_a start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT roman_cos ( - italic_ω italic_t + italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT over→ start_ARG italic_v end_ARG ⋅ over→ start_ARG italic_x end_ARG + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) (1)
ω=ma(1+12v2)𝜔subscript𝑚𝑎112superscript𝑣2\displaystyle\omega=m_{a}\left(1+\frac{1}{2}v^{2}\right)italic_ω = italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ( 1 + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_v start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) (2)
a0=2ρDMmasubscript𝑎02subscript𝜌𝐷𝑀subscript𝑚𝑎\displaystyle a_{0}=\frac{\sqrt{2\rho_{DM}}}{m_{a}}italic_a start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = divide start_ARG square-root start_ARG 2 italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG end_ARG start_ARG italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT end_ARG (3)

where masubscript𝑚𝑎m_{a}italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT is the mass of the axion, v𝑣vitalic_v is the typical velocity of the axion and ρDM0.3GeV/cm3similar-tosubscript𝜌𝐷𝑀0.3GeVsuperscriptcm3\rho_{DM}\sim 0.3\mathrm{GeV/cm^{3}}italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT ∼ 0.3 roman_GeV / roman_cm start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT is the local dark matter density. The Hamiltonian of the interaction between the axion and the nucleon is given byGraham and Rajendran (2013); Berlin et al. (2024)

H(t,x)=gaNNa(t,x)σ𝐻𝑡𝑥subscript𝑔𝑎𝑁𝑁𝑎𝑡𝑥𝜎H(t,\vec{x})=g_{aNN}\nabla a(t,\vec{x})\cdot\vec{\sigma}italic_H ( italic_t , over→ start_ARG italic_x end_ARG ) = italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT ∇ italic_a ( italic_t , over→ start_ARG italic_x end_ARG ) ⋅ over→ start_ARG italic_σ end_ARG (4)

where gaNNsubscript𝑔𝑎𝑁𝑁g_{aNN}italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT is the coulping constant and σ𝜎\vec{\sigma}over→ start_ARG italic_σ end_ARG is the nuclear-spin operator.

As it is known that the magnetic moment μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG can be expressed in terms of the nuclear spin σ𝜎\vec{\sigma}over→ start_ARG italic_σ end_ARG

μ=γσ𝜇𝛾𝜎\vec{\mu}=\gamma\vec{\sigma}over→ start_ARG italic_μ end_ARG = italic_γ over→ start_ARG italic_σ end_ARG (5)

where γ𝛾\gammaitalic_γ is the gyromagnetic ratio of the nuclear spin, related to the property of the nucleon. Define the effective magnetic field induced by the axionGarcon et al. (2017b):

BALP(t)=gaNN2ρDMvγsin(ωt+ϕ0)subscript𝐵𝐴𝐿𝑃𝑡subscript𝑔𝑎𝑁𝑁2subscript𝜌𝐷𝑀𝑣𝛾𝜔𝑡subscriptitalic-ϕ0\vec{B}_{ALP}(t)=g_{aNN}\sqrt{2\rho_{DM}}\frac{\vec{v}}{\gamma}\sin\left(-% \omega t+\phi_{0}\right)over→ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_A italic_L italic_P end_POSTSUBSCRIPT ( italic_t ) = italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG 2 italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG divide start_ARG over→ start_ARG italic_v end_ARG end_ARG start_ARG italic_γ end_ARG roman_sin ( - italic_ω italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) (6)

therefore we can write Eq.(4) into the form below, similar to the Hamiltonian of the particle possessing the magnetic moment in an external magnetic field:

H(t)=μBALP(t)𝐻𝑡𝜇subscript𝐵𝐴𝐿𝑃𝑡H(t)=-\vec{\mu}\cdot\vec{B}_{ALP}(t)italic_H ( italic_t ) = - over→ start_ARG italic_μ end_ARG ⋅ over→ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_A italic_L italic_P end_POSTSUBSCRIPT ( italic_t ) (7)

Now we generalize the above microscopic interaction to the macroscopic astrophysical object with a huge magnetic moment and the effective axion dark matter magnetic field.

Classically, the pulsar’s magnetic moment μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG and spin angular momentum L𝐿\vec{L}over→ start_ARG italic_L end_ARG are given by:

μ=R32Bp𝜇superscript𝑅32subscript𝐵𝑝\displaystyle\vec{\mu}=\frac{R^{3}}{2}\vec{B}_{p}over→ start_ARG italic_μ end_ARG = divide start_ARG italic_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG start_ARG 2 end_ARG over→ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT (8)
L=IΩ0=25MR2Ω0𝐿𝐼subscriptΩ025𝑀superscript𝑅2subscriptΩ0\displaystyle\vec{L}=I\vec{\Omega}_{0}=\frac{2}{5}MR^{2}\vec{\Omega}_{0}over→ start_ARG italic_L end_ARG = italic_I over→ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG 5 end_ARG italic_M italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT over→ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT (9)

where Bpsubscript𝐵𝑝\vec{B}_{p}over→ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT is the magnetic field of the pulsar’s pole, M𝑀Mitalic_M represents the pulsar’s mass, R𝑅Ritalic_R stands for the radius of the pulsar and Ω0subscriptΩ0\vec{\Omega}_{0}over→ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is pulsar’s angular velocity.

Substitute Eq.(8), (9) into Eq.(5), and we can get the gyromagnetic ratio of the pulsar:

γ=R3Bp/22MR2Ω0/5=5RBp4MΩ0𝛾superscript𝑅3subscript𝐵𝑝22𝑀superscript𝑅2subscriptΩ055𝑅subscript𝐵𝑝4𝑀subscriptΩ0\gamma=\frac{{R^{3}B_{p}}/{2}}{{2MR^{2}\Omega_{0}}/{5}}=\frac{5RB_{p}}{4M% \Omega_{0}}italic_γ = divide start_ARG italic_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_B start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT / 2 end_ARG start_ARG 2 italic_M italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT / 5 end_ARG = divide start_ARG 5 italic_R italic_B start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT end_ARG start_ARG 4 italic_M roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG (10)

Then we can arrive at the energy change of the pulsar induced by the motion of the axion

H(t)=vgaNN2ρDM25MR2Ω0B^p(t)v^sin(ωt+ϕ0)𝐻𝑡𝑣subscript𝑔𝑎𝑁𝑁2subscript𝜌𝐷𝑀25𝑀superscript𝑅2subscriptΩ0subscript^𝐵𝑝𝑡^𝑣𝜔𝑡subscriptitalic-ϕ0H\left(t\right)=-vg_{aNN}\sqrt{2\rho_{DM}}\frac{2}{5}MR^{2}\Omega_{0}\hat{B}_{% p}\left(t\right)\cdot\hat{v}\sin\left(-\omega t+\phi_{0}\right)italic_H ( italic_t ) = - italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG 2 italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG divide start_ARG 2 end_ARG start_ARG 5 end_ARG italic_M italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) ⋅ over^ start_ARG italic_v end_ARG roman_sin ( - italic_ω italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) (11)

where we have defined three unit vectors in the respective orientation:B^p(t)=Bp(t)/Bp,Ω^0=Ω0/Ω0,v^=v/v:absentformulae-sequencesubscript^𝐵𝑝𝑡subscript𝐵𝑝𝑡subscript𝐵𝑝formulae-sequencesubscript^Ω0subscriptΩ0subscriptΩ0^𝑣𝑣𝑣\colon\hat{B}_{p}\left(t\right)={\vec{B}_{p}\left(t\right)}/{B_{p}},\hat{% \Omega}_{0}={\vec{\Omega}_{0}}/{\Omega_{0}},\hat{v}={\vec{v}}/{v}: over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) = over→ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) / italic_B start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT , over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = over→ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT / roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , over^ start_ARG italic_v end_ARG = over→ start_ARG italic_v end_ARG / italic_v. Because the magnetic axis B^p(t)subscript^𝐵𝑝𝑡\hat{B}_{p}\left(t\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) rotates around the spin axis Ω^0subscript^Ω0\hat{\Omega}_{0}over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT with angular velocity Ω0subscriptΩ0\Omega_{0}roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, we can decompose B^p(t)subscript^𝐵𝑝𝑡\hat{B}_{p}\left(t\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) as::\colon:

B^p(t)subscript^𝐵𝑝𝑡\displaystyle\hat{B}_{p}\left(t\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) =[B^p(0)(B^p(0)Ω^0)Ω^0]cosΩ0t+absentlimit-fromdelimited-[]subscript^𝐵𝑝0subscript^𝐵𝑝0subscript^Ω0subscript^Ω0subscriptΩ0𝑡\displaystyle=\left[\hat{B}_{p}\left(0\right)-\left(\hat{B}_{p}\left(0\right)% \cdot\hat{\Omega}_{0}\right)\hat{\Omega}_{0}\right]\cos{\Omega_{0}t}+= [ over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) - ( over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ⋅ over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ] roman_cos roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t + (12)
[Ω^0×B^p(0)]sinΩ0t+(B^p(0)Ω^0)Ω^0delimited-[]subscript^Ω0subscript^𝐵𝑝0subscriptΩ0𝑡subscript^𝐵𝑝0subscript^Ω0subscript^Ω0\displaystyle\left[\hat{\Omega}_{0}\times\hat{B}_{p}\left(0\right)\right]\sin{% \Omega_{0}t}+\left(\hat{B}_{p}\left(0\right)\cdot\hat{\Omega}_{0}\right)\hat{% \Omega}_{0}[ over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT × over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ] roman_sin roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t + ( over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ⋅ over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT

Define:

P=[B^p(0)v^(B^p(0)Ω^0)(Ω^0v^)]𝑃delimited-[]subscript^𝐵𝑝0^𝑣subscript^𝐵𝑝0subscript^Ω0subscript^Ω0^𝑣\displaystyle P=\left[\hat{B}_{p}\left(0\right)\cdot\hat{v}-\left(\hat{B}_{p}% \left(0\right)\cdot\hat{\Omega}_{0}\right)\left(\hat{\Omega}_{0}\cdot\hat{v}% \right)\right]italic_P = [ over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ⋅ over^ start_ARG italic_v end_ARG - ( over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ⋅ over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ( over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⋅ over^ start_ARG italic_v end_ARG ) ] (13)
Q=[(Ω^0×B^p(0))v^]𝑄delimited-[]subscript^Ω0subscript^𝐵𝑝0^𝑣\displaystyle Q=\left[\left(\hat{\Omega}_{0}\times\hat{B}_{p}\left(0\right)% \right)\cdot\hat{v}\right]italic_Q = [ ( over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT × over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ) ⋅ over^ start_ARG italic_v end_ARG ] (14)
S=(B^p(0)Ω^0)(Ω^0v^)𝑆subscript^𝐵𝑝0subscript^Ω0subscript^Ω0^𝑣\displaystyle S=\left(\hat{B}_{p}\left(0\right)\cdot\hat{\Omega}_{0}\right)% \left(\hat{\Omega}_{0}\cdot\hat{v}\right)italic_S = ( over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) ⋅ over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ( over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⋅ over^ start_ARG italic_v end_ARG ) (15)

We then arrive at:

B^p(t)v^=PcosΩ0t+QsinΩ0t+Ssubscript^𝐵𝑝𝑡^𝑣𝑃subscriptΩ0𝑡𝑄subscriptΩ0𝑡𝑆\hat{B}_{p}\left(t\right)\cdot\hat{v}=P\cos{\Omega_{0}t}+Q\sin{\Omega_{0}t}+Sover^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( italic_t ) ⋅ over^ start_ARG italic_v end_ARG = italic_P roman_cos roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t + italic_Q roman_sin roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t + italic_S (16)

Assume that all the energy is converted into the rotational energy of the pulsar::\colon:

12IΩ2(t)12IΩ02=H(t)12𝐼superscriptΩ2𝑡12𝐼superscriptsubscriptΩ02𝐻𝑡\frac{1}{2}I\Omega^{2}\left(t\right)-\frac{1}{2}I\Omega_{0}^{2}=H\left(t\right)divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_I roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_t ) - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_I roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = italic_H ( italic_t ) (17)

A slight change in the rotation period can be defined,

Ω(t)=Ω01ϵ(t)Ω𝑡subscriptΩ01italic-ϵ𝑡\Omega\left(t\right)=\Omega_{0}\sqrt{1-\epsilon\left(t\right)}roman_Ω ( italic_t ) = roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT square-root start_ARG 1 - italic_ϵ ( italic_t ) end_ARG (18)

Dimensionless quantity ϵ(t)italic-ϵ𝑡\epsilon\left(t\right)italic_ϵ ( italic_t ) then becomes

ϵitalic-ϵ\displaystyle\epsilonitalic_ϵ (t)=2Ω0vgaNN2ρDM×\displaystyle\left(t\right)=\frac{2}{\Omega_{0}}vg_{aNN}\sqrt{2\rho_{DM}}\times( italic_t ) = divide start_ARG 2 end_ARG start_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG 2 italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG × (19)
(PcosΩ0t+QsinΩ0t+S)sin(ωt+ϕ0)𝑃subscriptΩ0𝑡𝑄subscriptΩ0𝑡𝑆𝜔𝑡subscriptitalic-ϕ0\displaystyle\left(P\cos{\Omega_{0}t}+Q\sin{\Omega_{0}t}+S\right)\sin\left(-% \omega t+\phi_{0}\right)( italic_P roman_cos roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t + italic_Q roman_sin roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t + italic_S ) roman_sin ( - italic_ω italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT )

Next, we turn to the determination of the pulsar’s timing residuals. Assume that dt𝑑𝑡dtitalic_d italic_t is the time pulsar takes to turn dθ𝑑𝜃d\thetaitalic_d italic_θ under the action of the axion field, dT𝑑𝑇dTitalic_d italic_T is the time pulsar takes to turn the same angle without the action of the axion field, we then can determine the timing residual of the pulsar per unit timeMaggiore (2018)::\colon:

z(t)=dtdTdT=dtdT1𝑧𝑡𝑑𝑡𝑑𝑇𝑑𝑇𝑑𝑡𝑑𝑇1z\left(t\right)=\frac{dt-dT}{dT}=\frac{dt}{dT}-1italic_z ( italic_t ) = divide start_ARG italic_d italic_t - italic_d italic_T end_ARG start_ARG italic_d italic_T end_ARG = divide start_ARG italic_d italic_t end_ARG start_ARG italic_d italic_T end_ARG - 1 (20)

According to the definition of the angular velocity, we can write dT𝑑𝑇dTitalic_d italic_T, dt𝑑𝑡dtitalic_d italic_t as::\colon:

dt=dθΩ(t)=dθΩ01ϵ(t)=dT1ϵ(t)𝑑𝑡𝑑𝜃Ω𝑡𝑑𝜃subscriptΩ01italic-ϵ𝑡𝑑𝑇1italic-ϵ𝑡\displaystyle dt=\frac{d\theta}{\Omega\left(t\right)}=\frac{d\theta}{\Omega_{0% }\sqrt{1-\epsilon\left(t\right)}}=\frac{dT}{\sqrt{1-\epsilon\left(t\right)}}italic_d italic_t = divide start_ARG italic_d italic_θ end_ARG start_ARG roman_Ω ( italic_t ) end_ARG = divide start_ARG italic_d italic_θ end_ARG start_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT square-root start_ARG 1 - italic_ϵ ( italic_t ) end_ARG end_ARG = divide start_ARG italic_d italic_T end_ARG start_ARG square-root start_ARG 1 - italic_ϵ ( italic_t ) end_ARG end_ARG (21)

dt/dT𝑑𝑡𝑑𝑇dt/dTitalic_d italic_t / italic_d italic_T can be written as:

dtdT=11ϵ(t)1+12ϵ(t)𝑑𝑡𝑑𝑇11italic-ϵ𝑡112italic-ϵ𝑡\frac{dt}{dT}=\frac{1}{\sqrt{1-\epsilon\left(t\right)}}\approx 1+\frac{1}{2}% \epsilon\left(t\right)divide start_ARG italic_d italic_t end_ARG start_ARG italic_d italic_T end_ARG = divide start_ARG 1 end_ARG start_ARG square-root start_ARG 1 - italic_ϵ ( italic_t ) end_ARG end_ARG ≈ 1 + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_ϵ ( italic_t ) (22)

The timing residual R(t)𝑅𝑡R(t)italic_R ( italic_t ) of the pulsar is defined to be the integration of z(t)𝑧𝑡z\left(t\right)italic_z ( italic_t ), measured with respect to a reference time t=0𝑡0t=0italic_t = 0:

R(t)𝑅𝑡\displaystyle R(t)italic_R ( italic_t ) =0t𝑑tz(t)absentsuperscriptsubscript0𝑡differential-d𝑡𝑧𝑡\displaystyle=\int_{0}^{t}dtz(t)= ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_t end_POSTSUPERSCRIPT italic_d italic_t italic_z ( italic_t ) (23)
=1Ω0vgaNN2ρDM[P(12(ωΩ0)cos((ωΩ0)t+ϕ0)\displaystyle=\frac{1}{\Omega_{0}}vg_{aNN}\sqrt{2\rho_{DM}}\left[P\left(\frac{% 1}{2\left(\omega-\Omega_{0}\right)}\cos\left(-\left(\omega-\Omega_{0}\right)t+% \phi_{0}\right)\right.\right.= divide start_ARG 1 end_ARG start_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG 2 italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG [ italic_P ( divide start_ARG 1 end_ARG start_ARG 2 ( italic_ω - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG roman_cos ( - ( italic_ω - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT )
ω(ωΩ0)(ω+Ω0)cosϕ0𝜔𝜔subscriptΩ0𝜔subscriptΩ0subscriptitalic-ϕ0\displaystyle-\frac{\omega}{\left(\omega-\Omega_{0}\right)\left(\omega+\Omega_% {0}\right)}\cos{\phi_{0}}- divide start_ARG italic_ω end_ARG start_ARG ( italic_ω - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ( italic_ω + roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG roman_cos italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT
+12(ω+Ω0)cos((ω+Ω0)t+ϕ0))\displaystyle\left.+\frac{1}{2\left(\omega+\Omega_{0}\right)}\cos\left(-\left(% \omega+\Omega_{0}\right)t+\phi_{0}\right)\right)+ divide start_ARG 1 end_ARG start_ARG 2 ( italic_ω + roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG roman_cos ( - ( italic_ω + roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) )
+Q(12(ωΩ0)sin((ωΩ0)t+ϕ0)\displaystyle+Q\left(\frac{1}{2\left(\omega-\Omega_{0}\right)}\sin\left(-\left% (\omega-\Omega_{0}\right)t+\phi_{0}\right)\right.+ italic_Q ( divide start_ARG 1 end_ARG start_ARG 2 ( italic_ω - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG roman_sin ( - ( italic_ω - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT )
Ω0(ωΩ0)(ω+Ω0)sinϕ0subscriptΩ0𝜔subscriptΩ0𝜔subscriptΩ0subscriptitalic-ϕ0\displaystyle-\frac{\Omega_{0}}{\left(\omega-\Omega_{0}\right)\left(\omega+% \Omega_{0}\right)}\sin{\phi_{0}}- divide start_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG ( italic_ω - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ( italic_ω + roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG roman_sin italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT
12(ω+Ω0)sin((ω+Ω0)t+ϕ0))\displaystyle\left.-\frac{1}{2\left(\omega+\Omega_{0}\right)}\sin\left(-\left(% \omega+\Omega_{0}\right)t+\phi_{0}\right)\right)- divide start_ARG 1 end_ARG start_ARG 2 ( italic_ω + roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG roman_sin ( - ( italic_ω + roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) )
+Sω(cos(ωt+ϕ0)cosϕ0)]\displaystyle\left.+\frac{S}{\omega}\left(\cos\left(-\omega t+\phi_{0}\right)-% \cos{\phi_{0}}\right)\right]+ divide start_ARG italic_S end_ARG start_ARG italic_ω end_ARG ( roman_cos ( - italic_ω italic_t + italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) - roman_cos italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ]

Notice that when we do the integral 0t𝑑tz(t)superscriptsubscript0𝑡differential-d𝑡𝑧𝑡\int_{0}^{t}dtz(t)∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_t end_POSTSUPERSCRIPT italic_d italic_t italic_z ( italic_t ), we do not include the case when ω=Ω0𝜔subscriptΩ0\omega=\Omega_{0}italic_ω = roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT which has to be computed separately, and the result is::\colon:

R𝑅\displaystyle Ritalic_R (t)=vgaNNρDM2Ω0(Psinϕ0Qcosϕ0)t𝑡𝑣subscript𝑔𝑎𝑁𝑁subscript𝜌𝐷𝑀2subscriptΩ0𝑃subscriptitalic-ϕ0𝑄subscriptitalic-ϕ0𝑡\displaystyle(t)=\frac{vg_{aNN}\sqrt{\rho_{DM}}}{\sqrt{2}\Omega_{0}}(P\sin{% \phi_{0}}-Q\cos{\phi_{0}})t( italic_t ) = divide start_ARG italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG end_ARG start_ARG square-root start_ARG 2 end_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ( italic_P roman_sin italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_Q roman_cos italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) italic_t (24)
+PvgaNNρDM22Ω02[cos(ϕ02Ω0t)cosϕ0]𝑃𝑣subscript𝑔𝑎𝑁𝑁subscript𝜌𝐷𝑀22superscriptsubscriptΩ02delimited-[]subscriptitalic-ϕ02subscriptΩ0𝑡subscriptitalic-ϕ0\displaystyle+P\frac{vg_{aNN}\sqrt{\rho_{DM}}}{2\sqrt{2}\Omega_{0}^{2}}[\cos(% \phi_{0}-2\Omega_{0}t)-\cos{\phi_{0}}]+ italic_P divide start_ARG italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG end_ARG start_ARG 2 square-root start_ARG 2 end_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ roman_cos ( italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - 2 roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t ) - roman_cos italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ]
+QvgaNNρDM2Ω02sinΩ0tcos(ϕ0Ω0t)𝑄𝑣subscript𝑔𝑎𝑁𝑁subscript𝜌𝐷𝑀2superscriptsubscriptΩ02subscriptΩ0𝑡subscriptitalic-ϕ0subscriptΩ0𝑡\displaystyle+Q\frac{vg_{aNN}\sqrt{\rho_{DM}}}{\sqrt{2}\Omega_{0}^{2}}\sin{% \Omega_{0}t}\cos(\phi_{0}-\Omega_{0}t)+ italic_Q divide start_ARG italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG end_ARG start_ARG square-root start_ARG 2 end_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG roman_sin roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t roman_cos ( italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t )
+SvgaNN2ρDMΩ02[cos(ϕ0Ω0t)cosϕ0]𝑆𝑣subscript𝑔𝑎𝑁𝑁2subscript𝜌𝐷𝑀superscriptsubscriptΩ02delimited-[]subscriptitalic-ϕ0subscriptΩ0𝑡subscriptitalic-ϕ0\displaystyle+S\frac{vg_{aNN}\sqrt{2\rho_{DM}}}{\Omega_{0}^{2}}[\cos(\phi_{0}-% \Omega_{0}t)-\cos{\phi_{0}}]+ italic_S divide start_ARG italic_v italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT square-root start_ARG 2 italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG end_ARG start_ARG roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ roman_cos ( italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_t ) - roman_cos italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ]

The first term in (24) has a linear dependence in t𝑡titalic_t and the timing residual will tend to infinity as time goes on. The physical meaning of this divergence is that the axion and the pulsar happen to have a resonance around 2.42 kHz(ma1011eV)2.42 kHzsubscript𝑚𝑎superscript1011eV2.42\text{\,kHz}\left(\frac{m_{a}}{10^{-11}\text{eV}}\right)2.42 kHz ( divide start_ARG italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 11 end_POSTSUPERSCRIPT eV end_ARG ). and the angular velocity of the pulsar is changed permanently which can be searched accordingly. Nonetheless, given the wide spectrum of the axion and different types of pulsars, it seems plausible for the neutron star to resonate with the dark matter.

If the dark matter is the QCD axion 242 MHz(ma106eV)242 MHzsubscript𝑚𝑎superscript106eV242\text{\,MHz}\left(\frac{m_{a}}{10^{-6}\text{eV}}\right)242 MHz ( divide start_ARG italic_m start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT eV end_ARG ) with pulsars Ω02π103s6.28×103Hzsimilar-tosubscriptΩ02𝜋superscript103ssimilar-to6.28superscript103Hz\Omega_{0}\sim\frac{2\pi}{10^{-3}\mathrm{s}}\sim 6.28\times 10^{3}\mathrm{Hz}roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ∼ divide start_ARG 2 italic_π end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT roman_s end_ARG ∼ 6.28 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT roman_Hz. Then we have ωΩ0much-greater-than𝜔subscriptΩ0\omega\gg\Omega_{0}italic_ω ≫ roman_Ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. Then, the axion background will induce a nutation of the magnetic axis of pulsars, which might be observed by closely looking at the neutron star’s movement. We will leave this and the resonant case to a future study.

For the other case, to determine the value of R(t)𝑅𝑡R(t)italic_R ( italic_t ), we choose:

Ω^0=(0,0,1)subscript^Ω0001\displaystyle\hat{\Omega}_{0}=\left(0,0,1\right)over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = ( 0 , 0 , 1 ) (25)
B^p(0)=(sinθ,0,cosθ)subscript^𝐵𝑝0𝜃0𝜃\displaystyle\hat{B}_{p}\left(0\right)=\left(\sin{\theta},0,\cos{\theta}\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) = ( roman_sin italic_θ , 0 , roman_cos italic_θ ) (26)
v^=(sinφ,0,cosφ)^𝑣𝜑0𝜑\displaystyle\hat{v}=\left(\sin{\varphi},0,\cos{\varphi}\right)over^ start_ARG italic_v end_ARG = ( roman_sin italic_φ , 0 , roman_cos italic_φ ) (27)

as the initial direction of the three axes.
Then we have

P=sinθsinφ𝑃𝜃𝜑\displaystyle P=\sin{\theta}\sin{\varphi}italic_P = roman_sin italic_θ roman_sin italic_φ (28)
Q=0𝑄0\displaystyle Q=0italic_Q = 0 (29)
S=cosθcosφ𝑆𝜃𝜑\displaystyle S=\cos{\theta}\cos{\varphi}italic_S = roman_cos italic_θ roman_cos italic_φ (30)

Consider a special case when B^p(0)subscript^𝐵𝑝0\hat{B}_{p}\left(0\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT ( 0 ) is parallel to v^^𝑣\hat{v}over^ start_ARG italic_v end_ARG, and is vertical to Ω^0::subscript^Ω0absent\hat{\Omega}_{0}\colonover^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT :

θ=φ=90𝜃𝜑superscript90\theta=\varphi=90^{\circ}italic_θ = italic_φ = 90 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT (31)

For the ultralight dark matter 2.42 nHz(m1023eV)2.42 nHz𝑚superscript1023eV2.42\text{\,nHz}\left(\frac{m}{10^{-23}\text{eV}}\right)2.42 nHz ( divide start_ARG italic_m end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 23 end_POSTSUPERSCRIPT eV end_ARG ), then we have

max𝑚𝑎𝑥\displaystyle maxitalic_m italic_a italic_x |R(t,x)|𝑅𝑡𝑥\displaystyle\left|R(t,\vec{x})\right|| italic_R ( italic_t , over→ start_ARG italic_x end_ARG ) | (32)
8.2×1017sec(gaNN109GeV1)(ρDM0.3GeV/cm3)absent8.2superscript1017secsubscript𝑔𝑎𝑁𝑁superscript109superscriptGeV1subscript𝜌𝐷𝑀0.3GeVsuperscriptcm3\displaystyle\approx 8.2\times 10^{-17}\mathrm{sec}\left(\frac{g_{aNN}}{10^{-9% }\mathrm{GeV^{-1}}}\right)\left(\sqrt{\frac{\rho_{DM}}{0.3\mathrm{GeV/cm^{3}}}% }\right)≈ 8.2 × 10 start_POSTSUPERSCRIPT - 17 end_POSTSUPERSCRIPT roman_sec ( divide start_ARG italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 9 end_POSTSUPERSCRIPT roman_GeV start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_ARG ) ( square-root start_ARG divide start_ARG italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG start_ARG 0.3 roman_GeV / roman_cm start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG end_ARG )

Consider another special case without losing generality:

θ=30φ=60formulae-sequence𝜃superscript30𝜑superscript60\theta=30^{\circ}\qquad\varphi=60^{\circ}italic_θ = 30 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT italic_φ = 60 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT (33)

The maximum value of the timing residual of the ultralight dark matter is::\colon:

max𝑚𝑎𝑥\displaystyle maxitalic_m italic_a italic_x |R(t,x)|𝑅𝑡𝑥\displaystyle\left|R(t,\vec{x})\right|| italic_R ( italic_t , over→ start_ARG italic_x end_ARG ) | (34)
1.9×107sec(gaNN1012GeV1)(ρDM0.3GeV/cm3)absent1.9superscript107secsubscript𝑔𝑎𝑁𝑁superscript1012superscriptGeV1subscript𝜌𝐷𝑀0.3GeVsuperscriptcm3\displaystyle\approx 1.9\times 10^{-7}\mathrm{sec}\left(\frac{g_{aNN}}{10^{-12% }\mathrm{GeV^{-1}}}\right)\left(\sqrt{\frac{\rho_{DM}}{0.3\mathrm{GeV/cm^{3}}}% }\right)≈ 1.9 × 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT roman_sec ( divide start_ARG italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 12 end_POSTSUPERSCRIPT roman_GeV start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_ARG ) ( square-root start_ARG divide start_ARG italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT end_ARG start_ARG 0.3 roman_GeV / roman_cm start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG end_ARG )

This is within the reach of current PTA observation. In Fig.1 we plot the timing residual induced by the axion together with the current Nanograv15 data. We can see that PTA data can already put constraints on the axion-nucleon coupling with some assumptions.

Refer to caption
Figure 1: Upper limits on the timing residual of pulsar generated by the oscillating dark matter, as a function of frequency/axion’s mass. The gray line is the recently released NanoGrav dataAgazie et al. (2023a). The blue, black, and red lines are the upper limits for the coupling constant gaNNsubscript𝑔𝑎𝑁𝑁g_{aNN}italic_g start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT equals 1012GeV1,1012.5GeV1,1013GeV1superscript1012superscriptGeV1superscript1012.5superscriptGeV1superscript1013superscriptGeV110^{-12}\mathrm{GeV^{-1}},10^{-12.5}\mathrm{GeV^{-1},10^{-13}\mathrm{GeV^{-1}}}10 start_POSTSUPERSCRIPT - 12 end_POSTSUPERSCRIPT roman_GeV start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , 10 start_POSTSUPERSCRIPT - 12.5 end_POSTSUPERSCRIPT roman_GeV start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , 10 start_POSTSUPERSCRIPT - 13 end_POSTSUPERSCRIPT roman_GeV start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT respectively.

III Timing residual the two-point correlation

The cross-correlation Hellings-Downs curve between the noise-independent data is essential to detect the gravitational wave backgroundHellings and Downs (1983); Armaleo et al. (2021). The cross-correlation between the timing residuals of the pulsar a and b can be calculated through za(t)zb(t)delimited-⟨⟩subscript𝑧𝑎𝑡subscript𝑧𝑏𝑡\left\langle z_{a}\left(t\right)z_{b}\left(t\right)\right\rangle⟨ italic_z start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ( italic_t ) italic_z start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ( italic_t ) ⟩, where delimited-⟨⟩\left\langle\cdots\right\rangle⟨ ⋯ ⟩ represents a time average or an ensemble average. Having the explicit form of the redshift z(t)𝑧𝑡z\left(t\right)italic_z ( italic_t ) in hand, we can then calculate,

za(t)zb(t)=1Ω0aΩ0bv2gaNN2ρDMSaSbcos(ϕ0aϕ0b)delimited-⟨⟩subscript𝑧𝑎𝑡subscript𝑧𝑏𝑡1subscriptΩ0𝑎subscriptΩ0𝑏superscript𝑣2subscriptsuperscript𝑔2𝑎𝑁𝑁subscript𝜌𝐷𝑀delimited-⟨⟩subscript𝑆𝑎subscript𝑆𝑏subscriptitalic-ϕ0𝑎subscriptitalic-ϕ0𝑏\left\langle z_{a}\left(t\right)z_{b}\left(t\right)\right\rangle=\frac{1}{% \Omega_{0a}\Omega_{0b}}v^{2}g^{2}_{aNN}\rho_{DM}\left\langle S_{a}S_{b}\right% \rangle\cos{\left(\phi_{0a}-\phi_{0b}\right)}⟨ italic_z start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ( italic_t ) italic_z start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ( italic_t ) ⟩ = divide start_ARG 1 end_ARG start_ARG roman_Ω start_POSTSUBSCRIPT 0 italic_a end_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT 0 italic_b end_POSTSUBSCRIPT end_ARG italic_v start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_g start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_a italic_N italic_N end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_D italic_M end_POSTSUBSCRIPT ⟨ italic_S start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_S start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ⟩ roman_cos ( italic_ϕ start_POSTSUBSCRIPT 0 italic_a end_POSTSUBSCRIPT - italic_ϕ start_POSTSUBSCRIPT 0 italic_b end_POSTSUBSCRIPT ) (35)

To simplify the calculation below, we then choose another frame that contains the pulsar location unit vector which points from the Earth towards the pulsar as the initial direction:

n^a=(0,0,1)subscript^𝑛𝑎001\displaystyle\hat{n}_{a}=\left(0,0,1\right)over^ start_ARG italic_n end_ARG start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT = ( 0 , 0 , 1 ) (36)
n^b=(sinζ,0,cosζ)subscript^𝑛𝑏𝜁0𝜁\displaystyle\hat{n}_{b}=\left(\sin{\zeta},0,\cos{\zeta}\right)over^ start_ARG italic_n end_ARG start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT = ( roman_sin italic_ζ , 0 , roman_cos italic_ζ ) (37)

Without losing generality we set the initial direction of the magnetic axis of the pulsars a and b as:

B^pa(0)=(0,0,1)subscript^𝐵𝑝𝑎0001\displaystyle\hat{B}_{pa}\left(0\right)=\left(0,0,-1\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p italic_a end_POSTSUBSCRIPT ( 0 ) = ( 0 , 0 , - 1 ) (38)
B^pb(0)=(sinζ,0,cosζ)subscript^𝐵𝑝𝑏0𝜁0𝜁\displaystyle\hat{B}_{pb}\left(0\right)=\left(-\sin{\zeta},0,-\cos{\zeta}\right)over^ start_ARG italic_B end_ARG start_POSTSUBSCRIPT italic_p italic_b end_POSTSUBSCRIPT ( 0 ) = ( - roman_sin italic_ζ , 0 , - roman_cos italic_ζ ) (39)

We also set the orientation of the axion’s velocity and the rotational axes of pulsars a and b:

v^=(sinαcosβ,sinαsinβ,cosα)^𝑣𝛼𝛽𝛼𝛽𝛼\displaystyle\hat{v}=\left(\sin{\alpha}\cos{\beta},\sin{\alpha}\sin{\beta},% \cos{\alpha}\right)over^ start_ARG italic_v end_ARG = ( roman_sin italic_α roman_cos italic_β , roman_sin italic_α roman_sin italic_β , roman_cos italic_α ) (40)
Ω^0a=(sinθacosφa,sinθasinφa,cosθa)subscript^Ω0𝑎subscript𝜃𝑎subscript𝜑𝑎subscript𝜃𝑎subscript𝜑𝑎subscript𝜃𝑎\displaystyle\hat{\Omega}_{0a}=\left(\sin{\theta_{a}}\cos{\varphi_{a}},\sin{% \theta_{a}}\sin{\varphi_{a}},\cos{\theta_{a}}\right)over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 italic_a end_POSTSUBSCRIPT = ( roman_sin italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT roman_cos italic_φ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT , roman_sin italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT roman_sin italic_φ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT , roman_cos italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) (41)
Ω^0b=(sinθbcosφb,sinθbsinφb,cosθb)subscript^Ω0𝑏subscript𝜃𝑏subscript𝜑𝑏subscript𝜃𝑏subscript𝜑𝑏subscript𝜃𝑏\displaystyle\hat{\Omega}_{0b}=\left(\sin{\theta_{b}}\cos{\varphi_{b}},\sin{% \theta_{b}}\sin{\varphi_{b}},\cos{\theta_{b}}\right)over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 italic_b end_POSTSUBSCRIPT = ( roman_sin italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT roman_cos italic_φ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT , roman_sin italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT roman_sin italic_φ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT , roman_cos italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ) (42)

According to Eq.(15), we can easily calculate S𝑆Sitalic_S of pulsar a and b:

Sa=cosθasubscript𝑆𝑎subscript𝜃𝑎\displaystyle S_{a}=-\cos{\theta_{a}}italic_S start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT = - roman_cos italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT (cosαcosθa+sinαcosβsinθacosφa\displaystyle(\cos{\alpha}\cos{\theta_{a}}+\sin{\alpha}\cos{\beta}\sin{\theta_% {a}}\cos{\varphi_{a}}( roman_cos italic_α roman_cos italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT + roman_sin italic_α roman_cos italic_β roman_sin italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT roman_cos italic_φ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT (43)
+sinαsinβsinθasinφa)\displaystyle+\sin{\alpha}\sin{\beta}\sin{\theta_{a}}\sin{\varphi_{a}})+ roman_sin italic_α roman_sin italic_β roman_sin italic_θ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT roman_sin italic_φ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT )
Sbsubscript𝑆𝑏\displaystyle S_{b}italic_S start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT =(cosζcosθbsinζsinθbcosφb)absent𝜁subscript𝜃𝑏𝜁subscript𝜃𝑏subscript𝜑𝑏\displaystyle=(-\cos{\zeta}\cos{\theta_{b}}-\sin{\zeta}\sin{\theta_{b}}\cos{% \varphi_{b}})= ( - roman_cos italic_ζ roman_cos italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT - roman_sin italic_ζ roman_sin italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT roman_cos italic_φ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ) (44)
×(cosαcosθb+sinαcosβsinθbcosφb\displaystyle\times(\cos{\alpha}\cos{\theta_{b}}+\sin{\alpha}\cos{\beta}\sin{% \theta_{b}}\cos{\varphi_{b}}× ( roman_cos italic_α roman_cos italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT + roman_sin italic_α roman_cos italic_β roman_sin italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT roman_cos italic_φ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT
+sinαsinβsinθbsinφb)\displaystyle+\sin{\alpha}\sin{\beta}\sin{\theta_{b}}\sin{\varphi_{b}})+ roman_sin italic_α roman_sin italic_β roman_sin italic_θ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT roman_sin italic_φ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT )

The pulsars we observe are distributed all over the sky, we can assume that their relative angle to the direction of the axion’s velocity is random. This means that the average of the cross-correlation signals among many pulsars effectively takes an average over the direction of the axion’s velocity. The rotational axis can also change its direction as a result of the motion of the galaxy, so we also take an average over the rotational axis and perform the integral over v^,Ω^0a,Ω^0b^𝑣subscript^Ω0𝑎subscript^Ω0𝑏\hat{v},\hat{\Omega}_{0a},\hat{\Omega}_{0b}over^ start_ARG italic_v end_ARG , over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 italic_a end_POSTSUBSCRIPT , over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 italic_b end_POSTSUBSCRIPT:

SaSb=dv^4πdΩ^0a4πdΩ^0b4πSaSb=127cosζdelimited-⟨⟩subscript𝑆𝑎subscript𝑆𝑏𝑑^𝑣4𝜋𝑑subscript^Ω0𝑎4𝜋𝑑subscript^Ω0𝑏4𝜋subscript𝑆𝑎subscript𝑆𝑏127𝜁\displaystyle\left\langle S_{a}S_{b}\right\rangle=\int\frac{d\hat{v}}{4\pi}% \frac{d\hat{\Omega}_{0a}}{4\pi}\frac{d\hat{\Omega}_{0b}}{4\pi}S_{a}S_{b}=\frac% {1}{27}\cos{\zeta}⟨ italic_S start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_S start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ⟩ = ∫ divide start_ARG italic_d over^ start_ARG italic_v end_ARG end_ARG start_ARG 4 italic_π end_ARG divide start_ARG italic_d over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 italic_a end_POSTSUBSCRIPT end_ARG start_ARG 4 italic_π end_ARG divide start_ARG italic_d over^ start_ARG roman_Ω end_ARG start_POSTSUBSCRIPT 0 italic_b end_POSTSUBSCRIPT end_ARG start_ARG 4 italic_π end_ARG italic_S start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_S start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 27 end_ARG roman_cos italic_ζ (45)

This curve is similar to the vector type modes in massive gravity types of studies Lee et al. (2008); Chamberlin and Siemens (2012); Cornish et al. (2018); Hellings and Downs (1983); Omiya et al. (2023) as expected, due to the one-dimensional directional nature of this modulation. We plot this axion-nucleon curve in Fig.2.

Refer to caption
Figure 2: The angle of separation dependence of the two-point correlation function, normalized to 1111 when θ=0𝜃0\theta=0italic_θ = 0. The dashed curve corresponds to the angular correlation of the timing residuals induced by the gravitational wave background(Hellings-Downs curve) and the solid line represents the correlation of the axion’s effect acting on the pulsars.

IV discussion

We consider a new detection possibility to employ the neutron stars as the axion wind detectors. Especially current PTA network observations can provide interesting sensitivity on the axion-nucleon couplings. This classical macroscopic detection mechanism is complementary to axion microscopic production and capture mechanism Wu et al. (2023c); Beznogov et al. (2018); Dev et al. (2024); Berlin et al. (2024); Chao et al. (2023); Gao et al. (2024); Caputo et al. (2020) to set the limits on the related coupling. Due to the complexity of the neutron star structure models, we expect noises from several different sources to get into the neutron star motion. Nevertheless, a close analysis of neutron star behavior might provide further insight into this type of coupling.

Acknowledgements

We thank Fei Gao for the discussion on the neutron stars. This work is supported by the National Key R&\&&D Program of China (grant 2023YFE0117200) and the National Natural Science Foundation of China (Nos. 12105013).

References