# Ultralight $(L_{\mu}-L_{\tau})$ vector dark matter interpretation of
NANOGrav observations

Debtosh Chowdhury<EMAIL_ADDRESS>Department of Physics, Indian Institute
of Technology Kanpur, Kanpur 208016, INDIA  Arpan Hait<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Kanpur, Kanpur 208016,
INDIA  Subhendra Mohanty<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Kanpur, Kanpur 208016, INDIA  Suraj Prakash
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Kanpur, Kanpur 208016, INDIA

The angular correlation of pulsar residuals observed by NANOGrav [1, 2, 3] and
other pulsar timing array (PTA) collaborations show evidence in support of the
Hellings-Downs correlation [4] expected from stochastic gravitational waves
(SGW). In this paper, we offer a non-gravitational wave explanation of the
observed pulsar timing correlations as caused by an ultra-light
$L_{\mu}-L_{\tau}$ gauge boson dark matter (ULDM). ULDM can affect the pulsar
correlations in two ways. The gravitational potential of vector ULDM gives
rise to a Shapiro time-delay of the pulsar signals and a non-trivial angular
correlation (as compared to the scalar ULDM case). In addition, if the pulsars
have a non-zero charge of the dark matter gauge group then the electric field
of the local dark matter causes an oscillation of the pulsar and a
corresponding Doppler shift of the pulsar signal. We point out that pulsars
carry a significant charge of muons and thus the $L_{\mu}-L_{\tau}$ vector
dark matter contributes to both the Doppler oscillations and the time-delay of
the pulsar signals. Our analysis shows that the NANOGrav data has a better fit
to the $L_{\mu}-L_{\tau}$ ULDM scenario compared to the SGW or the SGW with
Shapiro time-delay hypotheses.

††preprint:

##  I Introduction

Precision timing of the signals from array of pulsars with observations of
over a decade provides a method for the detection of stochastic gravitational
waves in the nano-Hertz frequency range. In their analysis of pulsar timing
residuals, the NANOGrav [1, 2, 3], EPTA [5], InPTA [6], PPTA [7] and CPTA [8]
experimental collaborations have observed a common red noise among the pulsars
whose angular correlation agrees with the Hellings-Downs curve [4] expected
from an isotropic stochastic background of gravitational waves. From the
15-year pulsar-timing data set correlated among 67 pulsars the NANOGrav
collaboration[1] finds evidence for a stochastic signal at $(3.5-4.0)\sigma$
statistical significance. The amplitude and spectral index of the observed
stochastic gravitational wave signal can be tested with theoretical templates
of predictions from inflation, phase transitions or cosmic defects [2]. The
observations of the NANOGrav collaboration also agree with prediction from the
cumulative stochastic gravitational waves (SGW) signals from the merger of
super-massive black holes ( $\sim 10^{6}M_{\odot}$ ) over cosmological times
[3]. In theories of gravity beyond Einstein’s GR, the scalar and vector modes
of the metric can propagate and these modes will give angular correlations
which deviate from the Hellings-Downs shape [9, 10], and the NANOGrav data [1]
have been used for constraining the non-Einsteinian modes of the metric
perturbations [11, 12].

The NANOGrav observations can also have explanations without gravitational
waves. Ultra-light vector dark matter (ULDM), with mass $\sim 10^{-23}{\rm
eV}$, induces displacements of the Earth as well as the pulsar that cause a
periodic signal in timing residuals with frequency $f=m/2\pi\simeq 2.4\times
10^{-9}{\rm Hz}$. This frequency falls into the sensitive region of PTA
experiments with observation time of decades [13, 14]. Examples of such
ultralight dark photons which can couple to the Earth and/or pulsars are gauge
bosons arising from $U(1)_{B}$ or $U(1)_{B-L}$ [15, 16, 17, 18, 19, 20, 21,
22].

The local Newtonian potential of Ultra-light dark matter can also cause a
gravitational frequency shift in the pulsar signals which affects the angular
correlations and whose shape depends on the spin of the ULDM [23, 24, 25, 26,
27]. The gravitational time-delay from ULDM will occur even if the DM
particles have no coupling to the Earth or Pulsar.

In this paper, we study the $L_{\mu}-L_{\tau}$ ULDM. In this case both the
effects mentioned above are present. Neutron stars carry substantial muons
charges and are therefore oscillated by the background dark matter [28, 29,
30]. The oscillation of the Earth term is absent as the Earth is expected to
have zero muon charge. In addition to oscillating the pulsars, the
$L_{\mu}-L_{\tau}$ ULDM also gives rise to the time-delay due to the
gravitational potential perturbation. We find that the taking both the time-
delay and the source oscillations together provides a better fit to the shape
of the angular correlation function, as found by the NANOGrav collaboration,
compared to the Hellings-Downs curved from SGW.

The $L_{\mu}-L_{\tau}$ gauge theory is of interest as it is an anomaly free
extension of the standard model without the need of adding extra families of
fermions for anomaly cancellation [31, 32, 33, 34] . The $L_{\mu}-L_{\tau}$
interaction is described by the Lagrangian

$\displaystyle{\cal
L}=-g^{\prime}Z^{\prime}_{\alpha}\left(\bar{\mu}\gamma^{\alpha}\mu-\bar{\tau}\gamma^{\alpha}\tau+\bar{\nu}_{\mu}\gamma_{\alpha}\nu_{\mu}-\bar{\nu}_{\tau}\gamma_{\alpha}\nu_{\tau}\right).$
(1)

These interactions can give rise to long range forces between leptons which
may be probed via neutrino oscillations [35, 36, 37, 38, 39, 40, 41] or binary
pulsar timings [42, 43]. The $L_{\mu}-L_{\tau}$ gauge bosons if they are
ultra-light can also serve as the fuzzy dark matter of the universe [45, 46,
48].

