# Searching for Time-Dependent Axion Dark Matter Signals in Pulsars

R. A. Battye<EMAIL_ADDRESS>Jodrell Bank Centre for
Astrophysics, School of Natural Sciences, Department of Physics and Astronomy,
University of Manchester, Manchester, M13 9PL, U.K.  M. J. Keith
<EMAIL_ADDRESS>Jodrell Bank Centre for Astrophysics, School of
Natural Sciences, Department of Physics and Astronomy, University of
Manchester, Manchester, M13 9PL, U.K.  J. I. McDonald
<EMAIL_ADDRESS>Centre for Cosmology, Particle Physics and
Phenomenology, Université Catholique de Louvain, Chemin du cyclotron 2,
Louvain-la-Neuve B-1348, Belgium  S. Srinivasan
<EMAIL_ADDRESS>Jodrell Bank Centre for Astrophysics,
School of Natural Sciences, Department of Physics and Astronomy, University of
Manchester, Manchester, M13 9PL, U.K.  B. W. Stappers
<EMAIL_ADDRESS>Jodrell Bank Centre for Astrophysics, School of
Natural Sciences, Department of Physics and Astronomy, University of
Manchester, Manchester, M13 9PL, U.K.  P. Weltevrede
<EMAIL_ADDRESS>Jodrell Bank Centre for Astrophysics,
School of Natural Sciences, Department of Physics and Astronomy, University of
Manchester, Manchester, M13 9PL, U.K.

###### Abstract

Axion dark matter can be converted into photons in the magnetospheres of
neutron stars leading to a spectral line centred on the Compton wavelength of
the axion. Due to the rotation of the star and the plasma effects in the
magnetosphere the signal is predicted to be periodic with significant time
variation that persists across phase but is narrow in frequency - a unique
smoking gun for axion dark matter. As a proof of principle and to develop the
methodology, we carry out the first time domain search of the signal using
data from PSR J2144$-$3933 taken as part of the MeerTIME project on MeerKAT
telescope. We search for specific signal templates using a matched filter
technique and discuss when a time-domain analysis (as is typically the case in
pulsar observations) gives greater sensitivity to the axion-coupling to
photons in comparison to a simple time-averaged total flux study. We do not
find any candidate signals and, hence, impose an upper limit on the axion-to-
photon coupling of $g_{a\gamma\gamma}<5.5\times 10^{-11}(D/0.165\,{\rm
pc})\,{\rm GeV}^{-1}$ where $D$ is the pulsar distance, over the mass range
$m_{\rm a}=3.9-4.7\,\mu{\rm eV}$ using this data. This limit relies on PSR
J2144$-$3933 not being an extremely aligned rotator, as strongly supported by
simple arguments based on the observed pulse profile width. We discuss the
possibilities of improving this limit using future observations with MeerKAT
and also SKA1-mid and the possibility of using other objects. Finally, to
evade modelling uncertainties in axion radio signals, we also carry out a
generic “any periodic-signal search” in the data, finding no evidence for an
axion signal.

Axions; Dark matter; Neutron stars

###### pacs:

95.35.+d; 14.80.Mz; 97.60.Jd

††preprint: IRMP-CP3-23-14

##  I Introduction

The search for dark matter in the form of axions [1, 2, 3, 4, 5, 6, 7, 8, 9,
10] is continuing to gather momentum at a dramatic pace. Of particular
interest is the mass range $m_{\rm a}\approx 0.1-1000\,\mu{\rm eV}$ where
plausible scenarios [11, 12, 13, 14, 15, 16, 17, 18] have been proposed to
realise the Cold Dark Matter abundance $\Omega_{\rm c}h^{2}\approx 0.12$ that
is found to be compatible with cosmological observations, for example, those
of the CMB [19]. There are a number of haloscope experiments [20] which have
placed constraints on the axion coupling to photons. These include early
cavity experiments [21, 22, 23] and their modern incarnations including ADMX
[24, 25, 24, 26, 27, 28] and its various upgrades and pathfinders [29, 30,
31]. This has spawned a plethora of active and proposed experiments including
CAPP [32, 33, 34, 35], HAYSTAC [36, 37, 38] QUAX [39, 40], ORGAN [41, 42],
CAST-RADES [43], TASEH [44] and GrAHal [45] aiming to detect axions using
similar techniques. Complementing this, there are also a number of proposed
experiments which go beyond the cavity paradigm, allowing laboratory searches
to achieve broad frequency coverage across previously challenging frequency
ranges. These include plasma haloscope designs like ALPHA [46], broadband
reflectors envisaged in the BREAD [47] collaboration, and dielectric
haloscopes such as MADMAX [48]. There are also novel designs which aim to
detect magnetic fields induced by axion sources such as DMRadio [49], which
will operate at lower frequencies ($m_{\rm a}\lesssim\mu{\rm eV}$) and those
seeking to using novel materials at higher frequencies [50, 51, 52]

The largest magnetic fields currently used in laboratory searches for axion
dark matter typically do not exceed $\sim 10^{5}{\rm G}\,(10\,{\rm T})$, a
limiting factor in such searches. By contrast, astrophysical magnetic fields
in neutron stars can be as high as $\sim 10^{15}{\rm G}\,(10^{11}\,{\rm T})$,
making them excellent targets for indirect searches of axion dark matter [53,
54, 55]. In addition, neutron stars are surrounded by a magnetosphere whose
varying plasma frequency matches the axion mass across a broad range of
masses. This degeneracy leads to a dramatic resonant enhancement of the signal
emanating from regions with $m_{a}\simeq\omega_{\rm p}$, where $\omega_{\rm
p}$ is the plasma frequency. As a result, neutron stars can act as broadband
axion dark matter detectors.

Based on a simple but representative model for a neutron star magnetosphere
and the density of axions around the star, Refs. [56, 55, 54] predicted
signals that could be easily detected using current and future telescopes
operating in the radio-mm waveband, which corresponds to the Compton
wavelength of the dark matter axion scenarios referred to above. Spurred on by
this, great progress has recently been made in characterising the signal
properties using sophisticated ray-tracing methods [57, 58, 59] which are
capable of computing the line width induced from plasma effects and the
precise time-variation and angular dependence of the signal. Early attempts
have also been made to address axion-photon mixing in 3D [60, 61], though this
remains an ongoing area of research.

Various searches have been carried out to detect radio signals produced by
axion dark matter converting into photons in the magnetospheres of neutron
stars using the Goldreich-Julian (GJ) model [62] for the magnetosphere and
estimates of the local density of dark matter, extrapolated to the location of
the star in question. These searches have either looked for a background
excess near the Galactic centre from populations of neutron stars [63, 64], or
have focused on single objects such as the Galactic centre magnetar [65, 66,
67] or isolated neutron stars [68, 69], and have established bounds on the
axion-to-photon coupling, $g_{\rm a\gamma\gamma}$, which are better than the
bounds from the axion helioscope CAST [70].

Now that the first wave of searches has been carried out, a natural question
to ask is what improved observational strategies are available to increase our
sensitivity to the axion photon coupling and boost our chances of detecting
dark matter axions from neutron stars. In this vein, one might also ask how
our newly attained understanding of precise signal properties (including time
and frequency information) might be leveraged to increase the power of such
searches. To date, all the searches for axions using neutron stars have
focused on looking for a spectral line in the frequency domain. The goal of
the present work is twofold: (i) to establish a framework of time-domain
searches for axion dark matter signals in radio data and (ii) to demonstrate
this technique by carrying out time-domain observations of pulsars. Our goal
is to understand under what circumstances augmenting these searches to include
time-domain information of the signal can improve sensitivity to
$g_{a\gamma\gamma}$ and by injecting signal templates into radio data,
demonstrate a practical route to obtaining limits on $g_{a\gamma\gamma}$ which
out-perform a simple line-search in the frequency domain.

The structure of the paper is organised as follows. In sec II we review the
mechanism for photon production from axion dark matter, describe our ray-
tracing procedure for modelling the radio signal and discuss the time-
dependence of the expected signal. In sec. III we describe a procedure to
search for time-dependent signals in data using a matched filter, and use this
to outline what types of periodic signals lead to a strong gain from time-
domain information. In sec. IV we apply our pipeline to MeerKAT observations
of PSR J2144$-$3933 to search for axion dark matter. Our null result is used
to place limits on the axion coupling to photons. In sections V and VI we
explore possible future targets and perform a generalized search for periodic
signals. In section VII we offer our conclusions.

##  II Modelling the signal due to axions

The conversion between axions and photons in strong magnetic fields was laid
out in the classic reference [71]. It was pointed out in [53] that this mixing
could convert dark matter axions into radio photons111See [72] for a similar
mechanism with dark photons. in the strongly magnetised plasmas which surround
neutron stars. In recent years, as axions have moved to the forefront as dark
matter candidates, these ideas have been pursued with renewed vigour [55, 73]
and this programme has lead to a variety of observations [74, 75, 67, 76, 64]
searching for radio lines from dark matter axions. These observational efforts
have been accompanied by a more concerted effort to improve the modelling of
the signal itself [77, 59, 58, 61, 78] which consists primarily of developing
ray-tracing packages to precisely track the photons from their point of
emission to the observer, thereby allowing one to derive signal templates
which could be detected by a radio telescope. We now briefly review the basic
features of the production mechanism and ray-tracing routine. More details can
be found in [77, 59, 58, 79].

###  II.1 Ray-tracing photons in magnetised plasmas

Resonant conversion between axions and photons occurs at points where
$k^{\mu}_{\gamma}=k_{a}^{\mu}$ where $k^{\mu}_{\gamma}$ and $k_{a}^{\mu}$ are
the photon and axion 4-momentum, respectively. An attempt was made to
understand the conversion probability for axions to photons $p_{a\gamma}$ in
[61] leading to

$p_{a\gamma\gamma}=\frac{\pi}{2}\frac{g_{a\gamma\gamma}^{2}\sin^{2}\theta_{B}\left|\textbf{B}\right|^{2}}{\left|\textbf{k}_{\gamma}\right|\left|\omega_{p}^{\prime}\right|}\cdot\frac{m_{a}^{5}}{\left(\left|\textbf{k}_{\gamma}\right|^{2}+m_{a}^{2}\sin^{2}\theta_{B}\right)^{2}},$
(1)

which attempts to incorporate 3D effects into the conversion probability. This
characterises the ratio of the energy density between an axion wavepacket and
a photon wavepacket, the latter being subsequently transported out of the
magnetosphere along geodesics determined by the photon dispersion relation in
the strongly magnetised plasma. Note in the present work we include
simultaneously the effects of gravity (by incorporating the curved spacetime
metric of the neutron star) and strongly magnetised fields. This results in a
covariant dispersion relation for photons in a magnetised plasma which have a
dispersion relation [80, 81]

$\displaystyle{\cal D}(k)=g^{\mu\nu}k_{\mu}k_{\nu}$
$\displaystyle-(\omega^{2}-k_{\parallel}^{2})\sum_{s}\frac{4\pi
q_{s}^{2}n_{s}}{\gamma_{s}^{2}[\mu_{s}(\omega-
k_{\parallel})^{2}-c_{s}^{2}(\omega v_{s}-k_{\parallel})^{2}]}.$ (2)

Here, the sum $s$ is over different charge carrier species, $\gamma_{s}$ is a
generalised Lorentz factor, and $k_{\parallel}$ gives the 4-momentum projected
onto the magnetic field and $v_{s}$ corresponds to the velocity of charge
carriers and $\mu_{s}$ the energy per particle. The number density of each
species $s$ is given by $n_{s}$ and the charge by $q_{s}$. Full definitions
and detailed explanations of the various terms can be found in [80]. We will
take the non-relativistic limit, setting $v_{s}=0$, $\gamma=1$, $c_{s}=0$ and
$\mu_{s}=m_{s}$. We also consider a purely electron-positron plasma so that
$q_{s}=e$. Setting $D(k)=0$ then gives the dispersion relations for photons in
a non-relativistic plasma. The equations of motion for the photon rays are
then given by Hamilton’s equations

$\frac{dx^{\mu}}{d\lambda}=\frac{\partial{\cal D}}{\partial
k_{\mu}},\qquad\frac{dk_{\mu}}{d\lambda}=-\frac{\partial{\cal D}}{\partial
x^{\mu}}.$ (3)

To compute the power, we back-trace from the observer to the point of
emission, following the equations of motion (3). This is analogous to the
procedure used in [57, 58].

In addition, we also now include multiple axion-photon conversions arising
from multiple reflections off the critical surface, as happens within
“throats” of the plasma distribution around the neutron star. These throats
are partially enclosed regions near the charge separation gap off whose walls
the photon can be multiply-reflected due to plasma gradients. This can enhance
the power of the signal relative to not including such effects as was done in
[57, 58]. An extensive analysis of ray-tracing techniques which combines the
physical effects considered across [58] and [59] is currently underway and
will appear in a companion paper [79], where the full details of our scheme
will be presented. This will include a systematic study of anisotropic
plasmas. We do not consider the effects of so-called “de-phasing” conjectured
in [59] which awaits a more robust physical description to see if the effect
persists under more mathematically rigorous formulation. We do not need to
consider non-linear effects arising from very large conversion probabilities
where photons may convert back into axions. This is safe for PSR J2144$-$3933
on which we performed our observations, which has sufficiently low magnetic
fields that the conversion probability remains small.

###  II.2 Signal templates

Having outlined the basic details of our ray-tracing scheme. This can now be
used to begin deriving signal templates. In particular, these simulations
allow one to model the radio signal as a function of pulsar and axion input
parameters. In particular, one can compute the profile of the signal in
frequency and time. The frequency dependence of the profile is determined by
the mass of the axion - which sets the central frequency of the radio line.
The width is set by a combination of the velocity dispersion of dark matter
and by line-broadening induced by the time-dependent nature of the plasma,
which modifies photon frequencies as they move through the magnetosphere. The
time variation of the signal arises from the fact that the plasma surrounding
the neutron star is not axisymmetric. In the present approximation we assume
an electron-positron plasma which co-rotates with the star with regions of
positive and negative charge separated according to the Goldreich-Julian
density [62]

$\displaystyle n_{\rm
GJ}=\frac{B_{0}\Omega}{2\,e}\left(\frac{R}{r}\right)^{3}\Big{[}$
$\displaystyle\cos\alpha+3\cos\alpha\cos(2\theta)$
$\displaystyle+3\sin\alpha\cos(\phi-\Omega t)\sin 2\theta\Big{]}\,,$ (4)

