Searching for Time-Dependent Axion Dark Matter Signals in Pulsars

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 - 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}<4\times 10^{-11}\,{\rm GeV}^{-1}$ 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.

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 a ≲ µ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 ∼ 10 5 G (10 T), a limiting factor in such searches.By contrast, astrophysical magnetic fields in neutron stars can be as high as ∼ 10 15 G (10 11 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 ≃ ω p , where ω 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.[54][55][56] 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 aγγ , 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γγ and by injecting signal templates into radio data, demonstrate a practical route to obtaining limits on g aγγ 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 timedependence 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 timedomain 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 photons 1 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 [64,67,[74][75][76] 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 [58,59,61,77,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 [58,59,77,79].

A. Ray-tracing photons in magnetised plasmas
Resonant conversion between axions and photons occurs at points where k µ γ = k µ a where k µ γ and k µ a are the photon and axion 4-momentum, respectively.An attempt was made to understand the conversion probability for axions to photons p aγ in [61] leading to 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] .
1 See [72] for a similar mechanism with dark photons.
Here, the sum s is over different charge carrier species, γ s is a generalised Lorentz factor, and k ∥ gives the 4momentum projected onto the magnetic field and v s corresponds to the velocity of charge carriers and µ 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, γ = 1, c s = 0 and µ 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 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.

B. 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] with the plasma frequency given by ω p = 4π|n GJ |/m e .
Here, Ω is the frequency of the pulsar, α 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 θ and azimuthal angle ϕ are defined with respect to the rotational axis of the star.
It is obvious that whenever α ̸ = 0, the plasma is timedependent 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 ρ 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 Ṗ 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 Ṗ 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 ∼ 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 √ S ∝ B 0.4 (see sec.V) and s such, we leave a more sophisticated modelling of the magnetosphere for future work., but it can be substantial for intermediate angles.When θ = 90 • 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.
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 ≈ 0.29 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 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 ∝ 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 (α, θ).We, therefore, generate a simulated database 2 of periodic flux profiles as a function of (α, θ).Some of these profiles are displayed in Fig. 1 which indicates the (α, θ) dependence of the time-variability of the signal.
2 Formally we generate templates for discrete (α i , θ i ) and numerically interpolate to generate a template for a continuous range of α and θ.See appendix A for a description of this procedure and an illustration of its accuracy.
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.

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, q, defined below, to quantify the likelihood of detection.In this sense, q acts as a statistical test for whether the data is distinguishable from noise.Values of 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 exp =< q > (see eq. ( 15)) from an axion signal.By comparing this to the measured value q, we can exclude regions of axion parameter space.
This procedure provides a means to derive limits on g aγγ as a function of m a , and importantly allows us to take into account that for a fixed value of m 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 α between the magnetic axis and spin axis of the star, and the angle θ between the line of sight and the spin axis.As in previous attempts to derive constraints on g aγγ using neutron stars [64][65][66][67][68] 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 α and θ.For each value of m a , we therefore obtain a constraint on the value of g aγγ for every pair (α,θ).We can then exclude certain ranges of α and θ with further modelling and observations of the pulsar signal, notably the pulse width.We then use the value of (α, θ) from this remaining subset which gives the most conservative constraints on g aγγ .

A. 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 B we use it to quantify the pros and cons of time domain observations and in section III C 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", d.We will assume the data is the sum of a signal S = S 0 F(p) and some additive noise n Here, we have decomposed the signal according to so that S 0 , gives the root mean squared flux density of the signal The signal is further characterised by a "parameter vector" p.In section II we will compute a number of templates F(p) for the signal where p = (m a , α, θ).The noise vector n is assumed to be Gaussian with ⟨n⟩ = 0 and ⟨nn T ⟩ = C where ⟨..⟩ denotes an ensemble average of noise realizations and C ij is the covariance matrix with i, j = 1, .., n d where n d is the total number of data points.Assuming a Gaussian likelihood, −2 log L = χ 2 , one can calculate the maximum likelihood estimate Ŝ0 by minimizing In the above equation, we assume the data vector d contains the true signal S = S true so that minimising χ 2 above can be thought of as minimising a generalised leastsquares difference (relative to the noise in each channelhence the factor C −1 ) between the true signal and possible templates S 0 F. Viewing χ 2 as an unknown function of S 0 , we can find the minimising value, Ŝ0 given by Ŝ0 = One can also deduce the "matched filter noise" and the signal-to-noise estimate q may then be written as It is important to understand the difference between the noise in the data, characterised by C, and the "matched filter noise", σ MF .They are related, but σ 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 where δ ij is the Kronecker delta and σ 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 d = n 1 × ... × 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 d = n f × n t where n f is the number of frequency channels and n t is the number of time samples.In that case, the data vector could be written as where d ωi tj labels the data vector in the ith frequency bin and jth 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 that is, the dot product of the filter F with the data vector.When F matches that in the true signal, we have ⟨q⟩ = S 0 /σ N (from eq. ( 5)).Now assume that p true has entries which are the true parameters.The ensemble average of q for a filter with arbitrary parameter, p is Next, we note that as a trivial consequence of the Cauchy-Schwarz inequality, we have Assuming that the template is non-degenerate with respect to the values of p, this occurs only for p = p true for which q is then maximal.Thus 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 nonvanishing entries in only one frequency channel where ω = m a .This means filters with different values of m 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 α and θ are not orthogonal and will in general have overlap such that 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 (α, θ), thereby obtaining the most conservative constraints on g aγγ for allowed values of angles.
In order to turn our continuous axion signals into discrete vectors we must perform some kind of coarsegraining.We therefore define a binning scheme for N time-channels centered on the points t i = (i − 1/2)∆t where i = 1, .., N and ∆t = 1/N so that the discretized signal is given by where S(t, ω) is the flux density of the axion signal as a function of time (and frequency) that we derive using our radio signal models.

B. 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 timedomain 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 Then according to Eqs ( 5)-( 7) and ( 14), the matched filter will return where σ N is the noise on each of the n t = N timechannels.The noise amplitude σN averaged across all times is then given by σ = σ N / √ N .We can therefore re-write (18) as where µ S = Σ i S i /N and σ S = Σ i S 2 i /N − µ 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 σ given by averaging the noise over all integration time, the signal √ S • S is simply given by the mean µ S so that 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 σ 2 S ≥ 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 (σ S /µ S ≫ 1), the timedomain search provides a gain in sensitivity by increasing the signal to noise by a factor ≃ σ S /µ S , i.e. the squareroot of the relative variance.This is a simple consequence of the fact that 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 subtracted3 .Thus while we retain time-domain information, the baseline µ S of the signal is essentially renormalised to zero.In this case, the time-domain analysis gives a baseline subtracted (BS) value Clearly, when there is only a small time variation, (µ S ≫ σ S ), subtracting the baseline leads to a lower value of ⟨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, µ S ≪ σ 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.
C. Demonstration on simulated data In the previous subsections we have derived the main properties of the matched filter estimate of the signal-tonoise 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 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 A.
We will inject signals with a mass of m a = 4.2 µeV which corresponds to an observing frequency f obs ≃ 1.0GHz and g aγγ = 10 −10 GeV −1 .We will consider two different choices for the angles α and θ.Case A with (α, θ) = (0 • , 60 • ) is close to an aligned rotator and hence we would expect no time variation, whereas case B with (α, θ) = (40 • , 60 • ) has a very strong time variation.We will consider two observations: one in which the pulseaveraged 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 B 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 (α = 0 • ), 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 timedependence (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 α = 0 case, the filter successfully returns the maximal SNR for the input value of θ.

FIG.
2. An illustration of using the matched filter SNR to search for the signal.q is shown as function of observing angle θ 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 (α, θ) = (0 • , 60 • ) (case A, top panel) and another with (α, θ) = (40 • , 60 • ) (case B, bottom panel).In both cases we choose a noise amplitude of ∼ 2.5 mJy consistent with that observed in our data for ma = 4.2 µeV.We use an axion mass of 3.7 µeV and gaγγ = 10 −10 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] 4 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 where ρ 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 ρ DM ≈ 0.45 GeV 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−2900 5 , this is likely to give a reasonable estimate of the trade-off between ρ DM and D in the FOM.
The PSR J2144−3933 which has B 0 ≈ 2.1 × 10 12 G estimated from the P and Ṗ based on electromagnetic spin down, P = 8.51 sec and D = 0.12 kpc came third on the list 6 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 ρ a used in our predictions.
Within the GJ model there is a maximum axion mass [67] given by which is ≈ 4.7 µeV corresponding to f obs ≈ 1.15 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 ≳ 512 pulses from each pulsar, and hence the total observing duration was 4416 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.

A. Modelling the noise
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, µ S , and the standard deviation, σ 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 ∼ 0.95 GHz, and Global Navigation Satellite System (GNSS), f ∼ 1.2 GHz and f ∼ 1.6 GHz7 .
This residual RFI can be removed by excising any data with σ N > 3 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 σ 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 A. 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 ϕ.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, ϕ.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.

B. Constraining the magnetic orientation α and observing angle θ
We have already pointed out that the amplitude of the time dependence of the signal depends on the values of α and θ -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γγ .
Obtaining precise values for the pulsar geometry re- lies 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 hereit 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 θ ≈ α ≈ 0.
Let W be the pulse width corresponding to a fully illuminated circular radiation beam with half-opening angle ρ.These parameters can be related to the parameters in our misaligned rotator model for the neutron star (α, θ) using [101,102] cos ρ = cos α cos θ + sin α sin θ cos(W/2) . ( The width of the profile for PSR J2144−3933 is measured at 10% of the amplitude is (2.1±0.2) • [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 ( 24) 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 • -5.0 • for this object.We now turn to the estimation of ρ, whose uncertainty mainly stems from a lack of knowledge of the height h 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]) where we have used the small angle approximation for ρ.
In the second expression we have replaced the light cylinder radius R c with the pulse period P = 2πR c /c. Estimation of h em is complicated by the fact that the beam is not necessarily filled, but across the pulsar population h em at 1.4 GHz has been constrained to be in the range of 200-400 km irrespective of pulse period [105].We take a more conservative range of 100 ≤ h em ≤ 1000 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 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 α 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 α − θ degeneracy.
Since parameters for a given pulsar are uncertain, values of α all the way down to zero are allowed by (24), for which there would be no time-variation.We, therefore, appeal to further arguments which allow us to place a lower bound on α by excluding implausible geometries.
A problem with very small α 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 α geometries are unlikely, we construct an effective probability distribution for α and β which essentially measures the number of beam realisations associated to each pair (α, β) 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 • and 5.0 • in 100 steps.(ii) For each given W, scan over a discrete grid of (α, β) values between 0 and π/2.(iii) For each point (α, β) calculate ρ(W, α, β) with Eq. (24).For each (α, β) (for the specific W under consideration) if the resulting ρ satisfies 100 ≤ h em ≤ 1000 km and |β/ρ| < 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 β ≃ ρ, which is both unlikely and implausible for the reasons above.At the end of this process, for each W , one has an α-β 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 |β| needs to be small in order for the line of sight to intersect the beam and small α solutions are excluded.No solutions exist for α ≲ 10 • .Furthermore, one can see in the top panel that solutions α ≲ 20 • are unlikely geometrical solutions, which is because they require a fine tuned (large) β.It should be stressed that α ≲ 20 • geometries are not just unlikely, we have also argued them to be contrived 8 .In what follows we will impose α > 20 • and |β| < 4 • .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γγ .

C. Constraints on gaγγ 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 a , θ, α) 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 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 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 ≈ 1.45 somewhat higher than the expected value which points to the fact that the noise model we are using is not optimal 9 .We find that, out of the ≈ 3 × 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 ≈ 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 meas., appear to be randomly distributed with m 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 α = 0, but can be larger than one for a range of values of (α, θ), but that we have argued that regions with α < 20 • are unlikely and contrived in section IV B.
In the absence of a detection, one can derive constraints on g aγγ by comparing the measured value of q meas.for a particular template, with the expected value ⟨q⟩ for that same template.The measured value is defined by This must be compared to the expected value ⟨q(g aγγ , θ, α)⟩ if a signal were present in the data In the perturbative limit of the conversion probability, which we have checked is always the case here, this is ∝ g aγγ .Therefore, in order to impose a limit we calculate ⟨q(g fid αγγ , θ, α)⟩ for g fid aγγ = 10 −10 GeV −1 and exclude any value such that The factor of two in (28) corresponds to the observed noise level being twice the equivalent signal for g aγγ meaning that this will be a 2 (≈ 95% confidence) upper limit on coupling constant.We have explained above that the typical values of q 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 (28) are not optimal and hence are slightly conservative.
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 show the constraints from this procedure in the left panel of Fig. 6 for m a = 4 µ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 (α, θ).We do this by taking the ratio of theoretical signal-to-noise in the case where the baseline is subtracted, q time , with respect to the case where the time profile is averaged over the pulse-period and integrated, q freq. .When this ratio is larger than 1, the time-domain analysis leads to stronger constraints on g aγγ .
We have performed the same analysis as presented in Fig. 6 over the mass range 3.9 µeV ≤ m a ≤ 4.7 µeV and then searched for the weakest upper limit in the range α > 20 • and |β| < 4 • .The results, assuming that D = 0.165 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γγ > 5.5 × 10 −11 GeV −1 .If D = 0.29 kpc as suggested by electron density models, we find a weaker constraint of g aγγ > 9.6 × 10 −11 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.