Since planets and stars do not carry any $L_{\mu}-L_{\tau}$ charge, the long
range gauged forces associated with this charge cannot be probed from
planetary or stellar dynamics unlike other anomaly free combinations of
charges like $L_{e}-L_{\tau/\mu}$ or $B-L$ which are sourced by the electron
or baryon content of these bodies [47]. However, neutron stars carry a
substantial number of muons which can be a source of ultra-light gauge boson
radiation from binary neutron stars and these can be radiated in binary stars
and therefore their mass and couplings can be constrained from binary pulsar
timings [42, 43]. The best fit values deduced through our analysis, of the
mass and couplings of the $L_{\mu}-L_{\tau}$ vector field, fall within the
allowed parameter space based on the measured loss in the time periods of such
binary pulsars.

We introduce the vector dark matter model and evaluate the correlations of
pulsar timing residuals in section II. We then conduct a detailed comparison
of the angular correlations due to different sources, i.e. the stochastic
gravitational wave background, the gravitational potential of the dark matter
and from pulsar oscillations due to their coupling with the dark matter in
section III. We also identify the best fit to data, among the different
theoretical models, through a $\chi^{2}$ minimization. We then translate the
results of the fitting with data into constraints on the dark matter mass and
coupling parameters in section IV. Finally, we summarize the analysis and
present our conclusions in section V.

##  II Pulsar period residual from Proca dark matter

We consider the case when neutron stars carry non-zero charges of the gauge
group of the Proca field. This is true if the gauge group is $B$ or $B-L$.
This is also true for gauge groups like $L_{\mu}-L_{\tau}$ as neutron stars
carry large number of muons whose decay is prevented by Fermi-blocking of the
outgoing electrons due to electron degeneracy [42, 43].

###  II.1 Proca dark matter

The Lagrangian for a massive vector field is given by

$\displaystyle{\cal
L}=-\cfrac{1}{4}\,F_{\mu\nu}\,F^{\mu\nu}+\cfrac{1}{2}\,m^{2}A_{\mu}\,A^{\mu}+g\,J_{\mu}\,A^{\mu}\,.$
(2)

From this we obtain the equation of motion for the vector field
$A_{\mu}(\vec{x},t)$ in free space ($J^{\mu}=0$) given by

$\displaystyle\partial_{\mu}\,F^{\mu\nu}+m^{2}A^{\nu}=0,$ (3)

which reduces to the following Eqs.,

$\displaystyle\left(\partial_{t}^{2}-\nabla^{2}+m^{2}\right)A^{\mu}=0\,,$ (4)
$\displaystyle\partial_{t}\,A^{t}-\nabla\cdot\vec{A}=0\,.$ (5)

From Eq. (5), we see that $A^{t}=(k_{i}/m)A^{i}$. If the massive vector boson
constitutes dark matter than the momenta $k_{i}\sim mv_{i}$ and where the dark
matter velocity in the galaxy $v^{i}\sim 10^{-3}$. This implies that
$A_{t}=v_{i}A^{i}$. Thus the spatial components dominate $A_{i}\gg A_{t}$ and
we can take $A_{\mu}\simeq(0,\vec{A})$ for the dark matter vector field. The
energy density is given by the time component of the stress tensor,

$\displaystyle
T_{00}\simeq\cfrac{1}{2}\,m^{2}\,|\vec{A}|^{2}=\rho_{\text{dm}},$ (6)

where the dark matter density in the galaxy is $\rho_{\text{dm}}=0.4\,{\rm
GeV/cm^{3}}$. The vector field for the Proca dark matter can be written as

$\displaystyle\vec{A}(\vec{x},t)=\cfrac{\sqrt{2\rho_{\text{dm}}}}{m}\,\cos(mt+\vec{k}\cdot\vec{x}),$
(7)

with $|\vec{k}|=m\,v$.

###  II.2 Angular correlations of pulsar timing residuals

A neutron star of mass $M_{a}$ carrying the dark photon charge $Q_{a}$ in the
Proca dark matter will experience a time dependent force due to the dark
’electric field’ given by

$\displaystyle\vec{F}=\sum_{A=1,2}Q_{a}\,m\,\vec{\epsilon}_{A}(\vec{n})\,\cfrac{\sqrt{2\rho_{\text{dm}}}}{m}\,\cos(mt+\vec{k}\cdot\vec{x}),$
(8)

where $\vec{\epsilon}_{A}(\vec{n})$ is the polarisation vector of gauge field.
In time $\delta t$ the neutron star has a perturbation in the position given
by

$\displaystyle\delta\vec{x}=-\sum_{A=1,2}\,\vec{\epsilon}_{A}(\vec{n})\,\cfrac{Q_{a}}{M_{a}\,m}\,\cfrac{\sqrt{2\rho_{\text{dm}}}}{m}\,\cos(mt+\vec{k}\cdot\vec{x}).$
(9)

Due to the perturbation the time period shift $\Delta T_{a}$ of a pulsar
period $T_{a}$ which is located in the direction $\hat{n}_{a}$ and distance
$ct_{a}$ from the earth will have a period shift given by

$\displaystyle\Delta T_{a}(t)$ $\displaystyle=$
$\displaystyle\hat{n}_{a}\cdot\left[\,\delta\vec{x}(t-t_{a}+T_{a})-\delta\vec{x}(t-t_{a})\,\right]$
$\displaystyle\simeq$
$\displaystyle\sum_{A=1,2}\,T_{a}\,\hat{n}_{a}\cdot\vec{\epsilon}_{A}(\vec{n})\,\cfrac{Q_{a}}{M_{a}}\,\cfrac{\sqrt{2\rho_{\text{dm}}}}{m}\,\cos\left[m(t-t_{a})\right].$

Defining the pulsar period redshift and the pulsar timing residual as

$\displaystyle z_{a}(t)\,\equiv\,-\cfrac{\Delta
T_{a}}{T_{a}}\quad\text{and}\quad
R_{a}(t)=\int_{0}^{t}\,z_{a}(t^{\prime})\,dt^{\prime},$ (11)

we can compute the two point correlation of period residuals of two pulsars
with angular separation $\theta_{ab}$ as

$\displaystyle\langle
R_{a}(t)R_{b}(t+\tau)\rangle=\cfrac{Q_{a}Q_{b}}{M_{a}M_{b}}\,\cfrac{\rho_{\text{dm}}}{m^{4}}\,\cfrac{2}{3}\cos\theta_{ab}\cos(m\tau).$
(12)