with the plasma frequency given by $\omega_{\rm p}=\sqrt{4\pi|n_{\rm
GJ}|/m_{e}}$. Here, $\Omega$ is the frequency of the pulsar, $\alpha$ is the
angle between the magnetic axis of the co-rotating dipole and the rotation
axis of the star, $R$ is the stellar radius, and $B_{0}$ is the magnetic field
strength on the surface at the magnetic poles. The polar coordinate $\theta$
and azimuthal angle $\phi$ are defined with respect to the rotational axis of
the star. It is obvious that whenever $\alpha\neq 0$, the plasma is time-
dependent with respect to a non-rotating observer. This results in time-
dependent radio signals, as illustrated in Fig. 1. For a given pulsar, the
remaining input parameters to determine the axion dark matter radio signal are
then the distance to the pulsar $D$, and the dark matter density $\rho_{\rm
DM}$ at the position of the pulsar.

In our analysis, we will take $P$ and $B_{0}$ to be their quoted measured
values. In principle there are some extra uncertainties which would need to be
taken into account. Pulsar periods $P$ are one of the best-measured quantities
in astronomy, while the magnetic field strength of the pulsar at the pole is
inferred from measurements of $P$ and $\dot{P}$ the spin-down rate of the
pulsar combined with model-dependent parameters including the moment of
inertia of the pulsar and its radius [82]. This calculation is standard but
assumes that the energy released by the pulsar in the form of radio emission
comes from the loss of rotational energy calculated from the spin-down rate.
The comparatively large values of $P\dot{P}$ observed for magnetars form the
basis for their large inferred values of $B_{0}$, but it is also known that
the magnetars are known to emit large X-ray fluxes whose luminosity cannot be
explained by spin-down alone. The emission mechanism for some magnetars may be
powered magnetically in contrast to the rotation powered radio pulsars [83,
84]. In particular, magnetars can have a significant toroidal component to
their field whose magnitude one can constrain by imposing that the luminosity
from magnetic dissipation exceeds that from the traditional spin-down. Based
on observations of the X-ray flares/bursts from magnetars, the typical ratio
of the toroidal to poloidal component of the field is of the order $\sim 10$
[85, 86]. While it is possible to calculate this enhancement from running
magneto-hydrodynamic (MHD) simulations, the impact on the final constraints is
suppressed by $\sqrt{S}\propto B^{0.4}$ (see sec. V) and s such, we leave a
more sophisticated modelling of the magnetosphere for future work.

The distance to the pulsar, $D$ can be inferred from parallax measurements or
from the dispersion of the pulse as a function of frequency, since the photons
emitted in the main beam of the pulsar traverse through the galactic electron
density along the line-of-sight. Given a model for the galactic electron
density, one can estimate the distance to a pulsar. If we use the pyGEDM code
and the galactic coordinates as measured in the catalogue we obtain $D\approx
0.29\,{\rm kpc}$ [87, 88, 89, 90, 91, 92]. This illustrates that this
computation relies on the underlying model for the galactic electron density
distribution. A parallax measurement for this pulsar yields
$D=0.165^{+0.017}_{-0.014}\,{\rm kpc}$ [93]. This is a more direct measurement
and we will use it as our fiducial value when establishing limits, but we note
that the limit is $\propto D$. Galactic dark matter profiles allow one to
predict the dark matter density at the position of the pulsar, but these
models become highly uncertain as one gets closer to the galactic centre,
where some models predict a spike in the density, while others predict a more
cored profile.

Based on these arguments, we conclude that the completely unknown quantities
which parameterise the signal templates are $(\alpha,\theta)$. We, therefore,
generate a simulated database222Formally we generate templates for discrete
$(\alpha_{i},\theta_{i})$ and numerically interpolate to generate a template
for a continuous range of $\alpha$ and $\theta$. See appendix A for a
description of this procedure and an illustration of its accuracy. of periodic
flux profiles as a function of $(\alpha,\theta)$. Some of these profiles are
displayed in Fig. 1 which indicates the $(\alpha,\theta)$ dependence of the
time-variability of the signal.

In the next section, we describe how to harness the information and the larger
time-variability of the signal to improve the prospects of detecting axion
dark matter.

Figure 1: Time-dependence of radio templates from axion dark matter as a
function of pulse phase $\phi$ for different values of $\theta$ with fixed
value of $\alpha=20^{\circ}$ with $g_{a\gamma\gamma}=10^{-10}{\rm GeV}^{-1}$.
As expected there is little time variation for $\theta=10^{\circ}$, but it can
be substantial for intermediate angles. When $\theta=90^{\circ}$ the “throats”
of the GJ model never cross the line-of-sight so, although there is some
emission, it is relatively weak compared to cases where they are visible to
the observer.

##  III Searching for time-dependent signals

In order to search for time dependent signals we will employ a matched-filter
template-fitting approach similar to that used in Gravitational Wave Astronomy
to detect the waveforms of the late stages of binary black hole inspirals [94,
95]. In that case, once a detection of gravitational waves was made, this
allows estimates of physical parameters such as the black hole masses.
Formally, if an axion were detected, one could use axion radio signals to fit
model parameters of the pulsar magnetosphere. However, our ambition at this
stage is much more conservative: we will use the signal-to-noise estimate from
the matched filter, ${\hat{q}}$, defined below, to quantify the likelihood of
detection. In this sense, ${\hat{q}}$ acts as a statistical test for whether
the data is distinguishable from noise. Values of $\hat{q}$ above a threshold
then constitute a detection. Conversely, for values below this, by injecting
would-be signals into the data, we obtain the expected value $q_{\rm
exp}=<\hat{q}>$ (see eq. (15)) from an axion signal. By comparing this to the
measured value $\hat{q}$, we can exclude regions of axion parameter space.

This procedure provides a means to derive limits on $g_{a\gamma\gamma}$ as a
function of $m_{\rm a}$, and importantly allows us to take into account that
for a fixed value of $m_{\rm a}$ there are a wide range of templates for the
expected signal due to the parameters of the particular neutron star system
under consideration. These are the period of the pulsar, $P$, its surface
magnetic field flux density, $B_{0}$, the radius of the neutron star, $R$, the
angle $\alpha$ between the magnetic axis and spin axis of the star, and the
angle $\theta$ between the line of sight and the spin axis. As in previous
attempts to derive constraints on $g_{a\gamma\gamma}$ using neutron stars [68,
65, 66, 67, 64] for simplicity, as a demonstration of the filter, we do not,
for instance consider uncertainties in $B_{0}$ or $R$ (the period, $P$ is of
course measured with tremendous accuracy). We leave a computationally
intensive parameter scan for future work, but this would be a straightforward
extension of the existing framework.

Instead, our main focus here is on the sensitivity of the time-dependence of
the signal to pulsar parameters, which is especially sensitive to the value of
$\alpha$ and $\theta$. For each value of $m_{\rm a}$, we therefore obtain a
constraint on the value of $g_{a\gamma\gamma}$ for every pair
($\alpha$,$\theta$). We can then exclude certain ranges of $\alpha$ and
$\theta$ with further modelling and observations of the pulsar signal, notably
the pulse width. We then use the value of $(\alpha,\theta)$ from this
remaining subset which gives the most conservative constraints on
$g_{a\gamma\gamma}$.

###  III.1 Derivation of matched filter and mathematical properties

Matched filters are a standard technique in signal processing and they are
often used in astronomy to search for signals with a known, or parameterizable
profile. This can be done in the spatial, frequency or time domains. In this
section, we will derive the standard matched filter before discussing some of
its properties. In section III.2 we use it to quantify the pros and cons of
time domain observations and in section III.3 we will show how it can be used
to recover an injected signal in simulated radio data.

There are a number of ways of formulating the matched filter. Here, we will
use a discrete matched filter in which the data is represented by a vector of
finite length. This can be generalised to continuous functions in which the
vector inner products become convolution integrals over functions (see [95]).
The discrete formulation has the advantage of simplifying notation.

Our starting point is the so-called “data vector”, ${\bf d}$. We will assume
the data is the sum of a signal $\textbf{S}=S_{0}\textbf{F}(\textbf{p})$ and
some additive noise $\hat{\textbf{n}}$

$\textbf{d}=S_{0}\textbf{F}(\textbf{p})+\hat{\textbf{n}}\,.$ (5)

Here, we have decomposed the signal according to

$|{\bf F}|\equiv\sqrt{\textbf{F}\cdot\textbf{F}}=1$ (6)

so that $S_{0}$, gives the root mean squared flux density of the signal

$S_{0}=\sqrt{\textbf{S}\cdot\textbf{S}}\,.$ (7)

The signal is further characterised by a “parameter vector” ${\bf p}$. In
section II we will compute a number of templates ${\bf F}({\bf p})$ for the
signal where ${\bf p}=(m_{\rm a},\alpha,\theta)$. The noise vector
$\hat{\textbf{n}}$ is assumed to be Gaussian with $\langle{\bf n}\rangle=0$
and $\langle{\bf n}{\bf n}^{T}\rangle=C$ where $\langle..\rangle$ denotes an
ensemble average of noise realizations and $C_{ij}$ is the covariance matrix
with $i,j=1,..,n_{\rm d}$ where $n_{d}$ is the total number of data points.

Assuming a Gaussian likelihood, $-2\log{\cal L}=\chi^{2}$, one can calculate
the maximum likelihood estimate ${\hat{S}}_{0}$ by minimizing

$\chi^{2}=({\bf d}-S_{0}{\bf F})^{T}C^{-1}({\bf d}-S_{0}{\bf F})\,.$ (8)

In the above equation, we assume the data vector d contains the true signal
$\textbf{S}=\textbf{S}_{\rm true}$ so that minimising $\chi^{2}$ above can be
thought of as minimising a generalised least-squares difference (relative to
the noise in each channel - hence the factor $C^{-1}$) between the true signal
and possible templates $S_{0}\textbf{F}$. Viewing $\chi^{2}$ as an unknown
function of $S_{0}$, we can find the minimising value, $\hat{S}_{0}$ given by

${\hat{S}}_{0}=\frac{\textbf{F}^{T}C^{-1}\textbf{d}}{\textbf{F}^{T}C^{-1}\textbf{F}}\,.$
(9)

One can also deduce the “matched filter noise”

$\sigma_{\rm MF}=\left(\textbf{F}^{T}C^{-1}\textbf{F}\right)^{-1/2}$ (10)

and the signal-to-noise estimate ${\hat{q}}$ may then be written as

${\hat{q}}=\frac{{\hat{S}}_{0}}{\sigma_{\rm
MF}}=\frac{\textbf{F}^{T}C^{-1}\textbf{d}}{\left(\textbf{F}^{T}C^{-1}\textbf{F}\right)^{1/2}}\,.$
(11)

It is important to understand the difference between the noise in the data,
characterised by $C$, and the “matched filter noise”, $\sigma_{\rm MF}$. They
are related, but $\sigma_{\rm MF}$ also depends on the filter. We return to
this issue in the next section.

In order to understand properties of the matched filter we will assume
diagonal covariance matrix

$C_{ij}=\sigma_{N}^{2}\delta_{ij}$ (12)

where $\delta_{ij}$ is the Kronecker delta and $\sigma_{N}$ is the noise in
each channel - all of what said here can be adapted to the case of a general
covariance matrix, but it is less simple to see.

In general, the data has a number of dimensions, for example, space, time
and/or frequency, then d has $n_{\rm d}=n_{1}\times...\times n_{k}$ entries
where $n_{i}$ for $i=1,...,k$ are the number of points in each of the
dimensions. In our case we will search in the frequency and time directions so
the number of entries in the data vector will be $n_{\rm d}=n_{\rm f}\times
n_{\rm t}$ where $n_{\rm f}$ is the number of frequency channels and $n_{\rm
t}$ is the number of time samples. In that case, the data vector could be
written as

$\textbf{d}=\Big{(}d^{\omega_{1}}_{t_{1}},\dots,d^{\omega_{1}}_{t_{n_{t}}},\dots,d^{\omega_{n_{\rm
f}}}_{t_{1}},\dots,d^{\omega_{n_{\rm f}}}_{t_{n_{t}}}\Big{)},$ (13)

where $d^{\omega_{i}}_{t_{j}}$ labels the data vector in the $i{\rm th}$
frequency bin and $j{\rm th}$ time-channel. The covariance matrix of the form
(12) is then nothing more than the statement that the noise between all
possible pairs of time and frequency bins is totally uncorrelated.

Returning to our main discussion, it follows that for a covariance matrix of
the form (12), we have

${\hat{q}}=\frac{{\bf F}\cdot{\bf d}}{\sigma_{N}}\,,$ (14)

that is, the dot product of the filter F with the data vector. When F matches
that in the true signal, we have $\langle{\hat{q}}\rangle=S_{0}/\sigma_{N}$
(from eq. (5)).

Now assume that ${\bf p}_{\rm true}$ has entries which are the true
parameters. The ensemble average of ${\hat{q}}$ for a filter with arbitrary
parameter, ${\bf p}$ is

$\langle{\hat{q}}\rangle={\bf F}({\bf p})\cdot{\bf F}({\bf p}_{\rm
true}){S_{0}\over\sigma_{N}}\,.$ (15)

Next, we note that as a trivial consequence of the Cauchy-Schwarz inequality,
we have ${\bf F}({\bf p})\cdot{\bf F}({\bf p}_{\rm true})\leq\left|{\bf
F}({\bf p})\right|\left|{\bf F}({\bf p}_{\rm true})\right|$ with equality if
and only if ${\bf F}({\bf p})={\bf F}({\bf p}_{\rm true})$. Assuming that the
template is non-degenerate with respect to the values of p, this occurs only
for ${\bf p}={\bf p}_{\rm true}$ for which ${\hat{q}}$ is then maximal. Thus
$\hat{q}$ acts as a likelihood test for the values of p.

In what follows, since the line width of the axion is typically less than the
width of our frequency channels, the vector (13) is sparse for a given value
of $m_{a}$, with non-vanishing entries in only one frequency channel where
$\omega=m_{a}$. This means filters with different values of $m_{\rm a}$ are
orthogonal. More generally, with higher frequency resolution, we would expect
to be able to probe both the time and frequency structure of the signal. By
contrast, filters with different values of $\alpha$ and $\theta$ are not
orthogonal and will in general have overlap such that
$\textbf{d}(\theta,\alpha,m_{a})\cdot\textbf{d}(\theta^{\prime},\alpha^{\prime},m_{a})\neq
0$.

In principle, this means an axion detection would allow us to determine likely
values of the pulsar parameters in analogy to the way in which observable
gravitational wave signals allow inference of the mass and spin of their
associated black holes. However, our present approach will be to minimise the
expected signal over a conservative subset of values of $(\alpha,\theta)$,
thereby obtaining the most conservative constraints on $g_{a\gamma\gamma}$ for
allowed values of angles.