V. FUTURE SEARCHES
In the previous section we have shown how one can derive a limit on g aγγ 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 Hence, any improvement on g lim aγγ will come from two avenues.The first is better observations: lower system temperature T sys , larger collecting area A eff and increased observing time, t obs .The second is a better target: one with larger mean-signal power µ S or greater time variability as measured by σ 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. ( 29) implies that an increase in observation time from ≈ 1 hour, as we have now, to 100 hours would improve the limit by a factor ∼ 3.Moreover, going from MeerKAT with A eff /T sys ≈ 450 to A eff /T sys ≈ 1800 for SKA1-mid will yield a limit which is a factor ∼ 2 better.Combining these together one might be able to obtain a limit of g aγγ < 1.6 × 10 −11 GeV −1 .PSR J2144−3933 only allows constraints to be imposed for m a < 4.7 µ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 max a 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 ∝ B 2/3 0 P 7/3 , while the maximum mass probed is ∝ (B 0 /P ) 1/210 .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 max a 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 raytracing code to calculate the average signal, µ 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 θ and α, 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 α = 60 • we find that FOM ∝ B 0.8 0 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γγ 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 A and as is apparent from Eq. ( 29), the key parameter determining whether including a time-domain analysis can add value over just using the frequency domain is σ S /µ S .Clearly this will depend on the intrinsic pulsar parameters, B 0 , P and axion mass, m a as well as (α, θ).In Fig. 8 we show that σ S /µ S clearly increases with B 0 with the amount being sensitive to the choice of θ and α.We have also investigated the dependence on P , but this is typically much weaker for the relevant range of parameters.FIG. 6.In the left panel, we present the constraints on the axion-photon coupling for the average-subtracted case as a function of (θ, α) for ma = 4 µeV.As expected, the derived limit is weaker in the limit α → 0 since the time-dependence of the signal is negligibly small.In fact there is no limit for α ≡ 0 since there is no time dependence.Fortunately, the limits from the pulsar main beam modelling require α > 20 • and |β| < 4 • which are included as red lines on the plot, so we are able to achieve a limit gγγγ ≲ 9.6 × 10 −11 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.FIG. 7. 2 upper limits on gaγγ as a function of ma for d = 0.165 pc.This is determined by calculating the highest upper limit in the region of the α − θ plane allowed by |β| < 4 • and α > 20 • from the equivalent of Fig. 6.On the basis of this figure we quote an upper limit of gaγγ < 5.5 × 10 −11 GeV over the mass range 3.9 µeV ≤ ma ≤ 4.7 µeV.
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 FIG. 8.The relative time variance, σS/µS, of the profiles as a function of the magnetic field B0 at the surface of the neutron star.We fix P = 4 s and ma = 1 µeV and the value of gaγγ scales out in this ratio.Despite the time-variation between the maximum and minimum of the profiles increasing by orders of magnitude for B0 ∼ 10 14 G compared to B0 ∼ 10 12 G, σS/µS only goes up by a factor of a few.
in the GC, although this is mitigated somewhat by it being further away.
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.FIG. 9. Predicted pulse profiles for the GCM.In order to generate these profiles, we fix gaγγ = 10 −10 GeV −1 and ma = 1 µeV and α = 60 • 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 θ = 30 • and 60 • .Nonetheless these only correspond to a σS/µs ≈ 3 compatible with Fig. 8.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.
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 ≈ 1.4×10 14 G and a period of P = 3.76 sec and is a distance of D ≈ 8.3 kpc.In fact, the ratio of the FOM Eq. ( 22) from PSR J2144−3933 relative to the GCM is given by If one assumes a standard NFW profile for the galaxy, one obtains an enhancement of ∼ 10 5 with respect to the local value 11 which suggests that it will lead to a larger FOM by a factor of ∼ 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 σ S /µ 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 mJy, and this is presented in the top panel of Fig. 10 12 .
Making similar assumptions about constraints on (α, β) from the pulse profile, we obtain a projected upper limit of g aγγ < 4×10 −12 GeV −1 which is a factor ∼ 10 stronger than that we obtained from PSR J2144−3933.In addition, the GCM allows a much wider range masses with m max a ≈ 85µeV (f max obs ≈ 20 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 ∼ 6 using a 100 hour observation using SKA1-mid, which would lead to a projected limit of g aγγ < 6 × 10 −13 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 G ≤ B 0 ≤ 10 13 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 ρ DM , which translates to an uncertainty of ∼ 2 orders of magnitude in g aγγ 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.
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 (α, θ).We see that the values of q time /q freq larger, typically ≈ 2 − 3 compared to ≈ 1 for PSR J2144−3933.It seems reasonable to conclude that the factor ≈ 20 improvement in the constraining power seen in the case of the GCM comes from the combination of the FOM based on radial 12 We 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 Tsys.This means to achieve this noise figure one would need to observe for ∼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.
FIG. 10.In the top panel, we show the expected constraints on gaγγ from simulated observations of the GCM with an r.m.s.noise level of 3 mJy for ma = 1 µeV.We fix the pulsar parameters to be B0 = 1.4 × 10 14 G, P = 3.76 s and ρDM = 5.4×10 4 GeV cm −3 which is the value computed using a standard NFW profile for the galaxy.Note that, depending on the values of θ and α limits as low as gaγγ ≲ 10 −13 , GeV −1 might be possible.In the bottom panel, we quantify the effect of adding time-domain information over the parameter space (α, θ).Note that this has a slightly different morphology to that for PSR J2144−3933, but also the values of qtime/q freq are slightly larger approaching ≈ 3 as it is indicated might be the case in Fig. 8.
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γγ from PSR J2144−3933 are strongly dependent on constraints we have placed on α and β from the pulse profile.These are conservative and given the nature of the signal -there are regions of the α − θ plane where the constraints are very weak and even non-existent -meaning one is always dominated by the lower limit one imposes on α.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γγ ≲ 10 −13 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γγ 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γγ , 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 where P is the signal period.The coefficients a k are then given by 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 Note that in what follows we do not consider the k = 0 mode associated to the time-average µ 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 Im(a k ) = 0.In order to confirm our assumption that the profiles are dominated by low k modes we have computed the power spectrum 13 , |a k | 2 , of two profiles with varying levels of time-dependence, i.e., a relatively flat profile with (α = 10 • , θ = 10 • ) and one with a relatively large time variation with (α = 60 • , θ = 30 • ).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 cut .We have done this for n 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 cut = 3 and n 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.
The advantage of this approach is that we know a periodic signal search can be carried out using only the first n cut Fourier modes, reducing computational overhead when scanning.
Furthermore, in what follows, rather than scanning over every available template parametrised by (a 1 , • • • , a ncut ), 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 where Γ 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 scaleinvariant.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 con-  ) (red).The power spectrum has been normalised such that |a1| = 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 ncut = (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.FIG.12.In the top panel we present a histogram of the power spectrum computed from the data for ma = 4 µ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 χ 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 ma as in the top panel and have also included the theoretical PDFs which are χ 2 3 and χ 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.
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 cut modes, given by We then want to understand if the measured values of low-mode power are again consistent with the scaleinvariant spectrum.The PDF for Q(n cut ) is χ 2 ncut since it is the sum of n cut independent modes each ∼ χ 2 1 .This distribution can clearly be seen in the bottom two panels of Fig. 12 where for n cut = 3 and n cut = 7 we display simulated and analytic PDFs for the frequency channel corresponding to m a = 4 µeV.Note we have also included the measured value of Q(n 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 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 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 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 timedependent signals and have used this to show that timedomain 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 period 14 ) 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 θ between the stars rotation axis and the line of sight towards the pulsar, and α the angle between the stars magnetic axis and its rotation axis.Using the morphology of the observed pulsar mainbeam signal, we were able to exclude a range of values (α, θ), 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γγ < 9.6 × 10 −11 GeV −1 (D/0.29 pc) over the mass range 3.9 µeV ≤ m a ≤ 4.7 µ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.
There are a number of other limits in the literature from neutron stars [60,64,65,68,69,75].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 [58,59] and data [65,67,75].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 everincreasing 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.

FIG. 1 .
FIG. 1. Time-dependence of radio templates from axion dark matter as a function of pulse phase ϕ for different values of θ with fixed value of α = 20 • with gaγγ = 10 −10 GeV −1 .As expected there is little time variation for θ = 10• , but it can be substantial for intermediate angles.When θ = 90 • 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.

FIG. 3 .
FIG. 3. Measurements of the mean, µS and σ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 σS and µS for the full data-set.In the bottom panel, we show the equivalent after employing a cut of σN = 3.0 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 µ ≈ 0 and σN ≈ 2.8 mJy.

FIG. 4 .
FIG. 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 ϕ = 0 and ϕ = 1) removed.The top panel is our noise model described in the text, while the bottom panel is the actual data.

FIG. 5 .
FIG. 5.The constraint on the α − β plane based on the geometry of the beam where β = θ − α.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 hem we have allowed.Improbable solutions with |β/ρ| > 0.95 are rejected.The top panel shows the vertical integration of the bottom panel and it demonstrates that solutions with α < 20 • are very unlikely because they would require fine tuning of β.Therefore, we rule out α < 20 • and |β| > 4 • as statistically unlikely, indicated by the red-dashed lines, when we calculate our limit on gaγγ.
FIG. 11.In the top panel we present the power spectrum of two profiles with with (α = 60 • , θ = 30 • ) (blue) and (α = 10 • , θ = 10 • ) (red).The power spectrum has been normalised such that |a1| = 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 ncut = (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.