Here the angular brackets denote averaging over time and the directions of the
DM Proca field propagation direction $\vec{n}$. We have made use of the
relation

$\displaystyle\int\,\cfrac{d^{2}\vec{n}}{4\pi}\,\sum_{A}\,\epsilon_{A}^{i}(\vec{n})\,\epsilon_{A}^{j}(\vec{n})\,n^{i}_{a}\,n^{j}_{b}=\cfrac{2}{3}\,\cos\theta_{ab}.$
(13)

which we have derived in appendix B.

##  III Stochastic Gravitational Waves vs Dark Matter

###  III.1 Angular correlations due to distinct sources

We compare three different theoretical models aimed at explaining the
experimentally reported angular correlations:

Case 1: When the stochastic gravitational wave background (SGWB) is the only
source of pulsar oscillations, the angular correlation of the timing residual
can be written as a harmonic sum of the SGW frequency$f_{i}$ [1]:

$\displaystyle\xi^{(1)}_{ab}(\theta_{ab})=\sum_{f_{i}}\Phi_{\rm
GW}(f_{i})\Gamma_{\text{HD}}(\theta_{ab})\cos(f_{i}\tau),$ (14)

where the amplitude due to stochstic gravitational waves can be written as [1]

$\displaystyle\Phi_{\rm GW}(f_{i})=\frac{A_{\rm
GW}^{2}}{12\pi^{2}}\frac{1}{T_{\rm obs}}\left(\frac{f_{i}}{f_{\rm
ref}}\right)^{-\gamma}f_{\rm ref}^{-3}\,,$ (15)

where $A_{\rm GW}$ is the average amplitude of the SGW and $\gamma$ is the
spectral index which depends on the nature of the source and for binary black-
hole mergers $\gamma=13/3$ [44].

The angular correlation for SGW is given by the Hellings-Downs function [4]

$\displaystyle\Gamma_{\rm HD}(\theta_{ab})$ $\displaystyle=$
$\displaystyle\cfrac{1}{2}+\cfrac{3}{2}\,\cfrac{(1-\cos\theta_{ab})}{2}\,\ln\cfrac{(1-\cos\theta_{ab})}{2}$
(16) $\displaystyle\qquad-\cfrac{1}{4}\cfrac{(1-\cos\theta_{ab})}{2}.$

Case 2: If, in addition to the pulsar oscillations due to the SGWB, time-delay
due to the gravitational potential of vector dark matter (the Shapiro time-
delay) is also taken into account, then the total angular correlation can be
obtained as [27, 13]:

$\displaystyle\xi^{(2)}_{ab}(\theta_{ab})=\Phi_{\text{VDM}}^{\text{GP}}\,\Gamma_{\text{VDM}}^{\text{GP}}(\theta_{ab})+\Phi_{\text{GW}}\,\Gamma_{\text{HD}}(\theta_{ab}).$
(17)

where $\Phi_{\text{VDM}}^{\text{GP}}$ is a function of the dark matter mass
and density and

$\displaystyle\Gamma_{\text{VDM}}^{\text{GP}}(\theta_{ab})$ $\displaystyle=$
$\displaystyle\cfrac{5}{138}\,P_{0}\left(\cos\theta_{ab}\right)+\cfrac{64}{138}\,P_{2}\left(\cos\theta_{ab}\right)$
(18) $\displaystyle=$
$\displaystyle\cfrac{1}{46}\,\left(32\cos^{2}\theta_{ab}-9\right).$

We define the ratio
$\beta=\cfrac{\Phi_{\text{GW}}}{\Phi_{\text{VDM}}^{\text{GP}}}$ and normalize
the angular correlation so as to ensure that $\xi^{(2)}_{ab}(0)=\cfrac{1}{2}$
irrespective of the choice of $\beta$. Therefore,

$\displaystyle\xi^{(2)}_{ab}(\theta_{ab})$ $\displaystyle=$
$\displaystyle\cfrac{\Phi_{\text{VDM}}^{\text{GP}}\,\Gamma_{\text{VDM}}^{\text{GP}}(\theta_{ab})+\Phi_{\text{GW}}\,\Gamma_{\text{HD}}(\theta_{ab})}{\Phi_{\text{VDM}}^{\text{GP}}+\Phi_{\text{GW}}}$
(19) $\displaystyle=$
$\displaystyle\cfrac{1}{1+\beta}\left[\Gamma^{\text{GP}}_{\text{VDM}}(\theta_{ab})+\beta\,\Gamma_{\text{HD}}(\theta_{ab})\right].$

Case 3: As an alternative to the effect of the SGWB, if the source of pulsar
oscillations is their coupling with the Proca field, on account of their muon
content (see appendix A), then also taking into account the Shapiro time-
delay, the total angular correlation has the form:

$\displaystyle\xi^{(3)}_{ab}(\theta_{ab})=\Phi_{\text{VDM}}^{\text{GP}}\,\Gamma_{\text{VDM}}^{\text{GP}}(\theta_{ab})+\Phi_{\text{VDM}}\,\Gamma_{\text{VDM}}(\theta_{ab}).$
(20)

Here, $\Phi_{\text{VDM}}^{\text{GP}}$ and $\Phi_{\text{VDM}}$ are functions of
the DM and pulsar parameters,

$\displaystyle\Phi_{\text{VDM}}$ $\displaystyle=$
$\displaystyle\cfrac{4}{3}\,\cfrac{Q_{a}\,Q_{b}}{M_{a}\,M_{b}}\,\cfrac{\rho_{\text{dm}}}{m^{4}},$
$\displaystyle\Phi^{\text{GP}}_{\text{VDM}}$ $\displaystyle=$
$\displaystyle\cfrac{69}{5\,(2m)^{2}}\,\cfrac{\pi^{2}G^{2}}{36}\,\cfrac{4\,\rho_{\text{dm}}^{2}}{m^{4}}.$
(21)

The angular correlation due to Doppler shift of pulsar signal from pulsar
oscillations in the VDM background is given in Eq. (13),