In order to turn our continuous axion signals into discrete vectors we must
perform some kind of coarse-graining. We therefore define a binning scheme for
$N$ time-channels centered on the points $t_{i}=(i-1/2)\Delta t$ where
$i=1,..,N$ and $\Delta t=1/N$ so that the discretized signal is given by

$S^{\omega}_{t_{i}}=\frac{1}{\Delta t}\int^{t_{i}+\Delta t/2}_{t_{i}-\Delta
t/2}dtS(t,\omega),$ (16)

where $S(t,\omega)$ is the flux density of the axion signal as a function of
time (and frequency) that we derive using our radio signal models.

###  III.2 Why do a time domain analysis?

In this section we will discuss the advantages of doing a time domain analysis
in terms of increasing the chances of detecting axions. We will also discuss
the issue of whether subtraction of the pulse-average from the time-domain
data (as is often the case in pulsar observations and as we have done in our
data) will have significant impact.

In order to do this we will investigate some properties of the matched filter.
Consider applying the matched filter to a given frequency channel whose signal
consists of $N$ time-channels is

$\textbf{S}=(S_{1},\dots,S_{N}).$ (17)

Then according to Eqs (5)-(7) and (14), the matched filter will return

$\braket{\hat{q}}=\frac{\sqrt{\textbf{S}\cdot\textbf{S}}}{\sigma_{N}},$ (18)

where $\sigma_{N}$ is the noise on each of the $n_{t}=N$ time-channels. The
noise amplitude $\bar{\sigma_{N}}$ averaged across all times is then given by
$\bar{\sigma}=\sigma_{N}/\sqrt{N}$. We can therefore re-write (18) as

$\braket{\hat{q}}_{\rm
time}=\frac{\left(\sigma_{S}^{2}+\mu_{S}^{2}\right)^{1/2}}{\bar{\sigma}}$ (19)

where $\mu_{S}=\Sigma_{i}S_{i}/N$ and
$\sigma_{S}=\sqrt{\Sigma_{i}S_{i}^{2}/N-\mu^{2}_{S}}$ are the average and
standard deviation of S, respectively.

Now let us consider carrying out a measurement with no time resolution. This
is the case when one simply uses the telescope to make a total flux
measurement over a long observing time. In this case the noise is again
$\bar{\sigma}$ given by averaging the noise over all integration time, the
signal $\sqrt{\textbf{S}\cdot\textbf{S}}$ is simply given by the mean
$\mu_{S}$ so that

$\braket{\hat{q}}_{\rm flux-avg.}=\frac{\mu_{S}}{\bar{\sigma}}.$ (20)

This can be thought of as a trivial matched filter with a single time channel,
in which any fine-grained time information has been lost. This is in effect
how all previous single pulsar observations for axion dark matter have been
carried out [65, 66, 67].

By comparing the cases of a time-domain analysis (19) with a total-flux
measurement (20) it becomes immediately apparent that since
$\sigma^{2}_{S}\geq 0$, the time-domain analyses will always equal or
outperform the time-averaged measurement. Thus time-domain information
increases the potential to detect axions. In particular, when the relative
time-variation is large ($\sigma_{S}/\mu_{S}\gg 1$), the time-domain search
provides a gain in sensitivity by increasing the signal to noise by a factor
$\simeq\sqrt{\sigma_{S}/\mu_{S}}$, i.e. the square-root of the relative
variance. This is a simple consequence of the fact that $\hat{q}$ is
proportional to the root-square of the signal, and so it implicitly encodes
information about its variability. The same observation was also made in [57]
but the matched filter allows this to be justified from first-principles.

Although this is in general not necessary, the observations used later in the
paper have had the average of the signal subtracted 333Observations of pulsars
are often taken in this way since they are probing the time-variable signal
from the rotation of the neutron star. Ultimately, there is nothing to prevent
the use of the averaged signal as well, but it can require extra work to
calibrate it. Here, we see that if the expected time variation is significant,
removing the average does not significantly affect sensitivity to axions..
Thus while we retain time-domain information, the baseline $\mu_{S}$ of the
signal is essentially re-normalised to zero. In this case, the time-domain
analysis gives a baseline subtracted (BS) value

$\braket{\hat{q}}_{\rm time-BS}=\frac{\sigma_{S}}{\bar{\sigma}}.$ (21)

Clearly, when there is only a small time variation, ($\mu_{S}\gg\sigma_{S}$),
subtracting the baseline leads to a lower value of $\braket{\hat{q}}$,
relative to retaining it as in eq. (19). This is for the simple reason that if
the average signal is large, one loses a lot of signal power by removing the
average. Conversely, in the regime where there is large time variation,
$\mu_{S}\ll\sigma_{S}$ the time-domain analysis with baseline subtraction
performs almost as well as eq. (19). This is again the statement that high
peaks above the noise in the time-domain still allow for good discrimination
from noise.

###  III.3 Demonstration on simulated data

In the previous subsections we have derived the main properties of the matched
filter estimate of the signal-to-noise ratio. This is not particularly new to
those with a background in astronomy, but may not be familiar to a general
reader. The matched filter is the optimal filter for a Gaussian likelihood and
would be the clearest way to identify a signal in the data with ${\hat{q}}$
being a proxy for the signal to noise of detection.

In this subsection, to aid understanding, we will demonstrate the performance
of the matched filter in a toy example by injecting a signal for a neutron
star with the same physical characteristics as PSR J2144$-$3933 into simulated
data with similar noise to the observations that we have in hand, described in
sections IV and IV.1.

We will inject signals with a mass of $m_{\rm a}=4.2\,\mu{\rm eV}$ which
corresponds to an observing frequency $f_{\rm obs}\simeq 1.0{\rm GHz}$ and
$g_{a\gamma\gamma}=10^{-10}{\rm GeV}^{-1}$. We will consider two different
choices for the angles $\alpha$ and $\theta$. Case A with
$(\alpha,\theta)=(0^{\circ},60^{\circ})$ is close to an aligned rotator and
hence we would expect no time variation, whereas case B with
$(\alpha,\theta)=(40^{\circ},60^{\circ})$ has a very strong time variation. We
will consider two observations: one in which the pulse-averaged signal power
is retained, and another where it is removed which is more common in pulsar
observations as we have explained earlier. We have already discussed the pros
and cons of the two approaches in section III.2 and this is just an
illustration of the specific point. The full results of these test cases are
shown in Fig. 2.

For a nearly aligned rotator ($\alpha=0^{\circ}$), the pulsar is axisymmetric
about its rotation axis, and the signal has no time-dependence. Therefore, in
this case, if one searches for the signal with the pulse-average removed,
there is by definition no signal present in the effective data vector, leading
to a non-detection. The filter is essentially scanning a particular noise
realisation with zero signal. In the bottom panel, the input signal contains
significant time-dependence (roughly an order of magnitude). Therefore, the
input signal is detected with a SNR of the same order of magnitude as in the
total power case, in accordance with the discussion comparing (19) and (20).
In all the cases except the baseline-subtracted $\alpha=0$ case, the filter
successfully returns the maximal SNR for the input value of $\theta$.

Figure 2: An illustration of using the matched filter SNR to search for the
signal. $\hat{q}$ is shown as function of observing angle $\theta$ returned by
the matched filter (14) with simulated Gaussian noise similar to that expected
for the observations of PSR J2144$-$3933 discussed in this paper. We display
the SNR for two scenarios for the input signal, one with
$(\alpha,\theta)=(0^{\circ},60^{\circ})$ (case A, top panel) and another with
$(\alpha,\theta)=(40^{\circ},60^{\circ})$ (case B, bottom panel). In both
cases we choose a noise amplitude of $\sim 2.5\,{\rm mJy}$ consistent with
that observed in our data for $m_{\rm a}=4.2\,\mu{\rm eV}$. We use an axion
mass of $3.7\,\mu{\rm eV}$ and $g_{\rm a\gamma\gamma}=10^{-10}\,{\rm
GeV^{-1}}$. As an additional sanity check, we repeat the analysis by inserting
the signal into the real pulse-subtracted data (see sec. IV for details on the
data), represented by the dotted lines. The fact these lines appear to agree
indicates that our noise model is representative of the situation in the real
data.

##  IV Observations of PSR J2144-3933 with MeerKAT

In order to test this idea we selected PSR J2144$-$3933\. We did this by
considering the list of observed neutron stars
[96]444https://www.atnf.csiro.au/research/pulsar/psrcat/ which provides
estimates for $B_{0}$, $P$ and the pulsar distance $D$. In order to make an
estimate of the strength of the signal expected for a particular pulsar we use
an analytic formula based on the radial trajectories approach [55]. Although
this assumption has been shown to be not sufficiently correct to provide
accurate predictions [58, 59] it is likely that it gives a reasonable figure
of merit since it should have the correct scaling with the important
parameters. We will discuss the issue of what is the optimal target in more
detail again in section V in light of what we have learnt. Specifically, we
have used the figure of merit

$\hbox{FOM}=\rho_{\rm DM}{B_{0}^{2/3}P^{7/3}\over D^{2}}\,,$ (22)

where $\rho_{\rm DM}$ is the density of dark matter expected in the vicinity
of the pulsar, to create a ranked list of pulsars. This formula can be derived
from results presented in [55, 60]

We presume that all the dark matter in the Galactic halo is in the form of
axions, the standard assumption when obtaining limits, and extrapolate the
local density of $\rho_{\rm DM}\approx 0.45\,{\rm GeV}\,{\rm cm}^{-3}$ using
an NFW profile for dark matter in the galaxy. Except in the very centre, near
the location of the Galactic Centre Magnetar (GCM) PSR J1745$-$2900555Most
attempts to constrain axions using neutron stars have used the GCM as their
target, attracted by the large magnetic field, and indeed we re-visit using it
for this type of analysis in section V. Depending on the parameters of the NFW
profile we used it varied from around 10th in our list to 1st. A clear reason
to not use it is that there is an additional uncertainty created by the lack
of knowledge of the dark matter density in the centre of galaxy., this is
likely to give a reasonable estimate of the trade-off between $\rho_{\rm DM}$
and $D$ in the FOM.

The PSR J2144$-$3933 which has $B_{0}\approx 2.1\times 10^{12}\,{\rm G}$
estimated from the $P$ and ${\dot{P}}$ based on electromagnetic spin down,
$P=8.51\,{\rm sec}$ and $D=0.12\,{\rm kpc}$ came third on the list 666The
objects PSR J0736$-$6394 and PSR J1856$-$3754 came first and second on the
list. The first is Rotating RAdio Transient (RRAT) which is not monitored
routinely in the radio waveband, while the second has only been detected in
the X-ray waveband. Both these objects are unsuitable for our study here. and
seems an ideal object. This object has a long period, and hence a strong axion
signal, but is otherwise unremarkable. The fact that it is very nearby is also
a significant advantage since it means that we can be more sure about the
local value of $\rho_{\rm a}$ used in our predictions.

Within the GJ model there is a maximum axion mass [67] given by

$\displaystyle m_{\rm a}^{\rm max}\approx 85\,\mu{\rm eV}\left({B_{0}\over
10^{14}{\rm G}}\right)^{1/2}\left({P\over 1\,{\rm sec}}\right)^{-1/2}$ (23)
$\displaystyle\times\left(1+{1\over 3}\cos\alpha\right)^{1/2}\,,$ (24)

which is $\approx 4.7\,\mu{\rm eV}$ corresponding to $f_{\rm obs}\approx
1.15\,{\rm GHz}$ for this object. It is not possible to use any observations
above this frequency in obtaining a limit using the GJ model predictions, but
we do use the data in our search for generalised periodic signals in section
VI.

The specific observation of PSR J2144$-$3933 used in this work was taken at
2020-07-13 02:20:47 as part of the MeerTime Large Survey Programme on MeerKAT.
The observation was recorded as part of the Thousand Pulsar Array [97] census
observations and hence used the ‘full’ MeerKAT array, specifically in this
instance 58 of the antennas were used to form a single tied array beam pointed
at the pulsar. For long-period pulsars the Thousand Pulsar Array census aims
to record $\gtrsim 512$ pulses from each pulsar, and hence the total observing
duration was $4416\,\mathrm{s}$, much longer than typical Thousand Pulsar
Array observations. The data produced by the MeerKAT beamformer are processed
in real time by the PTUSE instrument [98], folding the data with the known
period of the pulsar. Post processing, including initial automated cleaning of
radio frequency interference and flux calibration is carried out on the
Swinburne OzStar supercomputer using the MeerPipe pipeline developed by
MeerTime. The calibration and cleaning procedure used for the Thousand Pulsar
Array data is described in [99]. The output data have 1024 rotational phase
bins and 928 frequency channels, each of width $0.8359375\,$MHz (total
bandwidth $775.75\,$MHz), and centred at $1283.58203125\,$MHz.

###  IV.1 Modelling the noise

Figure 3: Measurements of the mean, $\mu_{S}$ and $\sigma_{S}$ standard
deviations of the observed off-pulse flux density of PSR J2144$-$3933 used in
this paper. The average is over pulsar phase for a fixed frequency. In the top
panel, we show the $\sigma_{S}$ and $\mu_{S}$ for the full data-set. In the
bottom panel, we show the equivalent after employing a cut of
$\sigma_{N}=3.0\,{\rm mJy}$ but with a different scale. As demonstrated in the
figure this cut allows us to excise the frequency channels that are dominated
by RFI contamination, with the remaining channels being compatible with
$\mu\approx 0$ and $\sigma_{N}\approx 2.8\,{\rm mJy}$.

Figure 4: Dynamical spectra (i.e. the data as function of frequency and pulsar
phase) with the RFI dominated channels excised and also the pulse (around
$\phi=0$ and $\phi=1$) removed. The top panel is our noise model described in
the text, while the bottom panel is the actual data.

The observed data has already been processed to remove the effects of Radio
Frequency Interference (RFI). This is sufficient to locate the main peak of
the pulsar pulse, which is typically much stronger than the axion signal. The
first thing we do is remove the pulsar main beam signal from the time-domain,
so that the remaining data is in the off-phase of the pulsar. We do this by
excising 20 time channels from our data. In the top panel of Fig. 3, for the
remaining data we present the average over the pulsar phase for, $\mu_{S}$,
and the standard deviation, $\sigma_{S}$ of the data as a function of the
frequency. It is clear that, despite this procedure, there remains some low
amplitude RFI in certain frequency channels and it is clear that the frequency
channels affected by this must be discarded for the purposes of locating axion
signals. This RFI is typically due to mobile phone, $f\sim 0.95\,{\rm GHz}$,
and Global Navigation Satellite System (GNSS), $f\sim 1.2\,{\rm GHz}$ and
$f\sim 1.6\,{\rm GHz}$777We note that Karoo site where the MeerKAT telescope
is situated, is one of the cleanest radio observation sites in the world and
still there is low-level RFI in these bands which would likely be extremely
difficult to remove, meaning that there are some ranges of axion mass that
will remain unattainable to this technique, and indeed any other searching for
astrophysical axion signals..

This residual RFI can be removed by excising any data with $\sigma_{N}>3\,{\rm
mJy}$ as seen Fig. 3 with the excised data presented using a narrower flux
scale. This data appears to be relatively clean and free from obvious
terrestrial RFI since it removes the regions known to be affected by known
irreducible interference. Once this is done it seems reasonable to try to
model the noise in the data to be an uncorrelated Gaussian random process with
zero mean within each channel and standard deviation $\sigma_{N}(f)$ which is
given by that measured in a given channel. This is a small variation - by
allowing the standard deviation to vary with frequency - to the approach we
have used in section III.1. It is unlikely that the assumption of exact
Gaussianity and zero-correlations is entirely perfect, and indeed in the
subsequent sections we find some evidence to suggest that there may be weak
correlations in the noise as a function of the pulsar phase $\phi$.
Nonetheless, we will argue that this just leads to conservative constraints
and therefore, we will proceed to use this model. We note that if we are able
to accurately model the correlations in the data, this could be handled by the
match filter approach.

In what follows we will use this noise model to obtain limits on axion signal
and, therefore, we should examine to what extent our data resembles Gaussian
noise for the full dynamical spectrum, which is the term used to talk about
the data as a function of frequency, $f$, and pulsar phase, $\phi$. As a self-
consistency test of this noise model, we compare the real data set with that
randomly generated from this distribution and this is presented in Fig. 4.
Visually, the two datasets appear to be very similar and on that basis we
conclude that the models are compatible with each other.

###  IV.2 Constraining the magnetic orientation $\alpha$ and observing angle
$\theta$

We have already pointed out that the amplitude of the time dependence of the
signal depends on the values of $\alpha$ and $\theta$ \- this is also an issue
when using the time averaged signal (see [67], for example). In particular, we
have seen that there can be very little time variation when these angles are
small. Therefore, we will need some further information on the pulsar geometry
to enforce a constraint on $g_{a\gamma\gamma}$.

Obtaining precise values for the pulsar geometry relies in general on strong
assumptions on the observed properties of the neutron star’s radio pulse (e.g.
[100]). Fortunately, from the point of view of the present discussion, we only
need to rule out small angles, and when one observes a narrow pulse profile -
which is the case here - it is unlikely that the magnetic and observation axes
are aligned with the spin axis. In what follows we will describe a simple
model for the pulsar beam geometry with very conservative assumptions and use
it in the case of PSR J2144$-$3933 to argue that one can ignore the region of
parameter space around $\theta\approx\alpha\approx 0$.

Let $W$ be the pulse width corresponding to a fully illuminated circular
radiation beam with half-opening angle $\rho$. These parameters can be related
to the parameters in our misaligned rotator model for the neutron star
$(\alpha,\theta)$ using [101, 102]

$\displaystyle\cos\rho$
$\displaystyle=\cos\alpha\cos\theta+\sin\alpha\sin\theta\cos(W/2)\,.$ (25)

The width of the profile for PSR J2144$-$3933 is measured at 10% of the
amplitude is $(2.1\pm 0.2)^{\circ}$ [103]. It is possible that the profile is
asymmetric and hence assuming that the full open field line region is active
is not necessarily true. This means that the $W$ in (25) should be interpreted
as the pulse width that would be observed if the full beam is active. In rare
cases, the middle of the open field line region is centred in between one of
the profile peaks, and part of the otherwise maybe double profile is missing.
Based on these two caveats, we conservatively take $W$ to be in the range
$1.5^{\circ}$–$5.0^{\circ}$ for this object.

Figure 5: The constraint on the $\alpha-\beta$ plane based on the geometry of
the beam where $\beta=\theta-\alpha$. As explained in the text, we define the
likelihood to be the percentage of possible values of $W$ for which solutions
exist within the range of $h_{\rm em}$ we have allowed. Improbable solutions
with $\left|\beta/\rho\right|>0.95$ are rejected. The top panel shows the
vertical integration of the bottom panel and it demonstrates that solutions
with $\alpha<20^{\circ}$ are very unlikely because they would require fine
tuning of $\beta$. Therefore, we rule out $\alpha<20^{\circ}$ and
$|\beta|>4^{\circ}$ as statistically unlikely, indicated by the red-dashed
lines, when we calculate our limit on $g_{\rm a\gamma\gamma}$.

We now turn to the estimation of $\rho$, whose uncertainty mainly stems from a
lack of knowledge of the height $h_{\mathrm{em}}$ at which the emission
occurs. The beam is bounded by the tangents to the last open magnetic field
lines. Assuming that the field is dipolar, one finds that (e.g. [104])

$\displaystyle\rho$
$\displaystyle\approx\frac{3}{2}\sqrt{\frac{h_{\mathrm{em}}}{R_{\rm
c}}}=\sqrt{\frac{9\pi h_{\mathrm{em}}}{2cP}}.$ (26)

where we have used the small angle approximation for $\rho$. In the second
expression we have replaced the light cylinder radius $R_{c}$ with the pulse
period $P=2\pi R_{\rm c}/c$. Estimation of $h_{\rm em}$ is complicated by the
fact that the beam is not necessarily filled, but across the pulsar population
$h_{\rm em}$ at 1.4 GHz has been constrained to be in the range of
$200$–$400\,{\rm km}$ irrespective of pulse period [105]. We take a more
conservative range of $100\leq h_{\rm em}\leq 1000\,{\rm km}$, so as to ensure
that we are not strongly wedded to the modelling assumptions in the pulse-beam
simulations.

Note that, in principle, it is possible to compute the emission height $h_{\rm
em}$ from estimates of the swing of the polarisation angle (PA) (see for
example [106]) where theoretical PA profile is derived assuming the rotating
vector model [107]. However, the detailed applicability of the rotating vector
model is a debated topic in the literature and not all pulsars exhibit the
characteristic banana shape that the model predicts (see for example [108]).
In addition, it may not always be possible to measure the polarisation
properties of candidate pulsars, particularly magnetars that do not emit at
radio frequencies. We remark that while measuring $\alpha$ can significantly
impact the sensitivity of our method, it remains to be seen whether this
method can be considered an guaranteed way of breaking the $\alpha-\theta$
degeneracy.

Since parameters for a given pulsar are uncertain, values of $\alpha$ all the
way down to zero are allowed by (25), for which there would be no time-
variation. We, therefore, appeal to further arguments which allow us to place
a lower bound on $\alpha$ by excluding implausible geometries.

A problem with very small $\alpha$ geometries is that in order to explain the
very narrow $W$ the line-of-sight needs to graze the very outer part of the
beam such that most of the beam is invisible to us. This requires fine tuning
which is not only unlikely [109], but is also contrived for two reasons. First
of all, a small change in emission height as expected for different observing
frequencies [110] would lead to a drastic change in the observed pulse width,
which is not observed [103]. Secondly, a narrow pulse not only requires a
grazing line of sight, but also a circular beam with a hard edge. In reality
the pulsar beam does not have a hard edge, and hence the observed pulse shape
from a grazing line of sight will be dominated by the intrinsic smoothness of
the beam which will be much wider than predicted from the circular model.

To quantify why small $\alpha$ geometries are unlikely, we construct an
effective probability distribution for $\alpha$ and $\beta$ which essentially
measures the number of beam realisations associated to each pair
$(\alpha,\beta)$ assuming a uniform distribution for $W$. This is shown in
Fig. 5. We constructed this distribution according to the following algorithm
(i) uniformly sample W between $1.5^{\circ}$ and $5.0^{\circ}$ in 100 steps.
(ii) For each given W, scan over a discrete grid of $(\alpha,\beta)$ values
between $0$ and $\pi/2$. (iii) For each point $(\alpha,\beta)$ calculate
$\rho(W,\alpha,\beta)$ with Eq. (25). For each $(\alpha,\beta)$ (for the
specific $W$ under consideration) if the resulting $\rho$ satisfies $100\leq
h_{\rm em}\leq 1000\,{\rm km}$ and $|\beta/\rho|<0.95$ record a value of 1.
Otherwise assign it $0$. Note the second constraint is designed to exclude a
line-of-sight with an impact parameter $\beta\simeq\rho$, which is both
unlikely and implausible for the reasons above. At the end of this process,
for each $W$, one has an $\alpha$-$\beta$ grid with entries that are 1
(implying an acceptable beam geometry exists) or 0. Since $W$ is sampled with
100 steps, there are 100 grids. (iv) Fig. 5 then shows the sum of these grids
(appropriately normalised).

The main features that stand out are that $|\beta|$ needs to be small in order
for the line of sight to intersect the beam and small $\alpha$ solutions are
excluded. No solutions exist for $\alpha\lesssim 10^{\circ}$. Furthermore, one
can see in the top panel that solutions $\alpha\lesssim 20^{\circ}$ are
unlikely geometrical solutions, which is because they require a fine tuned
(large) $\beta$. It should be stressed that $\alpha\lesssim 20^{\circ}$
geometries are not just unlikely, we have also argued them to be contrived888A
less conservative limit on the acceptable solutions such that
$|\beta/\rho|<0.90$ would lead to $\alpha\gtrsim 25^{\circ}$, demonstrating
that the conclusions are relatively insensitive to the precise limits chosen..
In what follows we will impose $\alpha>20^{\circ}$ and $|\beta|<4^{\circ}$. We
would expect to be able to apply similar arguments to a large fraction of
pulsars that we might want to use to constrain $g_{a\gamma\gamma}$.

###  IV.3 Constraints on $g_{a\gamma\gamma}$ from PSR J2144$-$3933

We have shown in section III how one can compute the signal-to-noise (SNR)
parameter $q$ as a function of the input parameters $(m_{\rm
a},\theta,\alpha)$ using the matched filter. In order to now derive
constraints, we have carried out a parametric search for all possible profiles
in our interpolated library (see Fig. 1) using the pulsar data presented in
Fig. 4. This procedure then gives us a distribution of SNR values $q_{\rm
meas.}$ associated to each profile. We could repeat this process, but this
time using our noise model which by definition has no signal present. Remember
that it assumes uncorrelated Gaussian noise which simplifies the matched
filter. This means that the values of $q_{\rm meas.}$ should be Gaussian
distributed with zero mean and unit variance.

We do find that they are compatible with a Gaussian distribution that has zero
mean. However, we find that the standard deviation is $\approx 1.45$ somewhat
higher than the expected value which points to the fact that the noise model
we are using is not optimal 999We have assumed that he noise is uncorrelated
in phase which is is unlikely to be precisely true and this could easily lead
to more structure than would be expected for a totally random realization of
the noise. It is clear that one could achieve slightly tighter constraints by
improvement of the noise model.. We find that, out of the $\approx 3\times
10^{4}$ templates, there are three that have $q>5$ which, if the noise model
were perfect would suggest candidate detections, but in order to assess their
statistical significance they should probably be scaled down by $1/1.45$
reducing them to $\approx 3.5$. This suggests that they are chance alignments
with the templates; a conclusion that is further strengthened by the
observation that they, and indeed the other higher values of $q_{\rm meas.}$,
appear to be randomly distributed with $m_{\rm a}$. We are satisfied that our
data are compatible with a null detection.

In the case where the baseline has been subtracted from the data, which is the
case for the data being considered here, the constraining power of the matched
filter is determined by (21). For the pulsar we have chosen, Fig. 1 shows the
time-dependence of the profiles. Note the relative variance vanishes at
$\alpha=0$, but can be larger than one for a range of values of
$(\alpha,\theta)$, but that we have argued that regions with
$\alpha<20^{\circ}$ are unlikely and contrived in section IV.2.

In the absence of a detection, one can derive constraints on
$g_{a\gamma\gamma}$ by comparing the measured value of $q_{\rm meas.}$ for a
particular template, with the expected value $\braket{q}$ for that same
template. The measured value is defined by

$q_{\rm
meas.}(\theta,\alpha)=\frac{\textbf{F}(\theta,\alpha)\cdot\textbf{d}}{\sigma_{N}}.$
(27)

This must be compared to the expected value $\braket{q(g_{\rm
a\gamma\gamma},\theta,\alpha)}$ if a signal were present in the data

$\braket{q(g_{\rm
a\gamma\gamma},\theta,\alpha)}=\frac{\sqrt{\textbf{S}(g_{a\gamma\gamma},\theta,\alpha)\cdot\textbf{S}(g_{a\gamma\gamma},\theta,\alpha)}}{\sigma_{N}}.$
(28)

In the perturbative limit of the conversion probability, which we have checked
is always the case here, this is $\propto g_{a\gamma\gamma}$. Therefore, in
order to impose a limit we calculate $\braket{q(g_{\alpha\gamma\gamma}^{\rm
fid},\theta,\alpha)}$ for $g_{a\gamma\gamma}^{\rm fid}=10^{-10}\,{\rm
GeV}^{-1}$ and exclude any value such that

$g_{a\gamma\gamma}>g_{a\gamma\gamma}^{\rm fid}\sqrt{2q_{\rm
meas.}(\theta,\alpha)\over\braket{q(g_{\alpha\gamma\gamma}^{\rm
fid},\theta,\alpha)}}.$ (29)

The factor of two in (29) corresponds to the observed noise level being twice
the equivalent signal for $g_{a\gamma\gamma}$ meaning that this will be a $2$
$(\approx 95\%$ confidence) upper limit on coupling constant.

We have explained above that the typical values of $q_{\rm meas.}$ are higher
than they should be due to the noise model not being perfect. This will mean
that the upper limits we compute using (29) are not optimal and hence are
slightly conservative.

We show the constraints from this procedure in the left panel of Fig. 6 for
$m_{\rm a}=4\,\mu{\rm eV}$ along with the regions preferred by our analysis of
the orientation angles. In the right panel, we quantify the advantage of the
purely time-domain analysis compared to the purely frequency-domain analysis
as a function of $(\alpha,\theta)$. We do this by taking the ratio of
theoretical signal-to-noise in the case where the baseline is subtracted,
$q_{\rm time}$, with respect to the case where the time profile is averaged
over the pulse-period and integrated, $q_{\rm freq.}$. When this ratio is
larger than 1, the time-domain analysis leads to stronger constraints on
$g_{a\gamma\gamma}$.