$\displaystyle\Gamma_{\text{VDM}}(\theta_{ab})$ $\displaystyle=$
$\displaystyle\cfrac{1}{2}\cos\theta_{ab}.$ (22)

where we have absorbed a factor of $4/3$ in the amplitude so as to follow the
HD normalisation, see Eq. (16), whereby we define the shape parameters such
that $\Gamma(\theta_{ab}=0)=\frac{1}{2}$.

We define $\alpha=\cfrac{\Phi_{\text{VDM}}}{\Phi_{\text{VDM}}^{\text{GP}}}$
and normalize the angular correlation so as to ensure that
$\xi^{(3)}_{ab}(0)=\cfrac{1}{2}$, irrespective of the choice of $\alpha$.
Hence, we rewrite Eq. (20) as

$\displaystyle\xi_{ab}^{(3)}(\theta_{ab})$ $\displaystyle=$
$\displaystyle\cfrac{\Phi_{\text{VDM}}^{\text{GP}}\,\Gamma_{\text{VDM}}^{\text{GP}}(\theta_{ab})+\Phi_{\text{VDM}}\,\Gamma_{\text{VDM}}(\theta_{ab})}{\Phi_{\text{VDM}}^{\text{GP}}+\Phi_{\text{VDM}}}$
(23) $\displaystyle=$
$\displaystyle\cfrac{1}{1+\alpha}\left[\Gamma_{\text{VDM}}^{\text{GP}}(\theta_{ab})+\alpha\,\Gamma_{\text{VDM}}(\theta_{ab})\right].$

We ascertain the best fit values of $\beta$ and $\alpha$ by means of
$\chi^{2}$ minimization in the next section.

###  III.2 Fit with NANOGrav data

We have conducted separate fits with respect to the amplitude and shape of the
reported signal. This is facilitated by the fact that the angular correlations
have been normalized as $\xi^{(i)}_{ab}(0)=1/2$, $i=1,2,3$, independent of the
amplitude corresponding to the SGW, the Shapiro time-delay or the Doppler
shift of pulsar signals on account of the VDM background.

The $\chi^{2}$ function corresponding to the angular correlations can be
constructed in terms of the experimentally reported correlations
$\xi^{\text{exp}}_{ab}(\theta_{i})$, the theoretically computed correlations
$\xi^{\text{exp}}_{ab}(\theta_{i},\vec{\delta})$ and the experimentally
reported error $\sigma^{\text{exp}}(\theta_{i})$ as

$\displaystyle\chi^{2}(\vec{\delta})=\sum_{i}\left(\cfrac{\xi^{\text{exp}}_{ab}(\theta_{i})-\xi^{\text{th}}_{ab}(\theta_{i},\vec{\delta})}{\sigma^{\text{exp}}(\theta_{i})}\right)^{2},$
(24)

where $\theta_{i}$ designate the data points and $\vec{\delta}$ denotes the
set of parameters to be fitted. For the three cases considered in the previous
section, we find that,

  1. 1.

For pulsar oscillations due to the SGW, the Hellings-Downs function describes
the angular correlations. As evident from Eq. (16), there are no unknown
parameters and one can simply evaluate the $\chi^{2}$ function as

$\displaystyle\chi^{2}=\sum_{i}\left(\cfrac{\xi^{\text{exp}}_{ab}(\theta_{i})-\xi^{(1)}_{ab}(\theta_{i})}{\sigma^{\text{exp}}(\theta_{i})}\right)^{2}=94.303.$
(25)

  2. 2. Figure 1: A plot of $\chi^{2}$ as a function of the fraction $\beta$ for the model that takes into account the combined effects of the stochastic gravitational background and the gravitational potential of vector dark matter. The highlighted region corresponds to $\beta\in(0,1)$.

For the combined effect of oscillations due to the SGW and the Shapiro time-
delay in the signals emanating from the pulsars, owing to the gravitational
potential of the VDM, the $\chi^{2}$ is a function of only one unknown
variable $\beta$, i.e.,

$\displaystyle\chi^{2}(\beta)$ $\displaystyle=$
$\displaystyle\sum_{i}\left(\cfrac{\xi^{\text{exp}}_{ab}(\theta_{i})-\xi^{(2)}_{ab}(\theta_{i},\beta)}{\sigma^{\text{exp}}(\theta_{i})}\right)^{2}$
(26) $\displaystyle=$
$\displaystyle\cfrac{5.068+35.052\,\beta+94.303\,\beta^{2}}{(1+\beta)^{2}}.$

The minima is obtained for $\beta^{*}\rightarrow 0$ with
$\chi^{2}_{\text{min}}=5.068$, which is a notable improvement over the
$\chi^{2}$ corresponding to the Hellings-Downs correlation. The features of
$\chi^{2}(\beta)$ have been shown in Fig. 1.

Extrapolating the $\chi^{2}(\beta)$ curve beyond $\beta=1$, it can be seen
that the function keeps on increasing.

  3. 3. Figure 2: A plot of $\chi^{2}$ as a function of the fraction $\alpha$ for the theoretical model that takes into account the combined effects of the gravitational potential of Proca dark matter and its coupling to the pulsars. Here, $\alpha^{*}$ corresponds to the value that yields the minimum value of $\chi^{2}_{\text{min}}$ and the region highlighted in red corresponds to the $1\sigma$ confidence interval.

When the source of the oscillations is the coupling of the pulsars with an
$L_{\mu}-L_{\tau}$ ULDM background and also accounting for the Shapiro time-
delay, the $\chi^{2}$ can be expressed as a function of the variable $\alpha$,
i.e.,

$\displaystyle\chi^{2}(\alpha)$ $\displaystyle=$
$\displaystyle\sum_{i}\left(\cfrac{\xi^{\text{exp}}_{ab}(\theta_{i})-\xi^{(3)}_{ab}(\theta_{i},\alpha)}{\sigma^{\text{exp}}(\theta_{i})}\right)^{2}$
(27) $\displaystyle=$
$\displaystyle\cfrac{5.068+1.706\,\alpha+26.612\,\alpha^{2}}{(1+\alpha)^{2}}.$