Figure 6: In the left panel, we present the constraints on the axion-photon
coupling for the average-subtracted case as a function of $(\theta,\alpha)$
for $m_{\rm a}=4\,\mu{\rm eV}$. As expected, the derived limit is weaker in
the limit $\alpha\rightarrow 0$ since the time-dependence of the signal is
negligibly small. In fact there is no limit for $\alpha\equiv 0$ since there
is no time dependence. Fortunately, the limits from the pulsar main beam
modelling require $\alpha>20^{\circ}$ and $|\beta|<4^{\circ}$ which are
included as red lines on the plot, so we are able to achieve a limit
$g_{\gamma\gamma\gamma}\lesssim 9.6\times 10^{-11}\,{\rm GeV}^{-1}$. On the
right panel, we show the ratio of the signal-to-noise in the case where the
average has been subtracted (i.e., where only time-domain data is used) the
case where the time-averaged flux is used in frequency space. In other words,
this ratio quantifies the gain in working in the time-domain.

We have performed the same analysis as presented in Fig. 6 over the mass range
$3.9\,\mu{\rm eV}\leq m_{\rm a}\leq 4.7\,\mu{\rm eV}$ and then searched for
the weakest upper limit in the range $\alpha>20^{\circ}$ and
$|\beta|<4^{\circ}$. The results, assuming that $D=0.165\,{\rm kpc}$, are
presented in Fig. 7 and using this, we conclude that we can exclude dark
matter axions forming the entire Galactic halo with
$g_{a\gamma\gamma}>5.5\times 10^{-11}\,{\rm GeV}^{-1}$. If $D=0.29\,{\rm kpc}$
as suggested by electron density models, we find a weaker constraint of
$g_{\rm a\gamma\gamma}>9.6\times 10^{-11}\,{\rm GeV}^{-1}$. We leave the
detailed investigation of the impact of the galactic electron density models
and the veracity of the models for future work.

Figure 7: $2$ upper limits on $g_{a\gamma\gamma}$ as a function of $m_{\rm a}$
for $d=0.165\,{\rm pc}$. This is determined by calculating the highest upper
limit in the region of the $\alpha-\theta$ plane allowed by
$|\beta|<4^{\circ}$ and $\alpha>20^{\circ}$ from the equivalent of Fig. 6. On
the basis of this figure we quote an upper limit of
$g_{a\gamma\gamma}<5.5\times 10^{-11}\,{\rm GeV}$ over the mass range
$3.9\,\mu{\rm eV}\leq m_{\rm a}\leq 4.7\,\mu{\rm eV}$.

##  V Future Searches

In the previous section we have shown how one can derive a limit on
$g_{a\gamma\gamma}$ from baseline subtracted radio pulsar data using the
variation in the time domain calculated by our ray-tracing algorithm. There we
used a specific pulsar and telescope. An obvious question is what can be
gained in future observations. In general, the limits obtained will scale as

$g_{a\gamma\gamma}^{\rm
lim}\propto\frac{1}{(\mu_{S}^{2}+\sigma_{S}^{2})^{1/4}}\left({A_{\rm eff}\over
T_{\rm sys}}\right)^{-1/2}t_{\rm obs}^{-1/4}.$ (30)

Hence, any improvement on $g_{a\gamma\gamma}^{\rm lim}$ will come from two
avenues. The first is better observations: lower system temperature $T_{\rm
sys}$, larger collecting area $A_{\rm eff}$ and increased observing time,
$t_{\rm obs}$. The second is a better target: one with larger mean-signal
power $\mu_{S}$ or greater time variability as measured by $\sigma_{S}$.
Leveraging the latter is of course the key point of the present paper. Let us
now come to each of these factors in turn.

In terms of improved observations with the current target, we can imagine
future observations of PSR J2144$-$3933 with both MeerKAT or other similar
telescopes in the short term and the Square Kilometre Array (SKA) in the
longer term. All other things being equal, eq. (30) implies that an increase
in observation time from $\approx 1$ hour, as we have now, to 100 hours would
improve the limit by a factor $\sim 3$. Moreover, going from MeerKAT with
$A_{\rm eff}/T_{\rm sys}\approx 450$ to $A_{\rm eff}/T_{\rm sys}\approx 1800$
for SKA1-mid will yield a limit which is a factor $\sim 2$ better. Combining
these together one might be able to obtain a limit of
$g_{a\gamma\gamma}<1.6\times 10^{-11}\,{\rm GeV}^{-1}$.

PSR J2144$-$3933 only allows constraints to be imposed for $m_{\rm
a}<4.7\,\mu{\rm eV}$ within the framework of the GJ model. It might be
interesting to perform lower frequency observations of it, but perhaps it is
more interesting to find a source with a larger value of $m_{\rm a}^{\rm max}$
to probe mass ranges less accessible to terrestrial haloscopes. An interesting
point to clarify would be how the signal scales as a function of pulsar
parameters, in particular the period $P$ and the magnetic field $B_{0}$.

The figure of merit (FOM) based on radial trajectories is $\propto
B_{0}^{2/3}P^{7/3}$, while the maximum mass probed is $\propto(B_{0}/P)^{1/2}$
101010Note that this figure-of-merit is based on the calculation of power
radiated at the poles for an aligned rotator Hook et al [55]. We think that
this choice is a conservative estimate of the FOM. Indeed a time-varying FOM
can be computed in the same approach, but we think the resultant ranking for
targets from such choices for the FOM are subject to errors arising from the
radial trajectories assumption that has been superceded by the ray-tracing
approach. . At a fixed value of $P$ we see that it will always be best to
increase the value of $B_{0}$, while at fixed $B_{0}$ the FOM will increase
with $P$, but $m_{\rm a}^{\rm max}$ will decrease meaning that there is a
trade-off between the two and the optimal target will be a compromise between
the strength of the signal and the range of mass probed.

We have checked the scaling of the FOM using the ray-tracing code to calculate
the average signal, $\mu_{S}$, summed over all frequencies for a range of
values of $B_{0}$ and $P$. For $B_{0}$ these appear to broadly confirm that
the approximate scaling of the FOM with relatively weak dependence $\theta$
and $\alpha$, although there can be significant deviations for extreme values.
However, the dependence on $P$ seems to be somewhat weaker than that predicted
from the FOM from radial trajectories. For the specific choice of
$\alpha=60^{\circ}$ we find that $\hbox{FOM}\propto B_{0}^{0.8}P^{1.2}$. Given
the complications in constructing a FOM which depends on the orientation
angle, we conclude that it is reasonable to continue to use the FOM based on
radial trajectories as a rule of thumb, but we should not expect it to give
quantitatively accurate predictions for the increase in constraining power.
The argument above applies to all attempts to constrain $g_{a\gamma\gamma}$
using neutron stars, not just a search for time dependent signals.

Let us come now to the question of other targets, focusing in particular on
the time-dependence of their signals. Note that in our original target
selection, we used only the mean power, however as demonstrated in Sec. III.1
and as is apparent from Eq. (30), the key parameter determining whether
including a time-domain analysis can add value over just using the frequency
domain is $\sigma_{S}/\mu_{S}$. Clearly this will depend on the intrinsic
pulsar parameters, $B_{0}$, $P$ and axion mass, $m_{a}$ as well as
$(\alpha,\theta)$. In Fig. 8 we show that $\sigma_{S}/\mu_{S}$ clearly
increases with $B_{0}$ with the amount being sensitive to the choice of
$\theta$ and $\alpha$. We have also investigated the dependence on $P$, but
this is typically much weaker for the relevant range of parameters. Based on
this, we can expect that time-domain observations offer the greatest
enhancement over a total flux measurement for larger values of $B_{0}$. Hence,
objects such as magnetars stand to gain most from a time-domain versus total
flux analysis.

The GCM in particular has already been a popular target for searching for
axion signals. In Fig. 9 we present some GCM profiles analogous to those in
Fig. 1. The first thing to notice is that the signal is quite a bit larger,
mainly due to the enhanced dark matter density assumed in the GC, although
this is mitigated somewhat by it being further away.

Figure 8: The relative time variance, $\sigma_{S}/\mu_{S}$, of the profiles as
a function of the magnetic field $B_{0}$ at the surface of the neutron star.
We fix $P=4\,{\rm s}$ and $m_{\rm a}=1\,{\rm\mu eV}$ and the value of
$g_{a\gamma\gamma}$ scales out in this ratio. Despite the time-variation
between the maximum and minimum of the profiles increasing by orders of
magnitude for $B_{0}\sim 10^{14}\,{\rm G}$ compared to $B_{0}\sim
10^{12}\,{\rm G}$, $\sigma_{S}/\mu_{S}$ only goes up by a factor of a few.

Clearly, at least for some choices of the orientation angles the profiles are
substantially more localised in pulsar phase - they are almost pulse-like, but
they are still much narrower than the width main peak in the radio pulse
profile for this object. The strong dependence on pulse phase in this case is
due to the effects of the “throats” in the magnetosphere where axion
production is enhanced. This is likely to also enhance the constraining power.
Obviously, this is a model-dependent assumption that arises from assuming the
GJ model. We emphasise the point that the strength of this technique is
directly proportional to the size of the time-dependence, and by extension the
presence of these throats in the charge density.

Figure 9: Predicted pulse profiles for the GCM. In order to generate these
profiles, we fix $g_{a\gamma\gamma}=10^{-10}\,{\rm GeV}^{-1}$ and $m_{\rm
a}=1\,\mu{\rm eV}$ and $\alpha=60^{\circ}$ while varying the observing angle.
As in Fig. 1 we have removed the average signal from the profiles. We can see
very clearly that in some cases the profiles are very time variable - they
almost look pulse-like for $\theta=30^{\circ}$ and $60^{\circ}$. Nonetheless
these only correspond to a $\sigma_{S}/\mu_{s}\approx 3$ compatible with Fig.
8.

Further to these points, the GCM offers both large relative time-variance and
it lies in a region of larger dark matter density. Furthermore the axion-to-
photon conversion is enhanced due to large magnetic fields. More specifically
it has a large magnetic field $B_{0}\approx 1.4\times 10^{14}\,{\rm G}$ and a
period of $P=3.76\,{\rm sec}$ and is a distance of $D\approx 8.3\,{\rm kpc}$.
In fact, the ratio of the FOM Eq. (22) from PSR J2144$-$3933 relative to the
GCM is given by

${\hbox{FOM}|_{\rm GCM}\over\hbox{FOM}|_{J2144}}\approx 1.3\times
10^{-4}\times{\rho_{\rm DM}|_{\rm GCM}\over\rho_{\rm DM}|_{\rm local}}\,.$
(31)

If one assumes a standard NFW profile for the galaxy, one obtains an
enhancement of $\sim 10^{5}$ with respect to the local value 111111Note that
without this assumption, the constraints from this pulsar are weaker than PSR
J2144$-$3933. which suggests that it will lead to a larger FOM by a factor of
$\sim 20$. If the dependence on $B_{0}$ and $P$ is slightly stronger than
using the radial trajectory based FOM, as we have suggested above, we might
expect a slightly stronger improvement than given by this simple argument. In
addition, due to the larger value of $B_{0}$ from Fig. 8 we would expect
$\sigma_{S}/\mu_{S}$ to be large enough to produce the most significant gains
from a time-domain study compared to other targets.

Prima-facie there seems to be an argument for reconsidering the GCM as a
target. As an indication of what could be achieved for the GCM with similar
observational resources to those currently available, we have simulated the
equivalent of Fig. 7 using the GCM as the source assuming a similar noise
level to the present data, that is $3\,{\rm mJy}$, and this is presented in
the top panel of Fig. 10121212We note that at 1.4 GHz there is a significant
increase in the sky temperature at the location of the GCM and this will
dominate $T_{\rm sys}$. This means to achieve this noise figure one would need
to observe for $\sim$10 h with MeerKAT. In addition the pulsar and the axion
signal will be scattered [111]. However, as the signal can be completely
modelled using the observed pulsar profile it can be combined with the axion
templates before the matched filtering is performed. Both these effects will
fall off rapidly as a function of increasing frequency.. Making similar
assumptions about constraints on $(\alpha,\beta)$ from the pulse profile, we
obtain a projected upper limit of $g_{a\gamma\gamma}<4\times 10^{-12}\,{\rm
GeV}^{-1}$ which is a factor $\sim 10$ stronger than that we obtained from PSR
J2144$-$3933\. In addition, the GCM allows a much wider range masses with
$m_{\rm a}^{\rm max}\approx 85\mu{\rm eV}$ $(f_{\rm obs}^{\rm max}\approx
20\,{\rm GHz})$. While there are many caveats to such an analysis, notably the
noise levels that one might achieve in the direction of the GC, this suggests
recording time domain information for this object is well motivated. Note that
above we have argued that it might be possible to improve limit by a factor
$\sim 6$ using a 100 hour observation using SKA1-mid, which would lead to a
projected limit of $g_{a\gamma\gamma}<6\times 10^{-13}\,{\rm GeV}^{-1}$.

It is worth noting that the GJ model is believed to be a reasonable
description of the pulsar magnetosphere for typical radio pulsars with
magnetic field strengths in the range $10^{11}\,{\rm G}\leq B_{0}\leq
10^{13}\,{\rm G}$. However, it breaks down in pulsars with large magnetic
field values, magnetars in particular for which the emission mechanism is
poorly understood. Since a detailed analysis of the magnetosphere and the
associated modelling uncertainties and propagating these uncertainties through
our pipeline requires detailed MHD simulations, we leave such efforts to
future work. It is a well-known fact that DM density at the galactic centre is
very poorly constrained. It is our considered opinion that the extrapolation
of the NFW profile to the position of the GCM in the galactic halo to predict
$\rho_{DM}$, which translates to an uncertainty of $\sim 2$ orders of
magnitude in $g_{\rm a\gamma\gamma}$ is by far the biggest uncertainty in our
GCM modelling. In order to address the former problem, we will now describe a
model-independent analysis of the radio pulsar data to look for generalised
periodic signals.

Figure 10: In the top panel, we show the expected constraints on $g_{\rm
a\gamma\gamma}$ from simulated observations of the GCM with an r.m.s. noise
level of $3\,{\rm mJy}$ for $m_{\rm a}=1\,\mu{\rm eV}$. We fix the pulsar
parameters to be $B_{0}=1.4\times 10^{14}\,{\rm G}$, $P=3.76\,{\rm s}$ and
$\rho_{\rm DM}=5.4\times 10^{4}\,{\rm GeV\,cm^{-3}}$ which is the value
computed using a standard NFW profile for the galaxy. Note that, depending on
the values of $\theta$ and $\alpha$ limits as low as
$g_{a\gamma\gamma}\lesssim 10^{-13},{\rm GeV}^{-1}$ might be possible. In the
bottom panel, we quantify the effect of adding time-domain information over
the parameter space $(\alpha,\theta)$. Note that this has a slightly different
morphology to that for PSR J2144$-$3933, but also the values of $q_{\rm
time}/q_{\rm freq}$ are slightly larger approaching $\approx 3$ as it is
indicated might be the case in Fig. 8.

In the bottom panel of Fig. 10, similar to the right panel of Fig. 6, we show
how much one can gain from doing the time-domain analysis as a function of
$(\alpha,\theta)$. We see that the values of $q_{\rm time}/q_{\rm freq}$
larger, typically $\approx 2-3$ compared to $\approx 1$ for PSR J2144$-$3933\.
It seems reasonable to conclude that the factor $\approx 20$ improvement in
the constraining power seen in the case of the GCM comes from the combination
of the FOM based on radial trajectories, a slightly stronger dependence of the
signal on $B_{0}$ and $P$ than predicted by radial trajectories and the use of
the time domain structure. Our conclusion is that there is a strong argument
for attempting to apply this technique to the GCM.

Finally, we comment that it goes without saying that improved knowledge of the
orientation angles will lead to improved constraints. The constraints which we
have imposed on $g_{a\gamma\gamma}$ from PSR J2144$-$3933 are strongly
dependent on constraints we have placed on $\alpha$ and $\beta$ from the pulse
profile. These are conservative and given the nature of the signal - there are
some regions of the $\alpha-\theta$ plane where the constraints are very weak
and even non-existent - meaning one is always dominated by the lower limit one
imposes on $\alpha$. However, if one were to know the actual angles and,
indeed, if they were in the region where the signal is predicted to be
strongest then more optimal limits can be imposed. For example, in the case of
the GCM if the actual angles are in the region where the strongest limit is,
$g_{a\gamma\gamma}\lesssim 10^{-13}\,{\rm GeV}^{-1}$ for a 1 hour observation
with MeerKAT. This could be a factor of 6 lower for 100 hours with the SKA.

##  VI Search for generalized periodic signals

In the previous sections we have deduced limits on $g_{a\gamma\gamma}$ using
PSR J2144-3933 and have discussed how this might be improved in the future
using the GCM. The key qualitative feature of the predicted signal profiles is
that their timescale is given by the pulse period $P$ \- this is what allows
us to use the pulsar data already folded at the pulse period. One might be
concerned that the precise predictions using the GJ model might be too
simplistic given the complicated nature of the pulsar magnetosphere. However,
taking the qualitative prediction that the signal is dominated by low
harmonics of the period one might perform a generalised search for the
periodic signals in the data. Of course, without specific connection to the
physics, such a search cannot yield an upper limit on $g_{a\gamma\gamma}$, but
does allow us to search more specifically for periodic would-be axion signals
which might otherwise be missed due to modelling uncertainties.

Any time-periodic data $d(t)$ can be written as a Fourier series

$d(t)=\sum_{k=-\infty}^{\infty}a_{k}e^{2\pi ikt/P}\,,$ (32)

where $P$ is the signal period. The coefficients $a_{k}$ are then given by

$a_{k}=\frac{1}{P}\int_{0}^{P}dt\,\,d(t)e^{-2\pi ikt/P}\,.$ (33)

Reverting to the discrete case considered in this text, where the data is
defined on $N$ discrete time-bins, with entries $d_{q}$, $q=0,..,N-1$, the
Fourier coefficients in the discrete limit can be written as

$a_{k}=\frac{1}{N}\sum_{q=0}^{N-1}d_{q}\exp{\left(-\frac{2\pi
iqk}{N}\right)}\,.$ (34)

Note that in what follows we do not consider the $k=0$ mode associated to the
time-average $\mu_{S}$ of the profiles which, in principle, has been already
removed by the data processing. In addition, we comment that the parity
properties of the profiles we predict for the axion signal imposes ${\rm
Im}(a_{k})=0$.

In order to confirm our assumption that the profiles are dominated by low $k$
modes we have computed the power spectrum131313We note that $a_{k}$ depends on
the relative phase of the main peak of the pulsar profile and axion signal,
but when one takes the power spectrum this information, which is encoded in
the complex phase of $a_{k}$, is removed. So our test using the power spectrum
is independent of this assumption., $|a_{k}|^{2}$, of two profiles with
varying levels of time-dependence, i.e., a relatively flat profile with
$(\alpha=10^{\circ},\theta=10^{\circ})$ and one with a relatively large time
variation with ($\alpha=60^{\circ},\theta=30^{\circ})$. This is presented in
the in the top panel of Fig. 11. We see that the power spectrum for both
decreases rapidly as a function of the mode number $k$. The point can be
further reinforced by computing the inverse-Fourier transform from the $a_{k}$
while neglecting all modes above some value of $k=n_{\rm cut}$. We have done
this for $n_{\rm cut}=1,3,7$ and present the outcome in the middle and bottom
panels of Fig. 11. It is clear that the profiles are reasonably well-
represented by $n_{\rm cut}=3$ and $n_{\rm cut}=7$ in each case. We find that
this is true for a wide range of predicted templates, but not all e.g. some of
those presented in Fig. 9.

Figure 11: In the top panel we present the power spectrum of two profiles with
with ($\alpha=60^{\circ},\theta=30^{\circ})$ (blue) and
$(\alpha=10^{\circ},\theta=10^{\circ})$ (red). The power spectrum has been
normalised such that $|a_{1}|=1$. It is clear that both the power spectra are
dominated by the lowest $n$ modes. In the middle and bottom panels we
illustrate this point by presenting profiles obtained from the inverse
transform of the FFT of each profile with $n_{\rm cut}=(1,3,7)$ (red, green
and blue, respectively) compared to the exact profile (black). The
significance of the low $n$ modes is further demonstrated by the fact that the
blue curves are remarkably close in shape to the black curves.

Figure 12: In the top panel we present a histogram of the power spectrum
computed from the data for $m_{\rm a}=4\,\mu{\rm eV}$ in orange and a single
Gaussian noise realization of the same variance and sample size. Each of the
individual mode labelled by $k$ is binned, so this corresponds to a single
mode search. We also plot the $\chi^{2}_{1}$ PDF associated expected for
$|a_{k}|^{2}$ and observe that both the data, noise realization and analytic
PDF are clearly compatible with each other. This appears to suggest that the
data are compatible with being pure noise. In the middle and bottom panels we
compare the distribution of the sum of first three modes (middle) and first
seven modes (bottom) from 300 realisations of simulated Gaussian noise with
the same standard deviation as the data. The measured value is represented
with a red line. We use the same value $m_{\rm a}$ as in the top panel and
have also included the theoretical PDFs which are $\chi^{2}_{3}$ and
$\chi^{2}_{7}$ respectively. By doing this, we are searching for signals that
have more structure than just a single mode in the time-domain. Since the data
value is not a significant outlier, we exclude the presence of such signals in
the data.

The advantage of this approach is that we know a periodic signal search can be
carried out using only the first $n_{\rm cut}$ Fourier modes, reducing
computational overhead when scanning.

Furthermore, in what follows, rather than scanning over every available
template parametrised by $(a_{1},\cdots,a_{n_{\rm cut}})$, with say, a fixed
total power (which would be very costly from a numerical point of view) we
instead compare the expected distribution of the $a_{k}$ that would follow if
d were pure Gaussian noise. In that case the expected PDF for the values of
the squared-amplitudes $|a_{k}|^{2}$ is a $\chi^{2}_{m}$ distribution

${\cal P}(x;m)=\frac{1}{2^{m/2}\Gamma\left(\textstyle{1\over
2}m\right)}x^{m/2-1}\exp\left(-\textstyle{1\over 2}x\right)\,,$ (35)

where $\Gamma$ is the standard gamma function. For a given mode, we have that
$m=1$ is the number of degrees of freedom since the noise is real, and hence
the real and imaginary parts are correlated. Note the distribution is the same
for all $k$ meaning the spectrum is scale-invariant. However, if there is an
axion signal in the data, we expect this scale invariance to be broken, and,
in particular, for there to be a greater portion of power contained in the low
Fourier modes, providing a model independent test of periodic axion signals in
the data. Therefore, since the scale-invariance implies all $|a_{k}|^{2}$ have
an identical distribution, by Fourier transforming the data $d(t)$ and binning
the corresponding values of $|a_{k}|^{2}$ across all $k$, we can see to what
extent they are $\chi_{1}^{2}$ distributed. This comparison is shown in the
top panel of Fig. 12, where we present histograms of the data (together with a
single noise realization) for a particular frequency channel corresponding to
a value of $m_{\rm a}=4\,\mu{\rm eV}$, both of which seem compatible with
theoretical $\chi^{2}_{1}$ distribution. The sample is drawn from Fourier
modes up to some $k_{\max}$ where the wavelength of the mode would be less
than the temporal bin-width, beyond which the description breaks down owing to
insufficient resolution.

Having extracted spectral information from the data, we now want to analyse to
what extent power is concentrated in the lower modes, as expected for an axion
signal. To do this, we look at the power contained in the sum of the first
$n_{\rm cut}$ modes, given by

$Q(n_{\rm cut})=\sum_{k=1}^{n_{\rm cut}}|a_{k}|^{2}.$ (36)

We then want to understand if the measured values of low-mode power are again
consistent with the scale-invariant spectrum. The PDF for $Q(n_{\rm cut})$ is
$\chi^{2}_{n_{\rm cut}}$ since it is the sum of $n_{\rm cut}$ independent
modes each $\sim\chi^{2}_{1}$. This distribution can clearly be seen in the
bottom two panels of Fig. 12 where for $n_{\rm cut}=3$ and $n_{\rm cut}=7$ we
display simulated and analytic PDFs for the frequency channel corresponding to
$m_{\rm a}=4\,\mu{\rm eV}$. Note we have also included the measured value of
$Q(n_{\rm cut})$ shown in red.

In order to make a more precise statistical statement about the presence of a
signal dominated by low-Fourier modes, we have calculated the probability,
using the appropriate PDF, for there to be a larger value of $Q(n_{\rm cut})$
than that which is measured. This gives a sense of whether or not the measured
value happens by chance in a way consistent with our sample size, or whether
it is a sufficient outlier to indicate the presence of a periodic signal in
the data. For the case of $n_{\rm cut}=3$ we find that there are three
frequency channels where this “probability to exceed” $<0.05$ and one where it
is $<0.02$. However, since there are 213 individual frequency channels it
seems likely that this is a chance outcome in a few frequency channels
compatible with being a random noise realization. For $n_{\rm cut}=7$ there
are none with a probability $<0.05$.

##  VII Discussion and Conclusions

In this paper we have presented a procedure for carrying out time-domain
searches for radio signals produced by axion dark matter converting into
photons in the magnetospheres of pulsars. We have developed a matched filter
formalism to define the signal-to-noise ratio of time-dependent signals and
have used this to show that time-domain searches always improve the signal to
noise ratio. By how much the SNR is improved is determined by the relative
variance of the signal (the ratio of the standard deviation to the mean of the
signal over the pulse period141414Note that in pulsar astronomy, when applied
to e.g. the main pulse, this quantity is referred to as the modulation index.)
and the matched filter formalism provides a robust framework to understand why
this is the case.

As a test case, we then applied the matched filter formalism to real data on
PSR J2144$-$3933 obtained using MeerKAT, searching for expected periodic
signal templates for the radio signatures produced by axion dark matter. This
was selected from a list of pulsars on the basis of a simple figure of merit
for axion detection. In the present analysis, for a fixed axion mass, these
templates form a two-parameter family for each observing direction set by the
angle $\theta$ between the stars rotation axis and the line of sight towards
the pulsar, and $\alpha$ the angle between the stars magnetic axis and its
rotation axis. Using the morphology of the observed pulsar main-beam signal,
we were able to exclude a range of values $(\alpha,\theta)$, narrowing down
the number of viable templates. Scanning over the allowed set of templates we
find no evidence for axion dark matter and obtain an upper limit of
$g_{a\gamma\gamma}<9.6\times 10^{-11}\,{\rm GeV}^{-1}(D/0.29\,{\rm pc})$ over
the mass range $3.9\,\mu{\rm eV}\leq m_{\rm a}\leq 4.7\,\mu{\rm eV}$. Given
the astrophysical uncertainties in modelling the templates, we also carried
out a generic periodic signal search independent of any modelling. This also
returned no significant signal from axion dark matter.