The minima for this function is obtained for $\alpha^{*}=0.16365$ with
$\chi^{2}_{\text{min}}\approx 4.47$, which improves over both the two previous
cases. The features of $\chi^{2}(\alpha)$ along with the minima have been
distinctly highlighted in Fig. 2, where we have also highlighted the range of
$\alpha$ values that correspond to the $1\sigma$ confidence interval. Since
negative values of $\alpha$ do not carry any physical meaning, the lower limit
of the $1\sigma$ interval is truncated at $\alpha=0$.

Figure 3:  A comparison of the three theoretical models describing the angular
correlation of pulsars through: (1) The Hellings-Downs angular correlation
with $\chi^{2}=94.303$, (2) The curve corresponding to the best fit scenario
involving the Hellings-Downs correlation along with the angular correlation
due to the time-delay caused by the gravitational potential of Proca dark
matter, with $\chi^{2}_{\text{min}}=5.068$, (3) Angular correlation due to the
gravitational potential as well as the coupling of the Proca dark matter with
the pulsars with $\chi^{2}_{\text{min}}=4.47$, and (3). We have also shown the
experimentally reported data along with the error bars [1].

The angular correlations corresponding to these three cases have been plotted
against the experimentally reported
$(\theta_{i},\xi^{\text{exp}}_{ab}(\theta_{i}))$ values [1] in Fig. 3.

##  IV Constraints on the Parameters

Figure 4: The allowed range of $g$ vs $m$ values based on constraints from the
amplitude as well as shape of the reported signal [1]. The solid lines
correspond to contours obtained from Eq. (30) after substituting $\alpha$
values corresponding to the boundaries of the $1\sigma$ region. The vertical
dashed lines highlight the constraints on the dark matter mass from the
amplitude of the signal. An enlarged view of the overlap region, along with
the best fit-point $(m=3.825\times 10^{-23}\text{ eV},g=3.607\times 10^{-26})$
corresponding to $\alpha=0.16365$ has been shown as an inset. The $(g,m)$
parameter space above the dashed-dotted curve is ruled out based on the
measurements of binary pulsar timings [42].

The ratio $\alpha$ can be related to the DM parameters, i.e. the coupling $g$
and the mass $m$ by noting that from Eq. (III.1),

$\displaystyle\alpha=\cfrac{\Phi_{\text{VDM}}}{\Phi_{\text{VDM}}^{\text{GP}}}=\cfrac{80}{23\pi^{2}}\,\cfrac{Q_{a}\,Q_{b}}{M_{a}\,M_{b}\,G^{2}}\,\cfrac{m^{2}}{\rho_{\text{dm}}}.$
(28)

After substituting the following numerical values [42],

$\displaystyle N=10^{55},\quad M_{a}=M_{b}=10^{57}\,\text{GeV},$
$\displaystyle\rho_{\text{dm}}=3.1\times 10^{-42}\,\text{GeV}^{4}.$ (29)

we obtain

$\displaystyle\alpha=2.518\times 10^{95}\,g^{2}\,m^{2}.$ (30)

Here $m$ is in units of electron-Volt (eV). Substituting the best fit value of
$\alpha$, yields a contour in the $g-m$ plane. The $1\sigma$ allowed region
defined by $\alpha\in(0,0.477]$ in Fig. 2 has been translated to upper and
lower111For the lower bound, we have chosen $\alpha=10^{-10}$, since
$\alpha=0$ implies either $g=0$ or $m=0$ from Eq. (30). bounds in the $g-m$
plane and highlighted by solid lines in Fig. 4.

We also obtain constraints on the mass of the vector dark matter by
calibrating the amplitude, corresponding to the combined effect of the Shapiro
time-delay and the oscillations of the pulsars in the VDM background, against
the amplitude of the SGW signal reported by NANOGrav, i.e.,

$\displaystyle(1+\alpha)\,\Phi^{\text{GP}}_{\text{VDM}}=\Phi_{\text{GW}}(f_{i}).$
(31)

Using the expressions for $\Phi_{\text{GW}}(f_{i})$,
$\Phi^{\text{GP}}_{\text{VDM}}$ given in Eqs. (15) and (III.1) respectively,
we substitute the following values,

$\displaystyle A_{\text{GW}}$ $\displaystyle=$ $\displaystyle 2.4\times
10^{-15},\qquad\gamma=13/3,$ $\displaystyle f_{\text{ref}}$ $\displaystyle=$
$\displaystyle 1\text{yr}^{-1},\,\,\quad\qquad T_{\text{obs}}=15\text{yrs}.$
(32)

Then, setting $f_{i}=(m/\pi)$ eV and solving for $m$ in Eq. (31) gives us
constant values of $m$ corresponding to a fixed choice of $\alpha$. These
constraints are highlighted using the vertical dashed lines in Fig. 4. The
parallelogram formed by the intersection of the upper and lower bounds
obtained from the correlation as well as the amplitude defines the $1\sigma$
allowed region in the $g-m$ plane. We note that this narrow range of $(g,m)$
values falls neatly within the broader experimentally allowed region based on
the measurements of the time periods of binary pulsars [42]. The boundary
separating the allowed and disallowed regions have been demarcated by the
crimson dashed-dotted line in Fig. 4, and the area below the line corresponds
to the permitted parameter space.

##  V Conclusions

In this paper, we proposed a non-gravitational explanation of NANOGrav pulsar
timing observations [1, 2]. We proposed that $L_{\mu}-L_{\tau}$ gauge bosons
as the ultralight dark matter, which can have two effects on pulsar timings
observed on earth. Firstly, pulsars carry about $\sim 10^{55}$ muons,
therefore they will oscillate with a frequency $f=m/(2\pi)$ due to the
$L_{\mu}-L_{\tau}$ ‘electric field’ of the dark matter. This will give rise to
an angular correlation between pulsar residuals
$\propto\frac{1}{2}\cos\theta_{ab}$ as compared to pulsar oscillations due to
stochastic gravitational waves which give the well-known Hellings-Downs
angular correlations [4]. The second effect of ultra-light vector dark matter
is to cause a gravitational time-delay which also contributes to the angular
correlation of pulsar residuals [13, 27]. We found that the combination of
these two effects of $L_{\mu}-L_{\tau}$ dark matter gives a better fit to the
angular correlation observations of NANOGrav [1], than the standard stochastic
GW explanation or the SGW combined with time-delay explanations. The mass and
couplings of the $L_{\mu}-L_{\tau}$ gauge bosons which give the best fit to
NANOGrav are not ruled out by other astrophysical constraints [42]. Future
observations of the angular-correlation of pulsar timings may be able to
enhance the statistical significance of the $L_{\mu}-L_{\tau}$ dark matter
explanation or rule it out.

###### Acknowledgements.

S.P. is supported by the MHRD, Government of India, under the Prime Minister’s
Research Fellows (PMRF) Scheme, 2020. A.H. would like to thank the MHRD,
Government of India for the research fellowship. This research of D.C. is
supported by an initiation grant IITK/PHY/2019413 at IIT Kanpur and by a DST-
SERB grant SERB/CRG/2021/007579.

##  Appendix A $L_{\mu}-L_{\tau}$ charge of pulsars

In the neutron stars, electrons are degenerate and there is a beta-equilibrium
in the processes $n\rightarrow p+e^{-}+\bar{\nu}_{e}$, $n\rightarrow
p+\mu^{-}+\bar{\nu}_{\mu}$ and $\mu^{-}\rightarrow
e^{-}+\bar{\nu}_{e}+\nu_{\mu}$. The charge neutrality condition constrains the
number densities of $p,\,e^{-}$ and $\mu^{-}$ to obey the relation
$n_{p}=n_{e}+n_{\mu}$. In neutron stars electrons are degenerate and have a
Fermi momentum $k_{f}=(3\pi^{2}n_{e})^{1/2}$ greater than the mass of muon
$m_{\mu}\sim 106\,{\rm MeV}$. There is a beta blocking of the outgoing
electrons in the muon decay and muons are stable inside neutron stars. The
equilibrium of the beta-decay processes in the neutron stars gives the
relations between the chemical equilibrium of the different species the
relations $\mu_{n}=\mu_{p}+\mu_{e}=\mu_{p}+\mu_{\mu}$ which implies that
$\mu_{e}=\mu_{\mu}$. The chemical potential of the degenerate electrons is the
fermi-energy given by

$\displaystyle\mu_{e}={\epsilon_{e}}_{f}={k_{e}}_{f}=\left(3\pi^{2}n_{e}\right)^{1/3}\,.$
(33)

The muon-chemical potential written in terms of muon number density is

$\displaystyle\mu_{\mu}={\epsilon_{\mu}}_{f}$ $\displaystyle=$
$\displaystyle\left(m_{\mu}^{2}+{{k_{\mu}}_{f}}^{2}\right)^{1/2}$ (34)
$\displaystyle=$
$\displaystyle\left(m_{\mu}^{2}+\left(3\pi^{2}n_{m}u\right)^{2/3}\right)^{1/2}\,.$

Equating Eqs. (33) and (34) we can solve for $\nu_{\mu}$ in terms of $\nu_{e}$
to obtain

$\displaystyle
n_{\mu}=\cfrac{1}{3\pi^{2}}\left[\left(3\pi^{2}n_{e}\right)^{2/3}-m_{\mu}^{2}\right]^{3/2}.$
(35)

Model calculations give the electron fraction $Y_{e}=n_{e}/n=0.052$ and the
nucleon density in neutron stars $n=0.238{\,\rm fm^{-3}}$ [28]. This gives the
number density of electrons in neutron stars $n_{e}=Y_{e}\,n=0.01237\,{\rm
fm^{-3}}=1.237\times 10^{37}\,{\rm cm^{-3}}$. Using this in Eq. (35) we obtain
the number density of muons in neutron stars to be $n_{\mu}=3.63\times
10^{36}{\,\rm cm^{-3}}$ which implies that a typical neutron star with of
$R=10{\,\rm km}$ radius will carry a total muon number

$\displaystyle N_{\mu}=n_{\mu}\left(\cfrac{4\pi}{3}R^{3}\right)=1.55\times
10^{55}\,.$ (36)

The $L_{\mu}-L_{\tau}$ charge of a neutron star is then the coupling constant
times the number of muons, $Q=gN_{\mu}$.

##  Appendix B Angular integration over polarizations

The angular integration in Eq. (13) can be performed by choosing an arbitrary
$\hat{n}$ = $(\sin\theta\,\cos\phi,\sin\theta\,\sin\phi,\cos\theta)$ and
$\epsilon_{1,2}(\vec{n})$ as follows:

$\displaystyle\epsilon_{1}(\vec{n})=\begin{pmatrix}-\sin\phi\\\ \cos\phi\\\
0\end{pmatrix},\quad\epsilon_{2}(\vec{n})=\begin{pmatrix}\cos\theta\,\cos\phi\\\
\cos\theta\,\sin\phi\\\ -\sin\theta\end{pmatrix}.$ (37)

This choice of $\epsilon_{1,2}$ ensures that

$\displaystyle\vec{n}\cdot\epsilon_{1}(\vec{n})=\vec{n}\cdot\epsilon_{2}(\vec{n})=\epsilon_{1}(\vec{n})\cdot\epsilon_{2}(\vec{n})=0.$
(38)

For the pulsars we can choose one of them to be oriented along the $z$-axis,
$\vec{n}_{a}=(0,0,1)$, and the other along the direction
$\vec{n}_{b}=(\sin\theta_{ab}\,\cos\phi_{ab},\,\sin\theta_{ab}\,\sin\phi_{ab},\,\cos\theta_{ab})$.
With these choices, the integrand in Eq. (13) becomes

$\displaystyle\sum_{A=1,2}\,\epsilon_{A}^{i}(\vec{n})\,\epsilon_{A}^{j}(\vec{n})\,n^{i}_{a}\,n^{j}_{b}=\sin^{2}\theta\,\cos\theta_{ab}$
$\displaystyle\qquad\qquad-\sin\theta\,\cos\theta\,\sin\theta_{ab}\,\cos(\phi-\phi_{ab}).$
(39)