Figure 13: We show our limit (dark red) on the axion-photon coupling from 1
hour observations of PSR J2144-3933 ($B_{0}\simeq 2\times 10^{12}\,{\rm
G},\rho_{\rm DM}\simeq 0.45\,{\rm GeV}{\rm cm}^{-3}$ and the projections from
the 100 hours of observations from he SKA (red). We also show projections of
what can be achieved from $\sim 1$ hour observations of the GCM ($B_{0}\simeq
1.6\times 10^{14}\,{\rm G}$, $\rho_{\rm DM}\simeq 6.5\times 10^{4}{\rm
GeV}{\rm cm}^{3})$ using MeerKAT at L-band (purple) presented in this paper
and 100 hours with the SKA (black) assuming a frequency coverage $0.9\leq
f_{\rm obs}\leq 15.3$ GHz. Note that our projections for the GCM are derived
assuming $\theta=\alpha=20^{\circ}$ similar to what we assumed for PSR
J2144-3933, which is conservative since the time-dependence can be larger for
larger values of $\alpha$. Throughout we assumed an NFW profile. The dark blue
exclusion region is from the CAST helioscope [70]. The green and purple
exclusions are from the ADMX [24, 25, 24, 26, 27, 28, 29, 30, 31] and HAYSTAC
[36, 37, 38] haloscopes, respectively. The yellow exclusions are from the CAPP
haloscope [33, 32, 112, 34, 35, 113] (Note that the high-frequency constraint
comes from the CAST-CAPP haloscope operating inside the CAST dipole magnet at
CERN). The orange are from the RBF/UF haloscopes [114, 115]. Finally, the dark
blue limits from the ORGAN collaboration [116, 117]. We have also included the
limits of the frequency only analysis of the pulsar population at the galactic
centre using data from the Breakthrough Listen project [64] as probably the
most reliable limit, using all the known physics [58, 59] from neutron stars
available prior to the analysis performed here.

In Fig. 13 we have placed the limits derived here in context of other limits
on $g_{a\gamma\gamma}$ from the CAST helioscope [70], haloscopes [24, 25, 24,
26, 27, 28, 29, 30, 31, 118, 37, 38, 33, 32, 112, 34, 35, 114, 115, 116, 117].
There are a number of other limits in the literature from neutron stars [68,
65, 75, 60, 64, 69]. We have included [64] as the present best limits obtained
from observations that only use frequency information, but have not included
the rest. For example, when it comes to GCM observations, the latest work by
three of the present authors [60], was the most up-to-date since that work
includes ray-tracing, however this was based on a previous version of ray
tracing-code which, unlike the present work, did not include multiple
reflections that increase the predicted signal, leading to overly-conservative
constraints. The GCM will be re-analysed in [79], combining the most up-to-
date account of the combined results on modelling from [59, 58] and data [65,
75, 67]. Ref. [69] also does not include ray-tracing so we do not include it.
Similarly [68] (which we again do not include) has been superseded by the
authors’ follow-up work [64] (shown in Fig. 13), which uses the most up-to-
date modelling. We have tried to be fair in displaying those results which use
the most up-to-date modelling and conservative assumptions.

In this work, we set out with the intention of examining to what extent
detailed time-domain information could be leveraged to increase the reach of
radio searches for axions relative to a simple radio-line search which simply
averages the flux over a long-time. We were able to quantify this precisely in
terms of the time-variance of the signal which we examined both for two
specific pulsars and a range of other pulsar magnetic field strengths. It
seems that for characteristic pulsar parameters there is a modest enhancement
to the signal to noise from including time-domain information, however this
marginal gain could be enough to tip a tentative detection in a total-flux
measurement into a signal-to-noise level above 5 when time-domain information
is included, making it well worthwhile to extract maximum leverage from pulsar
observations. This is especially relevant as we look to future telescopes such
as the SKA where we want to use all possible tools at our disposal to enhance
the prospects for detection. Furthermore, given the rich variety of ever-
increasing astrophysical probes of axions [119, 120, 121, 122, 123] events
which are sharply peaked in time (or other detailed features amenable to a
matched filter search) could benefit from more sophisticated search
strategies.

Acknowledgements

JIM thanks Sam Witte for useful discussions and is supported by an FSR
Fellowship. SS was supported formerly by a George Rigg Scholarship and more
recently by the UK Science and Technology Facilities Council (STFC). Pulsar
research at the Jodrell Bank Centre for Astrophysics is supported by a
consolidated grant from the STFC. The data used in this project was obtained
through the MeerTIME project. We would like to thank that collaboration for
making this data available to us.

##  Appendix A Interpolation of the ray-tracing results

Figure 14: Examples of the interpolated profiles in the analysis. Each of
these is for PSR J2144$-$3933 using the parameters discussed in the main text.
The top panel is for $(\alpha,\theta)=(10^{\circ},30^{\circ})$, the middle for
$(\alpha,\theta)=(60^{\circ},40^{\circ})$ and the bottom for
$(\alpha,\theta)=(80^{\circ},60^{\circ})$. It is clear that the interpolation
is not perfect, but it is accurate to within the level required for our
analysis.

Due to the computational cost of the ray-tracing simulations, that require
$\sim 24$ hours to produce a pulse-profile when parallelized over 32 CPU
cores, we require a faster alternative to predict the time-dependence of the
signal for arbitrary input angles. Therefore, we generate a simulated database
of flux profiles as a discrete function of $(\theta,\alpha)$, represented by
the circular data points in the top panels of Fig. 14 with
$\Delta\theta=10^{\circ}$ and $\Delta\alpha=5^{\circ}$. Based on these
datasets we generate an interpolation routine (where we use the SciPy package
scipy.interpolate) that can then predict the signal for arbitrary values of
$\alpha$ and $\theta$, the performance of which can be seen in the bottom
panels of Fig. 14, where we compare the prediction of our interpolation
routine with the data points. For the purpose of our analysis, the level of
agreement is sufficient.

## References

  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977).
  * Kim [1979] J. E. Kim, Physical Review Letters 43, 103 (1979).
  * Shifman  _et al._ [1980] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nuclear Physics B 166, 493 (1980).
  * Dine  _et al._ [1981] M. Dine, W. Fischler, and M. Srednicki, Physics Letters B 104, 199 (1981).
  * Zhitnitsky [1980] A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980), [Yad. Fiz.31,497(1980)].
  * Dine and Fischler [1983] M. Dine and W. Fischler, Physics Letters B 120, 137 (1983).
  * Abbott and Sikivie [1983] L. Abbott and P. Sikivie, Physics Letters B 120, 133 (1983).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Physics Letters B 120, 127 (1983).
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Svrcek and Witten [2006] P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hep-th/0605206 .
  * Davis and Shellard [1989] R. L. Davis and E. P. S. Shellard, Nucl. Phys. B 324, 167 (1989).
  * Battye and Shellard [1994] R. A. Battye and E. P. S. Shellard, Phys. Rev. Lett. 73, 2954 (1994), [Erratum: Phys.Rev.Lett. 76, 2203–2204 (1996)], arXiv:astro-ph/9403018 .
  * Battye and Shellard [1999] R. A. Battye and E. P. S. Shellard, ArXiv Astrophysics e-prints (1999), astro-ph/9909231 .
  * Wantz and Shellard [2010] O. Wantz and E. P. S. Shellard, Phys. Rev. D 82, 123508 (2010), arXiv:0910.1066 [astro-ph.CO] .
  * Bae  _et al._ [2008] K. J. Bae, J.-H. Huh, and J. E. Kim, JCAP 09, 005 (2008), arXiv:0806.0497 [hep-ph] .
  * Marsh [2016] D. J. E. Marsh, Physics Reports 643, 1 (2016), arXiv:1510.07633 .
  * Borsanyi  _et al._ [2016] S. Borsanyi _et al._ , Nature 539, 69 (2016), arXiv:1606.07494 [hep-lat] .
  * Choi  _et al._ [2021] K. Choi, S. H. Im, and C. Sub Shin, Ann. Rev. Nucl. Part. Sci. 71, 225 (2021), arXiv:2012.05029 [hep-ph] .
  * Planck Collaboration VI [2018] Planck Collaboration VI, ArXiv e-prints (2018), arXiv:1807.06209 .
  * Sikivie [1983] P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
  * De Panfilis  _et al._ [1987] S. De Panfilis, A. C. Melissinos, B. E. Moskowitz, J. T. Rogers, Y. K. Semertzidis, W. Wuensch, H. J. Halama, A. G. Prodell, W. B. Fowler, and F. A. Nezrick, Phys. Rev. Lett. 59, 839 (1987).
  * Wuensch  _et al._ [1989] W. Wuensch, S. De Panfilis-Wuensch, Y. K. Semertzidis, J. T. Rogers, A. C. Melissinos, H. J. Halama, B. E. Moskowitz, A. G. Prodell, W. B. Fowler, and F. A. Nezrick, Phys. Rev. D 40, 3153 (1989).
  * Hagmann  _et al._ [1990a] C. Hagmann, P. Sikivie, N. S. Sullivan, and D. B. Tanner, Phys. Rev. D 42, 1297 (1990a).
  * Asztalos  _et al._ [2010] S. J. Asztalos _et al._ (ADMX), Phys. Rev. Lett. 104, 041301 (2010), arXiv:0910.5914 [astro-ph.CO] .
  * Du  _et al._ [2018a] N. Du _et al._ (ADMX Collaboration), Phys. Rev. Lett. 120, 151301 (2018a).
  * Du  _et al._ [2018b] N. Du _et al._ (ADMX), Phys. Rev. Lett. 120, 151301 (2018b), arXiv:1804.05750 [hep-ex] .
  * Braine  _et al._ [2020] T. Braine _et al._ (ADMX), Phys. Rev. Lett. 124, 101303 (2020), arXiv:1910.08638 [hep-ex] .
  * Bartram  _et al._ [2021a] C. Bartram _et al._ (ADMX), Phys. Rev. Lett. 127, 261803 (2021a), arXiv:2110.06096 [hep-ex] .
  * Boutan  _et al._ [2018] C. Boutan _et al._ (ADMX), Phys. Rev. Lett. 121, 261302 (2018), arXiv:1901.00920 [hep-ex] .
  * Bartram  _et al._ [2021b] C. Bartram _et al._ (ADMX), (2021b), arXiv:2110.10262 [hep-ex] .
  * Crisosto  _et al._ [2020] N. Crisosto, P. Sikivie, N. S. Sullivan, D. B. Tanner, J. Yang, and G. Rybka, Phys. Rev. Lett. 124, 241101 (2020), arXiv:1911.05772 [astro-ph.CO] .
  * Lee  _et al._ [2020a] S. Lee, S. Ahn, J. Choi, B. R. Ko, and Y. K. Semertzidis, Phys. Rev. Lett. 124, 101802 (2020a), arXiv:2001.05102 [hep-ex] .
  * Jeong  _et al._ [2020] J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, J. E. Kim, and Y. K. Semertzidis, Physical Review Letters 125 (2020), 10.1103/physrevlett.125.221302.
  * Lee  _et al._ [2022] Y. Lee, B. Yang, H. Yoon, M. Ahn, H. Park, B. Min, D. Kim, and J. Yoo, Physical Review Letters 128 (2022), 10.1103/physrevlett.128.241805.
  * Kim  _et al._ [2022] J. Kim, O. Kwon, C. Kutlu, W. Chung, A. Matlashov, S. Uchaikin, A. F. van Loo, Y. Nakamura, S. Oh, H. Byun, D. Ahn, and Y. K. Semertzidis, “Near-quantum-noise axion dark matter search at capp,”  (2022).
  * Brubaker  _et al._ [2017] B. Brubaker, L. Zhong, S. Lamoreaux, K. Lehnert, and K. van Bibber, Physical Review D 96 (2017), 10.1103/physrevd.96.123008.
  * Zhong  _et al._ [2018] L. Zhong _et al._ (HAYSTAC), Phys. Rev. D 97, 092001 (2018), arXiv:1803.03690 [hep-ex] .
  * Backes  _et al._ [2021] K. M. Backes _et al._ (HAYSTAC), Nature 590, 238 (2021), arXiv:2008.01853 [quant-ph] .
  * Alesini  _et al._ [2019] D. Alesini _et al._ , Phys. Rev. D 99, 101101 (2019), arXiv:1903.06547 [physics.ins-det] .
  * Alesini  _et al._ [2021] D. Alesini _et al._ , Phys. Rev. D 103, 102004 (2021), arXiv:2012.09498 [hep-ex] .
  * McAllister  _et al._ [2017a] B. T. McAllister, G. Flower, J. Kruger, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, Phys. Dark Univ. 18, 67 (2017a), arXiv:1706.00209 [physics.ins-det] .
  * Quiskamp  _et al._ [2022a] A. P. Quiskamp, B. T. McAllister, P. Altin, E. N. Ivanov, M. Goryachev, and M. E. Tobar, Sci. Adv. 8, abq3765 (2022a), arXiv:2203.12152 [hep-ex] .
  * Melcón  _et al._ [2020] A. A. Melcón _et al._ (CAST), JHEP 21, 075 (2020), arXiv:2104.13798 [hep-ex] .
  * Chang  _et al._ [2022] H. Chang _et al._ (TASEH), Phys. Rev. Lett. 129, 111802 (2022), arXiv:2205.05574 [hep-ex] .
  * Grenet  _et al._ [2021] T. Grenet, R. Ballou, Q. Basto, K. Martineau, P. Perrier, P. Pugnat, J. Quevillon, N. Roch, and C. Smith, (2021), arXiv:2110.14406 [hep-ex] .
  * Millar  _et al._ [2022] A. J. Millar _et al._ , (2022), arXiv:2210.00017 [hep-ph] .
  * Liu  _et al._ [2022] J. Liu _et al._ (BREAD), Phys. Rev. Lett. 128, 131801 (2022), arXiv:2111.12103 [physics.ins-det] .
  * Beurthey  _et al._ [2020] S. Beurthey _et al._ , (2020), arXiv:2003.10894 [physics.ins-det] .
  * Brouwer  _et al._ [2022] L. Brouwer _et al._ (DMRadio), (2022), arXiv:2204.13781 [hep-ex] .
  * Marsh  _et al._ [2023] D. J. E. Marsh, J. I. McDonald, A. J. Millar, and J. Schütte-Engel, Phys. Rev. D 107, 035036 (2023), arXiv:2209.12909 [hep-ph] .
  * McDonald [2022] J. I. McDonald, Phys. Rev. D 105, 083010 (2022), arXiv:2108.13894 [hep-ph] .
  * Schütte-Engel  _et al._ [2021] J. Schütte-Engel, D. J. E. Marsh, A. J. Millar, A. Sekine, F. Chadha-Day, S. Hoof, M. N. Ali, K.-C. Fong, E. Hardy, and L. Šmejkal, JCAP 08, 066 (2021), arXiv:2102.05366 [hep-ph] .
  * Pshirkov and Popov [2009] M. S. Pshirkov and S. B. Popov, J. Exp. Theor. Phys. 108, 384 (2009), arXiv:0711.1264 [astro-ph] .
  * Huang  _et al._ [2018] F. P. Huang, K. Kadota, T. Sekiguchi, and H. Tashiro, Phys. Rev. D 97, 123001 (2018), arXiv:1803.08230 [hep-ph] .
  * Hook  _et al._ [2018] A. Hook, Y. Kahn, B. Safdi, and Z. Sun, Phys. Rev. Lett. 121, 241102 (2018).
  * Pshirkov [2009] M. S. Pshirkov, J. Exp. Theor. Phys. 108, 384 (2009).
  * Leroy  _et al._ [2020] M. Leroy, M. Chianese, T. D. P. Edwards, and C. Weniger, Phys. Rev. D 101, 123003 (2020), arXiv:1912.08815 [hep-ph] .
  * Battye  _et al._ [2021] R. A. Battye, B. Garbrecht, J. I. McDonald, and S. Srinivasan, JHEP 09, 105 (2021), arXiv:2104.08290 [hep-ph] .
  * Witte  _et al._ [2021] S. J. Witte, D. Noordhuis, T. D. P. Edwards, and C. Weniger, Phys. Rev. D 104, 103030 (2021), arXiv:2104.07670 [hep-ph] .
  * Battye  _et al._ [2020] R. Battye, B. Garbrecht, J. McDonald, F. Pace, and S. Srinivasan, Physical Review D 102 (2020), 10.1103/physrevd.102.023504.
  * Millar  _et al._ [2021] A. J. Millar, S. Baum, M. Lawson, and M. C. D. Marsh, Journal of Cosmology and Astro-Particle Physics 2021, 013 (2021), arXiv:2107.07399 [hep-ph] .
  * Goldreich and Julian [1969] P. Goldreich and W. H. Julian, ApJ 157, 869 (1969).
  * Safdi  _et al._ [2019] B. R. Safdi, Z. Sun, and A. Y. Chen, Physical Review D 99 (2019), 10.1103/PhysRevD.99.123021, arXiv:1811.01020 .
  * Foster  _et al._ [2022] J. W. Foster, S. J. Witte, M. Lawson, T. Linden, V. Gajjar, C. Weniger, and B. R. Safdi, arXiv e-prints , arXiv:2202.08274 (2022), arXiv:2202.08274 [astro-ph.CO] .
  * Darling [2020a] J. Darling, Phys. Rev. Lett. 125, 121103 (2020a), arXiv:2008.01877 [astro-ph.CO] .
  * Darling [2020b] J. Darling, Astrophys. J. Lett. 900, L28 (2020b), arXiv:2008.11188 [astro-ph.CO] .
  * Battye  _et al._ [2022] R. A. Battye, J. Darling, J. I. McDonald, and S. Srinivasan, Phys. Rev. D 105, L021305 (2022), arXiv:2107.01225 [astro-ph.CO] .
  * Foster  _et al._ [2020] J. W. Foster, Y. Kahn, O. Macias, Z. Sun, R. P. Eatough, V. I. Kondratiev, W. M. Peters, C. Weniger, and B. R. Safdi, Phys. Rev. Lett. 125, 171301 (2020), arXiv:2004.00011 [astro-ph.CO] .
  * Zhou  _et al._ [2022] Y.-F. Zhou _et al._ , Phys. Rev. D 106, 083006 (2022), arXiv:2209.09695 [astro-ph.CO] .
  * Anastassopoulos  _et al._ [2017] V. Anastassopoulos, S. Aune, Barth, _et al._ , Nature Physics 13, 584 (2017).
  * Raffelt and Stodolsky [1988] G. Raffelt and L. Stodolsky, Phys. Rev. D37, 1237 (1988).
  * Hardy and Song [2022] E. Hardy and N. Song, (2022), arXiv:2212.09756 [hep-ph] .
  * Huang  _et al._ [2018] F. P. Huang, K. Kadota, T. Sekiguchi, and H. Tashiro, Phys. Rev. D 97, 123001 (2018), arXiv:1803.08230 [hep-ph] .
  * Darling [2020a] J. Darling, Phys. Rev. Lett. 125, 121103 (2020a), arXiv:2008.01877 [astro-ph.CO] .
  * Darling [2020b] J. Darling, The Astrophysical Journal Letters 900, L28 (2020b), arXiv:2008.11188 [astro-ph.CO] .
  * Foster  _et al._ [2020] J. W. Foster, Y. Kahn, O. Macias, Z. Sun, R. P. Eatough, V. I. Kondratiev, W. M. Peters, C. Weniger, and B. R. Safdi, Phys. Rev. Lett. 125, 171301 (2020), arXiv:2004.00011 [astro-ph.CO] .
  * Leroy  _et al._ [2020] M. Leroy, M. Chianese, T. D. P. Edwards, and C. Weniger, Phys. Rev. D 101, 123003 (2020), arXiv:1912.08815 [hep-ph] .
  * Carenza and Marsh [2023] P. Carenza and M. C. D. Marsh, (2023), arXiv:2302.02700 [hep-ph] .
  * [79] J. McDonald and S. Witte, “In prepartion,” .
  * Gedalin and Melrose [2001] M. Gedalin and D. B. Melrose, Phys. Rev. E 64, 027401 (2001).
  * Turimov  _et al._ [2019] B. Turimov, B. Ahmedov, A. Abdujabbarov, and C. Bambi, Int. J. Mod. Phys. D 28, 2040013 (2019), arXiv:1802.03293 [gr-qc] .
  * Belvedere  _et al._ [2015] R. Belvedere, J. A. Rueda, and R. Ruffini, The Astrophysical Journal 799, 23 (2015).
  * Thompson and Duncan [1995] C. Thompson and R. C. Duncan, Monthly Notices of the Royal Astronomical Society 275, 255 (1995).
  * Thompson and Duncan [2001] C. Thompson and R. C. Duncan, Astrophys. J. 561, 980 (2001), arXiv:astro-ph/0110675 [astro-ph] .
  * Mondal [2021] T. Mondal, The Astrophysical Journal Letters 913, L12 (2021).
  * Pons and Perna [2011] J. A. Pons and R. Perna, The Astrophysical Journal 741, 123 (2011).
  * Platts  _et al._ [2020] E. Platts, J. X. Prochaska, and C. J. Law, The Astrophysical Journal Letters 895, L49 (2020), arXiv:2005.06256 [astro-ph.GA] .
  * Yamasaki and Totani [2020] S. Yamasaki and T. Totani, Astrophys. J. 888, 105 (2020), arXiv:1909.00849 [astro-ph.HE] .
  * Yao  _et al._ [2017] J. M. Yao, R. N. Manchester, and N. Wang, Astrophys. J. 835, 29 (2017), arXiv:1610.09448 [astro-ph.GA] .
  * Cordes and Lazio [2003] J. M. Cordes and T. J. W. Lazio, arXiv e-prints , astro-ph/0301598 (2003), arXiv:astro-ph/0301598 [astro-ph] .
  * Finkbeiner [2003] D. P. Finkbeiner, The Astrophysical Journal Supplementary Series 146, 407 (2003), arXiv:astro-ph/0301558 [astro-ph] .
  * Cordes and Lazio [2002] J. M. Cordes and T. J. W. Lazio, arXiv e-prints , astro-ph/0207156 (2002), arXiv:astro-ph/0207156 [astro-ph] .
  * Deller  _et al._ [2009] A. T. Deller, S. J. Tingay, M. Bailes, and J. E. Reynolds, Astrophys. J. 701, 1243 (2009), arXiv:0906.3897 [astro-ph.SR] .
  * Meyer and Christensen [2016] R. Meyer and N. Christensen, Significance 13, 20 (2016), https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1740-9713.2016.00896.x .
  * Allen [2005] B. Allen, Phys. Rev. D 71, 062001 (2005), arXiv:gr-qc/0405045 .
  * Manchester  _et al._ [2005] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astron. J. 129, 1993 (2005), arXiv:astro-ph/0412641 .
  * Johnston  _et al._ [2020] S. Johnston, A. Karastergiou, M. J. Keith, X. Song, P. Weltevrede, F. Abbate, M. Bailes, S. Buchner, F. Camilo, M. Geyer, B. Hugo, A. Jameson, M. Kramer, A. Parthasarathy, D. J. Reardon, A. Ridolfi, M. Serylak, R. M. Shannon, R. Spiewak, W. van Straten, V. Venkatraman Krishnan, F. Jankowski, B. W. Meyers, L. Oswald, B. Posselt, C. Sobey, A. Szary, and J. van Leeuwen, Monthly Notices of the Royal Astronomical Society 493, 3608 (2020).
  * Bailes  _et al._ [2020] M. Bailes, A. Jameson, F. Abbate, E. D. Barr, N. D. R. Bhat, L. Bondonneau, M. Burgay, S. J. Buchner, F. Camilo, D. J. Champion, I. Cognard, P. B. Demorest, P. C. C. Freire, T. Gautam, M. Geyer, J. M. Griessmeier, L. Guillemot, H. Hu, F. Jankowski, S. Johnston, A. Karastergiou, R. Karuppusamy, D. Kaur, M. J. Keith, M. Kramer, J. van Leeuwen, M. E. Lower, Y. Maan, M. A. McLaughlin, B. W. Meyers, S. Osłowski, L. S. Oswald, A. Parthasarathy, T. Pennucci, B. Posselt, A. Possenti, S. M. Ransom, D. J. Reardon, A. Ridolfi, C. T. G. Schollar, M. Serylak, G. Shaifullah, M. Shamohammadi, R. M. Shannon, C. Sobey, X. Song, R. Spiewak, I. H. Stairs, B. W. Stappers, W. van Straten, A. Szary, G. Theureau, V. Venkatraman Krishnan, P. Weltevrede, N. Wex, T. D. Abbott, G. B. Adams, J. P. Burger, R. R. G. Gamatham, M. Gouws, D. M. Horn, B. Hugo, A. F. Joubert, J. R. Manley, K. McAlpine, S. S. Passmoor, A. Peens-Hough, Z. R. Ramudzuli, A. Rust, S. Salie, L. C. Schwardt, R. Siebrits, G. Van Tonder, V. Van Tonder, and M. G. Welz, Publications of the Astronomical Society of Australia 37, e028 (2020).
  * Main  _et al._ [2022] R. A. Main, A. Parthasarathy, S. Johnston, A. Karastergiou, A. Basu, A. D. Cameron, M. J. Keith, L. S. Oswald, B. Posselt, D. J. Reardon, X. Song, and P. Weltevrede, Monthly Notices of the Royal Astronomical Society 518, 1086 (2022), https://academic.oup.com/mnras/article-pdf/518/1/1086/47139759/stac3149.pdf .
  * Johnston  _et al._ [2022] S. Johnston, M. Kramer, A. Karastergiou, M. Keith, L. Oswald, A. Parthasarathy, and P. Weltevrede, arXiv e-prints , arXiv:2212.03988 (2022), arXiv:2212.03988 [astro-ph.HE] .
  * Lorimer and Kramer [2012] D. R. Lorimer and M. Kramer,  _Handbook of Pulsar Astronomy_ (2012).
  * Gil  _et al._ [1984] J. Gil, P. Gronkowski, and W. Rudnicki, Astronomy and Astrophysics 132, 312 (1984).
  * Posselt  _et al._ [2021] B. Posselt, A. Karastergiou, S. Johnston, A. Parthasarathy, M. J. Keith, L. S. Oswald, X. Song, P. Weltevrede, E. D. Barr, S. Buchner, M. Geyer, M. Kramer, D. J. Reardon, M. Serylak, R. M. Shannon, R. Spiewak, and V. Venkatraman Krishnan, Monthly Notices of the Royal Astronomical Society 508, 4249 (2021), arXiv:2109.12007 [astro-ph.HE] .
  * Rookyard  _et al._ [2015] S. C. Rookyard, P. Weltevrede, and S. Johnston, Monthly Notices of the Royal Astronomical Society 446, 3367 (2015), arXiv:1410.3294 [astro-ph.HE] .
  * Johnston and Karastergiou [2019] S. Johnston and A. Karastergiou, Monthly Notices of the Royal Astronomical Society 485, 640 (2019), arXiv:1902.03996 [astro-ph.HE] .
  * Weltevrede and Johnston [2008] P. Weltevrede and S. Johnston, Monthly Notices of the Royal Astronomical Society 391, 1210 (2008), https://academic.oup.com/mnras/article-pdf/391/3/1210/17319973/mnras391-1210.pdf .
  * Radhakrishnan and Cooke [1969] V. Radhakrishnan and D. J. Cooke, Astrophysical Letters 3, 225 (1969).
  * Yuen and Melrose [2014] R. Yuen and D. B. Melrose, Publications of the Astronomical Society of Australia 31, e039 (2014).
  * Dyks and Pierbattista [2015] J. Dyks and M. Pierbattista, Monthly Notices of the Royal Astronomical Society 454, 2216 (2015), arXiv:1507.08454 [astro-ph.HE] .
  * Cordes [1978] J. M. Cordes, Astrophys. J. 222, 1006 (1978).
  * Spitler  _et al._ [2014] L. G. Spitler, K. J. Lee, R. P. Eatough, M. Kramer, R. Karuppusamy, C. G. Bassa, I. Cognard, G. Desvignes, A. G. Lyne, B. W. Stappers, G. C. Bower, J. M. Cordes, D. J. Champion, and H. Falcke, The Astrophysical Journal Letters 780, L3 (2014), arXiv:1309.4673 [astro-ph.HE] .
  * Lee  _et al._ [2020b] S. Lee, S. Ahn, J. Choi, B. Ko, and Y. Semertzidis, Physical Review Letters 124 (2020b), 10.1103/physrevlett.124.101802.
  * Adair  _et al._ [2022] C. M. Adair, K. Altenmüller, V. Anastassopoulos, S. A. Cuendis, J. Baier, K. Barth, A. Belov, D. Bozicevic, H. Bräuninger, G. Cantatore, F. Caspers, J. F. Castel, S. A. Çetin, W. Chung, H. Choi, J. Choi, T. Dafni, M. Davenport, A. Dermenev, K. Desch, B. Döbrich, H. Fischer, W. Funk, J. Galan, A. Gardikiotis, S. Gninenko, J. Golm, M. D. Hasinoff, D. H. H. Hoffmann, D. D. Ibáñez, I. G. Irastorza, K. Jakovčić, J. Kaminski, M. Karuza, C. Krieger, Ç. Kutlu, B. Lakić, J. M. Laurent, J. Lee, S. Lee, G. Luzón, C. Malbrunot, C. Margalejo, M. Maroudas, L. Miceli, H. Mirallas, L. Obis, A. Özbey, K. Özbozduman, M. J. Pivovaroff, M. Rosu, J. Ruz, E. Ruiz-Chóliz, S. Schmidt, M. Schumann, Y. K. Semertzidis, S. K. Solanki, L. Stewart, I. Tsagris, T. Vafeiadis, J. K. Vogel, M. Vretenar, S. Youn, and K. Zioutas, Nature Communications 13 (2022), 10.1038/s41467-022-33913-6.
  * DePanfilis  _et al._ [1987] S. DePanfilis, A. C. Melissinos, B. E. Moskowitz, J. T. Rogers, Y. K. Semertzidis, W. U. Wuensch, H. J. Halama, A. G. Prodell, W. B. Fowler, and F. A. Nezrick, Phys. Rev. Lett. 59, 839 (1987).
  * Hagmann  _et al._ [1990b] C. Hagmann, P. Sikivie, N. S. Sullivan, and D. B. Tanner, Phys. Rev. D 42, 1297 (1990b).
  * McAllister  _et al._ [2017b] B. T. McAllister, G. Flower, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, Physics of the Dark Universe 18, 67 (2017b).
  * Quiskamp  _et al._ [2022b] A. Quiskamp, B. T. McAllister, P. Altin, E. N. Ivanov, M. Goryachev, and M. E. Tobar, Science Advances 8 (2022b), 10.1126/sciadv.abq3765.
  * Brubaker and et al. [2017] B. M. Brubaker and et al., Phys. Rev. Lett. 118, 061302 (2017).
  * Witte  _et al._ [2022] S. J. Witte, S. Baum, M. Lawson, M. C. D. Marsh, A. J. Millar, and G. Salinas, (2022), arXiv:2212.08079 [hep-ph] .
  * Noordhuis  _et al._ [2022] D. Noordhuis, A. Prabhu, S. J. Witte, A. Y. Chen, F. Cruz, and C. Weniger, (2022), arXiv:2209.09917 [hep-ph] .
  * Chadha-Day  _et al._ [2022] F. Chadha-Day, B. Garbrecht, and J. McDonald, JCAP 12, 008 (2022), arXiv:2207.07662 [hep-ph] .
  * Escudero  _et al._ [2023] M. Escudero, C. K. Pooni, M. Fairbairn, D. Blas, X. Du, and D. J. E. Marsh, (2023), arXiv:2302.10206 [hep-ph] .
  * Day and McDonald [2019] F. V. Day and J. I. McDonald, JCAP 10, 051 (2019), arXiv:1904.08341 [hep-ph] .