Integrating first with respect to $\sin\theta\,d\theta$ within the limits
$\theta\in[0,\pi]$ gives the result $\cfrac{4}{3}\,\cos\theta_{ab}$. An
overall factor of $2\pi$ results from the integration with respect to $\phi$.
Thus, we arrive at

$\displaystyle\cfrac{1}{4\pi}\int\sin\theta\,d\theta\,d\phi\,\sum_{A}\,\epsilon_{A}^{i}(\vec{n})\,\epsilon_{A}^{j}(\vec{n})\,n^{i}_{a}\,n^{j}_{b}=\cfrac{2}{3}\,\cos\theta_{ab}.$

## References

  * [1] G. Agazie et al. [NANOGrav], “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,” Astrophys. J. Lett. 951, no.1, L8 (2023) [arXiv:2306.16213 [astro-ph.HE]].
  * [2] A. Afzal et al. [NANOGrav], “The NANOGrav 15 yr Data Set: Search for Signals from New Physics,” Astrophys. J. Lett. 951, no.1, L11 (2023) [arXiv:2306.16219 [astro-ph.HE]].
  * [3] G. Agazie et al. [NANOGrav], “The NANOGrav 15 yr Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational-wave Background,” Astrophys. J. Lett. 952, no.2, L37 (2023) [arXiv:2306.16220 [astro-ph.HE]].
  * [4] R. W. Hellings and G. S. Downs, Upper limits on the istotropic gravitational radiation background from pulsar timing analysis, Astrophys. J. , 265, pp. L39-L42 (1983).
  * [5] J. Antoniadis, P. Arumugam, S. Arumugam, S. Babak, M. Bagchi, A. S. B. Nielsen, C. G. Bassa, A. Bathula, A. Berthereau and M. Bonetti, et al. “The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,” [arXiv:2306.16214 [astro-ph.HE]].
  * [6] P. Tarafdar, N. K., P. Rana, J. Singha, M. A. Krishnakumar, B. C. Joshi, A. K. Paladi, N. Kolhe, N. D. Batra and N. Agarwal, et al. “The Indian Pulsar Timing Array: First data release,” Publ. Astron. Soc. Austral. 39, e053 (2022) [arXiv:2206.09289 [astro-ph.IM]].
  * [7] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane and J. Askew, et al. “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951, no.1, L6 (2023) [arXiv:2306.16215 [astro-ph.HE]].
  * [8] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Res. Astron. Astrophys. 23, no.7, 075024 (2023) [arXiv:2306.16216 [astro-ph.HE]].
  * [9] K. J. Lee, F. A. Jenet, and R. H. Price, Pulsar Timing as a Probe of Non-Einsteinian Polarizations of Gravitational Waves, Astrophys. J. , 685, pp. 1304-1319, (2008).
  * [10] N. J. Cornish, L. O’Beirne, S. R. Taylor and N. Yunes, “Constraining alternative theories of gravity using pulsar timing arrays,” Phys. Rev. Lett. 120, no.18, 181101 (2018) [arXiv:1712.07132 [gr-qc]].
  * [11] Z. C. Chen, Y. M. Wu, Y. C. Bi and Q. G. Huang, “Search for Non-Tensorial Gravitational-Wave Backgrounds in the NANOGrav 15-Year Data Set,” [arXiv:2310.11238 [astro-ph.CO]].
  * [12] G. Agazie et al. [NANOGrav], “The NANOGrav 15-year data set: Search for Transverse Polarization Modes in the Gravitational-Wave Background,” [arXiv:2310.12138 [gr-qc]].
  * [13] K. Nomura, A. Ito and J. Soda, “Pulsar timing residual induced by ultralight vector dark matter,” Eur. Phys. J. C 80, no.5, 419 (2020) [arXiv:1912.10210 [gr-qc]].
  * [14] X. Xue et al. [PPTA], “High-precision search for dark photon dark matter with the Parkes Pulsar Timing Array,” Phys. Rev. Res. 4, no.1, L012022 (2022) [arXiv:2112.07687 [hep-ph]].
  * [15] R. E. Marshak and R. N. Mohapatra, Phys. Lett. B 91, 222-224 (1980) doi:10.1016/0370-2693(80)90436-0
  * [16] S. Khalil, J. Phys. G 35, 055001 (2008) doi:10.1088/0954-3899/35/5/055001 [arXiv:hep-ph/0611205 [hep-ph]].
  * [17] L. Basso, A. Belyaev, S. Moretti and C. H. Shepherd-Themistocleous, Phys. Rev. D 80, 055030 (2009) doi:10.1103/PhysRevD.80.055030 [arXiv:0812.4313 [hep-ph]].
  * [18] S. Iso, N. Okada and Y. Orikasa, Phys. Lett. B 676, 81-87 (2009) doi:10.1016/j.physletb.2009.04.046 [arXiv:0902.4050 [hep-ph]].
  * [19] N. Okada and O. Seto, Phys. Rev. D 82, 023507 (2010) doi:10.1103/PhysRevD.82.023507 [arXiv:1002.2525 [hep-ph]].
  * [20] L. Basso, S. Moretti and G. M. Pruna, Phys. Rev. D 82, 055018 (2010) doi:10.1103/PhysRevD.82.055018 [arXiv:1004.3039 [hep-ph]].
  * [21] L. Basso, S. Moretti and G. M. Pruna, Phys. Rev. D 83, 055014 (2011) doi:10.1103/PhysRevD.83.055014 [arXiv:1011.2612 [hep-ph]].
  * [22] S. Kanemura, O. Seto and T. Shimomura, Phys. Rev. D 84, 016004 (2011) doi:10.1103/PhysRevD.84.016004 [arXiv:1101.5713 [hep-ph]].
  * [23] S. Sun, X. Y. Yang and Y. L. Zhang, “Pulsar timing residual induced by wideband ultralight dark matter with spin 0,1,2,” Phys. Rev. D 106, no.6, 066006 (2022) [arXiv:2112.15593 [astro-ph.CO]].
  * [24] C. Unal, F. R. Urban and E. D. Kovetz, “Probing ultralight scalar, vector and tensor dark matter with pulsar timing arrays,” [arXiv:2209.02741 [astro-ph.CO]].
  * [25] Y. M. Wu et al. [PPTA], “Constraining ultralight vector dark matter with the Parkes Pulsar Timing Array second data release,” Phys. Rev. D 106, no.8, L081101 (2022) [arXiv:2210.03880 [astro-ph.CO]].
  * [26] Z. Q. Xia, T. P. Tang, X. Huang, Q. Yuan and Y. Z. Fan, “Constraining ultralight dark matter using the Fermi-LAT pulsar timing array,” Phys. Rev. D 107, no.12, L121302 (2023) [arXiv:2303.17545 [astro-ph.HE]].
  * [27] H. Omiya, K. Nomura and J. Soda, “Hellings-Downs curve deformed by ultralight vector dark matter,” [arXiv:2307.12624 [astro-ph.CO]].
  * [28] J. M. Pearson, N. Chamel, A. Y. Potekhin, A. F. Fantina, C. Ducoin, A. K. Dutta and S. Goriely, “Unified equations of state for cold non-accreting neutron stars with Brussels–Montreal functionals – I. Role of symmetry energy,” Mon. Not. Roy. Astron. Soc. 481, no.3, 2994-3026 (2018) [erratum: Mon. Not. Roy. Astron. Soc. 486, no.1, 768 (2019)] [arXiv:1903.04981 [astro-ph.HE]].
  * [29] A. Y. Potekhin, A. F. Fantina, N. Chamel, J. M. Pearson and S. Goriely, “Analytical representations of unified equations of state for neutron-star matter,” Astron. Astrophys. 560, A48 (2013) [arXiv:1310.0049 [astro-ph.SR]].
  * [30] R. Garani and J. Heeck, “Dark matter interactions with muons in neutron stars,” Phys. Rev. D 100, no.3, 035039 (2019) [arXiv:1906.10145 [hep-ph]].
  * [31] R. Foot, “New Physics From Electric Charge Quantization?,” Mod. Phys. Lett. A 6, 527-530 (1991)
  * [32] X. G. He, G. C. Joshi, H. Lew and R. R. Volkas, “Simplest Z-prime model,” Phys. Rev. D 44, 2118-2132 (1991)
  * [33] R. Foot, X. G. He, H. Lew and R. R. Volkas, “Model for a light Z-prime boson,” Phys. Rev. D 50, 4571-4580 (1994) [arXiv:hep-ph/9401250 [hep-ph]].
  * [34] J. Heeck and W. Rodejohann, Phys. Rev. D 84, 075007 (2011) doi:10.1103/PhysRevD.84.075007 [arXiv:1107.5238 [hep-ph]].
  * [35] J. A. Grifols and E. Masso, “Neutrino oscillations in the sun probe long range leptonic forces,” Phys. Lett. B 579, 123-126 (2004) [arXiv:hep-ph/0311141 [hep-ph]].
  * [36] A. S. Joshipura and S. Mohanty, “Constraints on flavor dependent long range forces from atmospheric neutrino observations at super-Kamiokande,” Phys. Lett. B 584, 103-108 (2004) [arXiv:hep-ph/0310210 [hep-ph]].
  * [37] J. Heeck and W. Rodejohann, “Gauged $L_{\mu}-L_{\tau}$ and different Muon Neutrino and Anti-Neutrino Oscillations: MINOS and beyond,” J. Phys. G 38, 085005 (2011) [arXiv:1007.2655 [hep-ph]].
  * [38] Y. Farzan and J. Heeck, “Neutrinophilic nonstandard interactions,” Phys. Rev. D 94, no.5, 053010 (2016) [arXiv:1607.07616 [hep-ph]].
  * [39] A. Y. Smirnov and X. J. Xu, “Wolfenstein potentials for neutrinos induced by ultra-light mediators,” JHEP 12, 046 (2019) [arXiv:1909.07505 [hep-ph]].
  * [40] J. A. Dror, “Discovering leptonic forces using nonconserved currents,” Phys. Rev. D 101, no.9, 095013 (2020) [arXiv:2004.04750 [hep-ph]].
  * [41] M. Singh, M. Bustamante and S. K. Agarwalla, “Flavor-dependent long-range neutrino interactions in DUNE & T2HK: alone they constrain, together they discover,” JHEP 08, 101 (2023) [arXiv:2305.05184 [hep-ph]].
  * [42] T. Kumar Poddar, S. Mohanty and S. Jana, “Vector gauge boson radiation from compact binary systems in a gauged $L_{\mu}-L_{\tau}$ scenario,” Phys. Rev. D 100, no.12, 123023 (2019) [arXiv:1908.09732 [hep-ph]].
  * [43] J. Kopp, R. Laha, T. Opferkuch and W. Shepherd, “Cuckoo’s eggs in neutron stars: can LIGO hear chirps from the dark sector?,” JHEP 11, 096 (2018) [arXiv:1807.02527 [hep-ph]].
  * [44] E. S. Phinney, “A Practical theorem on gravitational wave backgrounds,” [arXiv:astro-ph/0108028 [astro-ph]].
  * [45] W. Hu, R. Barkana and A. Gruzinov, “Cold and fuzzy dark matter,” Phys. Rev. Lett. 85, 1158-1161 (2000) [arXiv:astro-ph/0003365 [astro-ph]].
  * [46] L. Hui, J. P. Ostriker, S. Tremaine and E. Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D 95, no.4, 043541 (2017) [arXiv:1610.08297 [astro-ph.CO]].
  * [47] T. Kumar Poddar, S. Mohanty and S. Jana, “Constraints on long range force from perihelion precession of planets in a gauged $L_{e}-L_{\mu,\tau}$ scenario,” Eur. Phys. J. C 81, no.4, 286 (2021) [arXiv:2002.02935 [hep-ph]].
  * [48] D. Brzeminski, S. Das, A. Hook and C. Ristow, “Constraining Vector Dark Matter with Neutrino experiments,” [arXiv:2212.05073 [hep-ph]].