{'fraction_non_alphanumeric': 0.03576881134133043, 'fraction_numerical': 0.020837178089086486, 'mean_word_length': 4.381545684362586, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Searching for dark-matter waves with PPTA and QUIJOTE pulsar polarimetry

The polarization of photons emitted by astrophysical sources might be altered as they travel through a dark matter medium composed of ultra light axion-like particles (ALPs). In particular, the coherent oscillations of the ALP background in the galactic halo induce a periodic change on the polarization of the electromagnetic radiation emitted by local sources such as pulsars. Building up on previous works, we develop a new, more robust, analysis based on the generalised Lomb-Scargle periodogram to search for this periodic signal in the emission of the Crab supernova remnant observed by the QUIJOTE MFI instrument and 20 galactic pulsars from the Parkes Pulsar Timing Array (PPTA) project. We also carefully take into account the stochastic nature of the axion field, an effect often overlooked in previous works. This refined analysis leads to the strongest limits on the axion-photon coupling for a wide range of dark matter masses spanning $10^{-23}\text{ eV}\lesssim m_a\lesssim10^{-19} \text{ eV}$. Finally, we survey possible optimal targets and the potential sensitivity to axionic dark-matter in this mass range that could be achieved using pulsar polarimetry in the future.


Introduction
The fundamental nature of dark matter (DM) remains one of the major mysteries of modern physics [1]. One of the most promising candidates of particle DM is the axion [2][3][4], whose first incarnations were predicted as a consequence of the Peccei-Quinn solution to the strong CP -problem [5][6][7]. Equally interesting is the larger framework of axion-like particles (ALPs), pseudoscalar particles with Standard Model (SM) couplings resembling those of the QCD axion 1 , but not necessarily related to the strong CP -problem. As a consequence, their mass and decaying constant can be decoupled, widely broadening the possible phenomenological consequences.
A characteristic feature of ALPs is that they can be extremely light (∼ 10 −22 eV), becoming the prototypical example of bosonic Ultra Light Dark Matter (ULDM), also called fuzzy dark matter or wave dark matter [8][9][10][11][12][13]. Remarkably, this type of ultra light particles is naturally predicted in some extensions of the SM, including string theory [14,15]. These dark matter candidates are characterized by a large de Broglie wavelength, which suppresses structure at small scales. This observation has been used as a further motivation for these models, as this suppression may address some difficulties of the cold dark matter paradigm at galactic scales [16]. A prominent example is the core-cusp problem, arising from the discrepancy between the inferred dark matter density profiles of low-mass galaxies and the density profiles predicted by cosmological N-body simulations [17]. Other related observables where ULDM may play a critical role include: the Lyman-α forest [18,19], the abundance of high-redshift objects [20], galactic rotation curves [21,22] and the formation features of the Milky Way disk [23].
Another remarkable property of these DM candidates is that they can be easily (even "naturally") produced at the observed cosmological levels via the misalignment mechanism [2,3,24]. This relies on assuming that the field strength has a "natural" amplitude in the early universe (provided by, e.g., a phase transition) and that after some time (depending on the mass of the DM particle) they start to oscillate and behave as a cold dark matter fluid.
In this work we will focus on the ULDM parameter space of ALPs. One of the most interesting handles in their phenomenology is the non-renormalizable interaction of ALPs with electromagnetic (EM) fields. This interaction might be studied in several kinds of experiments such as CAST [25], or astrophysical phenomena like the supernova SN 1987A [26]. It also affects the emission, propagation and detection of EM waves from astrophysical objects (see e.g. Refs. [27][28][29] for reviews). In this work, we will describe how the polarization properties of light are modified in its propagation through an external ALP (pseudoscalar) field. In particular, we focus on the modulation of the light polarization angle induced by ALPs in virial equilibrium in galactic structures [13,30,31]. This effect has been searched for in, for instance, pulsars [32], parsec-scale jets in active galaxies [33], protoplanetary disks [34] or in observations of the CMB [35,36]. In the present manuscript we improve and refine the analysis of Ref. [32] using pulsars as light sources. Our analysis is novel in several ways: • First, we use generalized Lomb-Scargle periodograms. These are standard and powerful techniques used in astronomy to search for periodic signals in non-homogeneously distributed time series. This allows us to perform a robust analysis of the time dependence of the pulsar's polarization measurements constraining ALP-birefringence effects.
(Similar techniques were applied to active galaxy sources in this context [33]).
• Second, we use a more accurate model of the ULDM configuration. This takes into account the intrinsic stochastic nature of the amplitude of the field, which is relevant for the mass range of interest and experimental setup. This is tackled in our work by a global analysis of the polarization measurements of 20 pulsars observed for∼ 4.5 years by the Parkes Pulsar Time Array (PPTA) [37]. Furthermore, we include data from the polarization studies of the Crab Pulsar (and supernova remnant) that span also ∼ 4.5 years, performed by the QUIJOTE MFI instrument [38,39] to calibrate their measurements.
• Finally, we investigate the prospects to improve the current sensitivity of this search strategy. In particular, we indicate the optimal features of potential observational targets, which might allow us to start probing a larger parameter space of ALP DM in the near future. Still, our current analysis presents the strongest constraints that have been obtained so far for the range 10 −23 eV m a 10 −19 eV of ultralight ALP masses.
Our work is organized as follows: In sec. 2, we introduce the concept of ALP, its relevance as a cold dark matter candidate, and the birefringence effect. In the last part of the same section, we discuss the impact of the stochastic nature of the ALP oscillations. In sec. 3 we describe the datasets while in sec. 4 we explain the techniques-based on periodograms-to search for ALP effects on the polarization of radio EM waves emitted in pulsars. In the final part of sec. 4, the individual and combined searches using PPTA and QUIJOTE are presented. Finally, in sec. 5 we summarize some prospects of our studies.
2 ALPs as Ultra Light Dark Matter and the birefringence of light 2.1 ALPs as cold dark matter As mentioned above, one of the most appealing features of ALP dark matter is that it can easily be produced at the right abundance with the misalignment mechanism. The equation of motion of the ALP field a in an isotropic and homogeneous expanding universe reads where m a is the ALP mass, H = (ρ tot /3M 2 Pl ) 1/2 is the Hubble parameter with ρ tot being the total energy density in the Universe. At early times, the ALP field is frozen. Then, when the Hubble parameter becomes comparable to m a , it starts to oscillate and soon after the ALP energy density redshifts as cold matter. At the onset of oscillations, which is characterized by the temperature T osc , the ALP number density is given by where f a is the axion decay constant and θ 0 the initial misalignment angle, defined by the initial value of the axion field a 0 = f a θ 0 . In the conventional cosmological scenario, the ALP starts to oscillate during radiation domination, and the evolution below T osc is assumed to be adiabatic. In this case one can write the axion density parameter as [11,[40][41][42]] The scaling above is easy to understand. A larger decay constant means a larger axion energy density at the start of the oscillations. Furthermore, the larger the mass, the earlier the axion starts redshifting as matter (i.e., slower than radiation), allowing for ρ a /ρ tot to grow over a longer period.

ALP-induced birefringence
The presence of an ALP relic can lead to important phenomenological consequences. For example, its presence will modify the propagation of EM waves. In particular, the ALP medium would rotate the polarization plane of linearly polarized EM waves, such as those emitted by pulsars [30][31][32]43]. The origin of this phenomenon can be understood starting from the Lagrangian for an ALP field a, including its most relevant interactions to light. The latter has the following form where F µν is the EM stress tensor andF µν = ε µνρσ F ρσ /2 is its dual. As follows from dimensional analysis, the coupling constant g aγ has units of inverse of mass. The field equation for the ALP field reads + m 2 a a + g aγ 4 aF µνF µν = 0 , (2.5) and goes together with the modified Maxwell's equations coupled to the ALP field The changes produced by the g aγ a term can be considered adiabatic at scales of the order of the EM wavelength, since we will work only in the limit ∂ t a ω and |∇a| k, ω being the frequency of the radio wave and k its momentum (wave-vector). Then, as far as the propagation of the EM wave is concerned, ∇a and ∂ t a are approximately constant. Terms with first spatial and temporal derivatives squared scale with |∇a|/k and ∂ t a/ω, respectively, so they are also small and neglected. This approach coincides with geometrical optics, where the EM wavelength does not change appreciably over the characteristic length scale of spatial variations of a field. In this context, the EM field admits plane wave solutions with angular frequencies for left and right circularly polarized waves, to lowest order in g aγ [30,31]. Hence, the ALP field acts as a birefringent medium due to the parity-violating nature of the light-axion interaction: a linearly polarized EM wave changes its polarization plane due to the phase shift induced between left-and right-handed circular polarization components in Eq. (2.8). When integrated along the trajectory of propagation of the EM wave, the previous effect yields a modification of the polarization angle by, Two important conclusions follow from this equation. First, this birefringent effect is independent of the frequency of the EM wave. Secondly, it is "conservative" in the sense that the net rotation of the polarization angle depends only on the difference of the values of the ALP field at the source and observer (labeled by times t s and t o respectively) and not on the values of the ALP field crossed by the EM wave in its trajectory 2 . We show in Fig. 1 a graphical representation of this phenomenon. We now describe the changes of the polarization angle of light caused by the difference in the ALP field in the Solar system as compared to the local environment of pulsars as a function of time [32,33,43]. In both cases, the dark matter distribution can be considered as a collection of free waves [11,46,47], where a( k) incorporates the virialized properties of the distribution (in particular its coldness) and includes random phases for each wavenumber. From the dispersion relation derived from Figure 1: Time-dependent birefringence induced by the axion DM field. Light is emitted by a given source (point 1) with a given polarization and received (point 2) with a different polarization ∆φ(t). The effect is "conservative" in the sense that it depends exclusively on the difference between the values of the axion field at points 1 and 2. The time dependence is given by the coherent variation of the axion field, with frequency ν c . (Artist: Ève Barlier.) (2.5) at first order in the background fields, and from the fact that the virialized distribution is non-relativistic (cold), one can write this field configuration at any time t and position x as, a( x, t) = a 0 ( x) cos (m a t + δ( x)) , (2.11) where a 0 ( x) and δ( x) are the amplitude and phase of the field, respectively, at the location x. This description is valid for times smaller than the coherence time where σ ∼ 10 −3 is the dispersion velocity of dark matter in our Galaxy [48] and for distances smaller than the coherence length around the position x. Both properties hinge on the virialized axion field (2.10) containing modes with different nonrelativistic velocities and random phases, which interfere and decohere the net field. To summarize, within the previous limits, the distribution behaves as a coherent field oscillating at the nominal Compton frequency (2.14) Finally, one can relate the axionic field to the local dark matter energy density as derived from the stress energy tensor of an oscillating scalar field [33], where the average corresponds to time averages at cosmological time scales, which are much longer than τ c in (2.12) for the masses of interest.
In the case of interest in our work, the relevant time scale is set by the duration of the observations, which is τ c . In this situation, the connection between the axionic field amplitude and ρ DM in (2.15) acquires a stochastic nature, with the amplitude of the axion field a 0,i following a Rayleigh distribution in a given coherence domain (labeled by i) [46,49]. For the sake of clarity, we write a 0,i = √ 2ρ i m −1 a α i , where we factored out the dependence on the local density of dark matter ρ i and the ALP mass, and α i is a random variable following the probability distribution function (PDF) [46] 3 With the previous ingredients, we can finally study the effect of the axion dark-matter field on the polarization of the light received from a given pulsar as a function of time, taking into account the stochastic nature of the field. From Eq. (2.9), one finds where t is the local time measured since the beginning of the observations (which defines t = 0), T is the time light takes to travel from the pulsar to Earth, and δ i , ρ i and α i are the phases, local dark-matter density and stochastic amplitude variable, respectively, at a given source (i = s) or observer (i = o). The variables δ i are random, with a flat PDF in [0, 2π]. The expression (2.17) can be compactly reorganized as with ∆ = m a T + δ o − δ s . The expression of the phase ϕ a in Eq. (2.18) is not relevant, since it can eventually be considered as a random phase with a flat distribution in [0, 2π] and uncorrelated with φ a 4 . The argument ∆ in Eq. (2.19) can also be treated as a random phase even if the pulsar and Earth belong to the same coherence domain (i.e. when δ s = δ o ). In this case, the phase only depends on m a T , which should also be treated as a random variable since the uncertainty in the measurement of the distance typically spans several cycles of the ALP field.
In summary, one can indirectly search for the interaction of light with the dark-matter waves of the ALP field by looking for periodic variations in ∆φ(t), Eq. (2.18). There are two important caveats to this strategy. First, periodic variations in polarization measurements may be caused by more prosaic mechanisms, such as seasonal effects at Earth (e.g., Faraday rotation in the ionosphere [37,51]), astrophysical phenomena or even instrumental systematic effects. Some of these contributions are, in principle, well known and can be corrected for, like the ionospheric corrections or some of instrumental origin. Others are potentially unknown and would need to be modelled and subtracted before using ∆φ(t) effectively to search for DM. Secondly, the analysis depends on unknown environmental variables such as ∆ or the stochastic amplitudes α i . This may cause a loss of sensitivity to the ALP birefringence effect in the unlucky situation of near-zero field amplitudes or cancellations in Eq. (2.19) during our measurements.
Both problems can be circumvented by performing a global analysis of several sources measured, ideally, at different frequencies and with independent experiments. Using several sources allows one to discard astrophysical periodic effects, as the ALP-induced periodic changes in the polarization of the EM wave should have the same period for all sources, while those of astrophysical origin are expected to have source-specific periods and patterns. On the other hand, different frequencies and experiments allow us to eliminate frequency-dependent changes (e.g. Faraday rotation [51]) and instrumental effects, respectively. Finally, a global analysis enables us to tackle the lack of knowledge about the environmental variables by effectively sampling them statistically from their PDF, such as Eq. (2.16) for the stochastic amplitudes.
Before moving to the bulk of the analysis, it is important to stress that the stochastic behavior of the axion field is a crucial ingredient to obtain realistic bounds, often overlooked in previous analysis. In the following, to further stress this point, we refer to results both assuming a fixed α i = 1 ("deterministic" scenario, applicable only for T obs τ c , where the experiment explores the whole distribution of the field values) or by sampling it from the Rayleigh distribution in Eq. (2.16) ("stochastic" scenario, correct for T obs τ c ).

Birrefringence searches and UV models of ALP DM
As discussed in Sec. 2.1 and shown in Eq. (2.3), the cosmological abundance of ALP DM is set by its mass and decay constant. The latter is, in turn, related to the axion-photon coupling via an equation of the type where α em is the EM fine structure constant and κ is a model-dependent coefficient. Therefore, once f a is adjusted for a given mass m a to provide for the relic DM abundance, one can interpret our search for ALP DM in terms of the sensitivity to different models. In typical QCD-axion models κ receives two contributions: an infrared (low-energy) contribution of O(1) coming from mixing between axion and QCD mesons, and an ultraviolet (high-energy) contribution which usually depends on the electromagnetic and QCD anomaly coefficients. The latter is usually also taken to be of O(1); however, it can be, in principle, much larger than the infrared contribution, leading to an enhanced axion-photon coupling. This occurs in several scenarios, including the photophilic axion of Ref. [52], where the axion has an exponentially large coupling to photons thanks to the "clockwork" mechanism, the generalized Kim-Nilles-Peloso (KNP) alignment scenarios discussed in Ref. [53], or through a kinetic mixing of abelian gauge fields, as considered in Ref. [54], where the axion-photon coupling inherits the (potentially large) dark photon gauge coupling. Therefore, searches of axion DM probing their coupling with photons g aγ can always be reinterpreted as measurements or upper limits of the κ coefficient (see [42] for a related discussion).

Data
As discussed in Sec. 2.2, in order to perform a robust search of ultra-light axion DM through its birefringence effects, it is crucial to combine many different sources and, if possible, different experiments observing at different wavelengths. In this work we use an ensemble of observations of the polarization angle of 21 Galactic sources measured by the Parkes Pulsar Time Array (PPTA) at about 1.4 GHz [37] and by the QUIJOTE experiment between 10 and 20 GHz with its MFI instrument.
In Table 1 we show the list of the pulsars used in this work and the parameters relevant to our analysis. In particular, knowing the location of the pulsars with respect to Earth and within the Milky Way allows us to estimate the local dark matter density at each source and also determine whether the distance between the source and the observation point is smaller than the coherence length l c . In Table 1 we show the distance of all the pulsars used in our analysis to Earth [55][56][57][58][59] and to the Galactic center. We also show the corresponding ρ DM at the pulsar estimated using a Navarro-Frenk-White (NFW) profile [60] with the parameters extracted from Ref. [61]. Note that the pulsar distances to Earth, which are crucial to locate them in the Galaxy, are determined by different methods and in some cases have large uncertainties. Nonetheless, even in the most extreme cases this is translated in an uncertainty never larger than 5% in ρ DM by taking the 1σ band in the corresponding distance. Given that the birefringence effect is ∝ ρ 1/2 DM , Eq. (2.17), this uncertainty leads to changes of ∼ 2 − 3% that we neglect in our analysis. We omit for the same reason assigning errors to the values of the distances quoted in Table 1.
Our goal is to combine long-time measurements of the polarization angle of these pulsars to search for a global modulation ∆φ(t) produced by the axion field oscillations, scanning its frequency ν c in a certain range that we will determine from the structure of the observation windows below. In the next two subsections we describe the PPTA and QUIJOTE measurements on these sources.

PPTA data
The Parkes Pulsar Timing Array (PPTA) started in 2005 and uses the Parkes 64 m antenna (Parkes, Australia) to perform continuous observations of a sample of 20 ms pulsars in three radio frequency bands centred at 0.7 GHz, 1.4 GHz and 3.1 GHz, with intervals from two to three weeks, with the primary goal to detect the background of gravitational waves in the pulsar timing band [62]. In this work we use PPTA data released in Ref. [37] that include polarization measurements of 20 pulsars at 1.4 GHz. These data span a total observation time of ∼ 4.5 years with a precision that can be as high as O(0.1 • ) for some pulsars. Further technical details about these observations can be found in [37]. At 1.4 GHz changes in the Faraday rotation arising in the Earth's ionosphere lead to variations of the polarisation angle in scales of one day that are way larger than the data statistical error bar. This effect must be modeled and subtracted from the data [37]. The PPTA datasets implement two models to treat ionospheric corrections: the 2007 International Reference Ionosphere (IRI) model-GETRM−IONO [63] and the FARROT (Faraday Rotation) software. The latter was developed by the Dominion Radio Astrophysical Observatory (DRAO), Penticton, Canada. FARROT uses three Chapman layers and a semi-empirical model for variations in these based on the 10.7-cm solar radio flux measured by DRAO. Instead, IRI uses the model from the International Union of Radio Science (URSI) [64]. FARROT and IRI models incorporate the International Geomagnetic Reference Field2 to take into account the effects from Earth's magnetic field [65]. These models allow one to study the changes in the rotation measure and in the polarization angle of our signal due to ionospheric variations. In this work we analyze the PPTA data corrected by the IRI model because it gives better results than FARROT [37] (see discussion below).

QUIJOTE data
PPTA data are supplemented by observations of the Tau A supernova remnant (Crab Nebula) performed with QUIJOTE [38,39] (hereafter Crab-QUIJOTE), an experiment dedicated to observe the cosmic microwave background polarization between 10 and 40 GHz from the Teide Observatory (2400 m a.s.l., Tenerife, Spain). In particular we use data at 11, 13, 17 and 19 GHz taken with the QUIJOTE Multi-Frequency Instrument (MFI) [66], which was operative between 2012 and 2018. Owing to it being the brightest polarised compact source on the sky in this frequency range, Tau A was the main reference calibrator of the MFI. As such, this source is observed on a regular (virtually daily) basis, using observations in raster-scan mode, lasting ∼ 25 min to produce maps of size ∼ 8 • × 8 • around the source to reach deep sensitivities. Tau A was also observed during the QUIJOTE-MFI wide-survey of the full northern sky [67], using the so-called "nominal mode". In this case the telescope is continuously rotated in azimuth at a constant elevation, resulting in a map of a large area of the Northern sky after a full day. Due to its declination (22.0 • ) being very close to Tenerife's latitude (28.3 • ) Tau A was detected in all these observations.
In Figure 2 we show two examples of QUIJOTE-MFI 11 GHz maps of Stokes Q and U of Tau A, using individual observations in raster and nominal mode. Due to differences in the two scanning strategies, the former leads to a higher integration time per unit area, which results in an improved map sensitivity that is apparent by eye. We gathered a total of 1363 and 685 Tau A observations between November 2012 and December 2018, respectively in these two modes, for a total of 2048 Tau A maps of Stokes Q and U. In each individual map we apply an aperture photometry technique, consisting in integrating the Q and U flux densities enclosed by the solid circle in Figure 2 and subtracting a median background level calculated in the region inside the two concentric dotted circles. The polarisation direction is then obtained as γ = 0.5 arctan(−U/Q) 5 . Analogous maps are produced for the other MFI frequencies, so that for each observation we have four estimates of the Tau A polarisation angle at 11, 13, 17 and 19 GHz. MFI has four horns, each of which observe a pair of frequency bands. In this analysis we have used the 11 and 13 GHz frequency bands of horn 3, and the 17 and 19 GHz frequency bands of horns 2 and 4. The 11 GHz provides the best sensitivity, while, due to the combined effect of the decay of the Tau A flux density with frequency (spectral index α ∼ −0.3) and to the proximity to the atmospheric 22 GHz water vapour line, the signal-to-noise gradually worsens at higher frequencies. We compute for each observation a weighted average of the four polarisation angle estimates, justified by the fact that the ALP birefringence effect is frequency independent (see Sec. 2.2). We reach sensitivities in the determination of the polarisation angle of O(1 − 2 • ) for each observation in raster mode, and a bit worse in the case of the nominal mode.
The QUIJOTE-MFI calibration on Tau A consists in setting the global polarization direction, by combining many individual observations. The instrument polarization direction is stable during large periods of time (∼ 1.5 years), and therefore variations of the observed polarization angle within these periods could only be due to real on-sky variations. Possible variations in time scales larger than these ∼ 1.5-year long periods were checked through comparison with the polarization direction measured by WMAP at 22.8 GHz and Planck at 28.4 GHz on the diffuse Galactic emission, and were found to be smaller than ∼ 0.5 • . Further details on this study will be provided in different upcoming papers ( [67] and [68]).
The QUIJOTE-MFI instrument underwent different technical modifications during its 6-year lifetime, leading to differences in the effective instrument sensitivity as a function of time. In this study we decided to remove data taken before April 2014, which was less precise due to the lack of a reference calibration diode, and we end up with a sample of 1099 QUIJOTE-MFI estimates of Tau A polarisation direction, spanning a period of 4.5 years.

Generalized Lomb-Scargle periodograms
The best approach for our endeavor and, in general, to search for periodic signals in unevenly distributed time series is the Lomb-Scargle (LS) periodogram [69,70]. This is a standard tool in astronomy based on a combination of techniques of signal analysis: Fourier methods, least squares, phase-folding, and Bayesian approaches [71]. The central element of the method, the LS periodogram, which we label P LS (ν) in the following, can be regarded as the analog of the power spectrum in continuous Fourier analysis (see Ref. [71] for subtleties of this interpretation). It can also be understood as the result of a least-square minimization process, which can be extended to include the uncertainties of the data. This leads to the generalized LS periodogram [71,72], which is the version that we use in our work. In Appendix A we  Table 1: List of pulsars used in our work including parameters relevant for our analysis: Central value of measured distances to Earth, D Earth , distance to the Galactic center, D GC , and local dark matter density, ρ DM , calculated with a NFW profile with parameters extracted from [61] and observation time (time difference between last and first observation), OT, for each pulsar used in the analysis. See the main text for a discussion related to the impact in our analysis of the uncertainties in the distance measurements.
provide a brief description of the generalized LS periodogram, the standard procedure we follow to calculate it and its relation to a statistical minimization procedure. The significance of any peak appearing in P LS (ν) can be assessed by a False Alarm Probability (FAP). The FAP quantifies the probability that a data set with no periodic signal leads to such peak by coincidental alignment of random fluctuations [71]. The most robust method to compute the FAP is the bootstrap approach, which resamples the time series by keeping the temporal coordinates but with observations randomly distributed using Monte-Carlo (MC) simulations [73]. The maximum peaks obtained in these MC simulations form a distribution, which can be used to directly compute the FAP.
To apply the LS technique and extract reliable bounds, it is important to consider the frequency limits and the grid spacing of our data [71]. For a pulsar observation spanning a total time T obs , a signal with frequency 1/T obs will complete exactly one oscillation cycle, providing a suitable minimum frequency ν min for the LS periodogram 6 . For the maximum frequency ν max we choose the inverse of the folding time (i.e., resolution time for observations) that we identify with the pseudo-Nyquist frequency (see [71] for more details). Taking as a reference the observations of J0437-4715 by PPTA [37], ν max 20 days −1 is determined by the folding time of ∼ 1 hour and ν min 6 × 10 −4 days −1 by T obs ∼ 4.4 yrs. As a result, the range of ALP masses that can be probed spans about four orders of magnitude, 10 −23 m a 10 −19 eV. Once the frequency interval is identified, one must choose a frequency resolution. The grid must have enough resolution to avoid the possibility of losing peak power. However, a large number of points may make the search for peaks-and hence the production of bootstrapsprohibitive in computational cost. To ensure a frequency grid of good resolution for each peak, we over-sample by a factor n 0 = 5 per peak, and use a grid of spacing ν min /n 0 . With this choice, the total number of required periodogram evaluations is N = n 0 T obs ν max [73,74].

Searching for peaks in the power spectrum
After adjusting the range and density of the grids, we compute the generalized LS periodograms corresponding to the polarization time series of the 20 PPTA pulsars (implementing the IRI ionospheric corrections) and of Crab-QUIJOTE. In each case we also compute the FAPs with false-probability rates of 32%, 5% and 1%, using 1000 MC bootstraps of the data. For the PPTA pulsars we observe that the highest values of P LS (ν), and those of the FAPs, are smaller when the data set is more precise. In general, we do not find peaks with FAPs smaller than 1% (or even smaller than 5% in most cases). There are two exceptions: J1730-2304 and Crab-QUIJOTE, for which we find peaks with false-rate probability similar or smaller than 1%. This is illustrated in Fig. 3, where we show the periodograms for two of the statiscally most powerful 7 PPTA pulsars, J0437-4715 and J1643-1224, and for the two pulsars where prominent peaks appear. It is important to note that the periodograms of the time series without ionospheric corrections show prominent peaks at ν ∼ 1 day −1 and ν ∼ 1 yr −1 for almost all the PPTA pulsars. This is precisely the type of seasonal pattern that is expected to be produced by ionospheric effects. Therefore, subtracting them is important for our analysis, albeit at the cost of introducing some dependence on the corresponding modelling. In particular, the peak found in J1730-2304 seems to be spuriously introduced by these models as it does not appear in the uncorrected data (see Appendix B for details about a LS analysis of these ionospheric corrections). In the case of Crab-QUIJOTE, we observe two peaks surpassing FAPs significance of 1%. This indicates the presence of periodic signals, the most significant of which corresponds to ∼ 5 day −1 .
The absence of similar peaks in the power-spectrum of all the other pulsars, especially in the statistically most powerful ones, such as J0437-4715 and J1643-1224, suggests that the peaks observed with low FAP are not produced by the oscillations of the ALP field. In order to address this question precisely we first need to translate the absence of peaks for a given source into an upper limit on the value of the amplitude of the oscillation. Finally, comparing and combining different pulsars requires expressing this bound as an upper limit on the axion-photon coupling g aγ .

Constraints from single sources
Next, we perform a new set of MC simulations generating 3000 pseudo-experiments, as done in the bootstrap method for the calculation of the FAPs, but injecting a harmonic oscillation in the time series ∆φ sim = φ cos (2πν t + ϕ) , (4.1) with a given frequency ν and amplitude φ, and a random phase ϕ sampled from [0, 2π] with a flat PDF. With this population of pseudo-experiments, we compute the PDF of the corresponding LS periodograms P LS (ν, φ) and find the value corresponding to the lower 5% tail of the distribution (i.e., such that the probability of obtaining a P LS (ν, φ) larger than that value is 95%) that we callP LS (ν, φ). We then find numerically the value φ 95 for which P LS (ν, φ 95 ) is equal to the experimental P LS (ν). Hence, φ 95 can be interpreted as an upper bound on φ at 95% confidence level for the given frequency ν 8 . To smear out the fluctuations related to the MC realizations over the LS periodograms, we smooth out the limits considering segmentation of the bound using a rolling mean method 9 : at a given frequency ν j , the new output for φ 95 (ν j ) is the average of the φ 95 within a window centered at ν j , where we take as 7 In the sense of smallest nominal errors and dispersion of the data around the mean. 8 Important aspects of this method are similar to the techniques used to detect tertiary companionsexoplanets in binary systems, see Ref. [75]. 9 In the context of nonlinear signal processing, Ref. [76] gives a general discussion of the rolling mean method (RMM). On the other hand, it is also possible to relate the fluctuations to statistical variances of the periodograms using Barttlet's method [77], which gives a smoothing of the LS periodograms through averaging with some filters over the neighbouring power spectrum frequencies ending in a method similar to RMM. In MC-LS periodograms lower points arise in different frequencies (windows) yielding stronger bounds of φ95, hence the most accurate method to treat the under-fluctuations is one with dynamical windows as given by the RMM. the reference the φ 95 at the frequencies ν j−1 and ν j+1 . Now the output for φ 95 (ν j+1 ) is the average of the φ 95 (ν j ) and φ 95 (ν j+2 ) within the window centered in ν j+1 . This procedure is repeated for all j with the same window size for all grid frequencies and all pulsars. Thus, at each frequency ν j , the rolling mean computes a new estimate for φ 95 . Figure 4 shows individual φ 95 bounds for seven pulsars from Table 1, corresponding to the four PPTA pulsars with the most competitive bounds (J0437-4715, J1600-3053, J1643-1224, and J1824-2452), two PPTA pulsars with the weakest limits (J1024-0719, J1730-2304) and Crab. At low frequencies, the pulsar J1600-3053 gives the strongest bound in the range of 6× 10 −4 day −1 ν 10 −3 day −1 . Instead, for frequencies greater than 10 −3 days −1 , the pulsar J0437-4715 bounds dominate over those of the other pulsars, giving the strongest constraint on φ 95 . The rest of pulsars from PPTA have constraints among these limits. Finally, the analysis of Crab-QUIJOTE leads to bounds in the ballpark of those from the weakest PPTA pulsars. Nonetheless, the difference of the bounds obtained from the statistically stronger and weaker pulsars is not dramatic and is generally not larger than a factor ∼ 5.
We emphasize at this point that, for each pulsar, φ 95 is just an experimental limit on the amplitude of a harmonic signal that is assumed to appear overlaid on a random resampling of the data. An upper limit on the axion-photon coupling g aγ can be obtained by assuming that this harmonic signal is entirely produced by the ALP DM field oscillations; namely, φ a = φ in Eq. (2.18). This is a conservative assumption at the level of setting constraints, as we are allowing the ALP to saturate any possible periodic signal in the data. We can then use Eq. (2.19) to extract a bound on g aγ from the one on φ a . The amplitude φ a is a random variable whose PDF can be extracted from the MC studies and the data as described above. For simplicity, we model this PDF for each source and frequency as a Rayleigh distribution with a scale parameter σ(ν) adjusted to match optimally the PDF obtained directly from MCs 10 . The PDF of g aγ then follows from the one of φ a and from those of the environmental variables: ∆ in the deterministic case or ∆, α o and α s in the stochastic case. Surprisingly, we find that in our setup the sensitivity to g aγ improves in the stochastic analysis compared to the deterministic one. We refer the reader to Appendix C, where we perform a toy MC study of this comparison. In the following, we derive the bounds on g aγ using the stochastic analysis.
For each PPTA pulsar and Crab-QUIJOTE, we generate 10 5 MC pseudo-experiments. In each one we draw a measurement of φ a , with the corresponding distributions adjusted above, and a value of the environmental variables ∆, α o and α s . Here we assume that the source and Earth are in different spatial coherence patches (D Earth l c ). By comparing the data in Tab. 1 and Eq. (2.13), one can see that this assumption is valid for all pulsars and the whole mass range except for J0437-4715 and J0711-6830 in the mass region close to the lower end m a 10 −23 eV. Given that these pulsars do not lead to the strongest constraints in that region, see Fig. 4, we neglect the effect of spatial Earth-source coherence in our analysis.
Using Eq. (2.19) and the values of ρ DM displayed in Tab. 1, we obtain the individual PDFs of g aγ as a function of m a . The strongest bounds are derived from the four statistically most powerful PPTA pulsars, J0437-4715, J1600-3053, J1643-1224, and J1824-2452 shown in Fig. 4. The relatively weaker bound derived on φ a from J1824-2452 is compensated by the higher ρ DM density estimated at source. In the end, these four pulsars lead to relatively similar bounds in all the considered mass range. For instance, at m a = 10 −22 eV one obtains g aγ 10 −12 GeV −1 at 95% C.L., which is in the same ballpark as the limits derived from other astrophysical probes (see below).
Finally, note that the discussion here is purely statistical and assumes not only that we can repeat the experiments for a given pulsar (φ a ), but also that we can draw different values for the environmental variables ∆ and the α i . Since these are unknown but fixed parameters, the only way to sample their distributions in our given Universe is by performing a global analysis including many sources, as done below.

Constraint from a combined analysis
We are now ready to derive the combined bounds from the PPTA and Crab-QUIJOTE pulsars. In order to do this we carry out a MC analysis, as in the previous section. In each pseudoexperiment (or MC run) we obtain a value of g aγ for each PPTA pulsar and Crab-QUIJOTE. The global value g aγ is obtained as a weighted mean. By running several pseudo-experiments, we are able to construct its PDF.
To be more precise: Under the Rayleigh model for the distributions of φ a obtained in the previous subsection, we compute the scale σ i (ν j ) for each pulsar (i = 1, ..., N pulsars ) and each frequency we are probing. We then generate 10 5 MC pseudo-experiments drawing a value of φ a,i for each pulsar. We also draw a value of the environmental variables α o (for the ALP amplitude at Earth) and the α s,i and ∆ i for the pulsars. Using Eq. (2.19) we obtain the values of g aγ,i for each pulsar. The weights for the weighted mean of the g aγ (or global average) are found by propagating the root squares of the variances of φ a,i , which are directly related to the scale factor in the Rayleigh distribution. We report our upper limit on the axion-photon coupling as the 95% C.L. value obtained from the resulting distribution of the global g aγ . As an illustration of the impact of combining several pulsars, in Fig. 5 we show the combination of the 4 pulsars with the strongest individual φ a bounds (see Fig. 4). We observe that as more pulsars of similar sensitivity are added into the combined analysis, the reach of the bound increases. The pulsar J1824-2452 contributes in a significant way to the combination due to its higher ρ DM , despite not being as statistically powerful as the others.
We are also now in the position to analyse the peaks in the periodogram discussed in Sec. 4.2, appearing in J1730-2304 and Crab-QUIJOTE with FAPs overpassing (or close to) 1% significance. In particular, for a given peak we compare the PDF of g aγ at the given frequency and pulsar with the global PDF of g aγ at that same frequency and obtained excluding that same pulsar. In so doing, we find that the hypothesis that the peak observed by J1730-2304 at ν 3 days −1 or by Crab-QUIJOTE at ν 5 days −1 is produced by the interactions of the photons with the ALP field is excluded by the global analysis in both cases with a significance larger than 99.999% 11 .
In Fig. 6 we show the final bound including all 20 pulsars from PPTA and Crab-QUIJOTE. The 95% C.L. upper combined limit for the ALP coupling is comparable to other bounds obtained by searching for the same effect either in parsec-scale jets from active galaxy nuclei by MOJAVE VLBA [33] or from CMB by BICEP-Keck [79] (for the origin of these acronyms see the corresponding references). The bound on g aγ derived in this paper is the strongest for all the masses in the range 3 × 10 −23 eV ≤ m a ≤ 2 × 10 −21 eV. In particular, in the mass regions around m a 3×10 −23 eV and m a 10 −22 eV, our limit is stronger than the other ones by a factor ∼ 5. For masses approximately above 2 × 10 −21 eV the strongest limit is set by supernova SN 1987A cooling via Primakoff effect, via a null search of γ rays [26]. The 11 More precisely, the value that we are excluding is the median of the PDF.  Figure 6: Bounds on the axion-photon coupling as a function of the ALP mass. The blue area indicates the excluded region at 95% C.L. obtained in this work from the combination of PPTA and Crab-QUIJOTE pulsars including stochastic effects. The darker gray dot-dashed line indicates the region excluded by CAST experiment [25]. The gray solid line is the limit from supernova SN 1987A [26]. The orange area indicates the region excluded by MOJAVE VLBA. The green region is the excluded zone by BICEP-Keck (smoothed bound) [79].
strongest limit by a direct experimental search is set by the CAST experiment, by searching for axions emitted by the Sun [25].
It is important to stress that neither the analysis of VLBA-MOJAVE nor BICEP-Keck incorporate the stochastic nature of the ULDM field amplitude. As we have discussed in our analysis, this might have a significant impact in the quoted upper limits. Finally, we also note that our bounds are considerably weaker (and extend over a shorter range of masses) as compared to those derived from J0437-4715 in [32] by some of the authors of this paper, using the same dataset. The main reasons for this difference are: (i) the robust statistical treatment provided by the FAPs in a periodogram analysis of time series with large dispersion of central values and small nominal errors, in contrast to a standard frequentist analysis; and (ii) a proper treatment of the environmental variables in a global analysis of data from different pulsars. In Appendix D we present a more detailed discussion of this comparison. In summary, the limits from the PPTA and Crab-QUIJOTE pulsars on the axion-photon couplings presented in this work update and improve (in terms of quality and robustness of analysis) those presented in [32].

Future prospects
Let us finally comment on future directions and possible ways in which our analysis can be further improved. One simple extension would be to consider data from pulsars closer to the for 20 Galactic pulsars, as presented in Fig. 6 (blue), and for the future prospect of observing 10 pulsars with similar precision to J0437-4715 during a span of 5 years located at the same distance to the Galactic center as J1745-2900 (green). The black dashed line is the limit from supernova SN 1987A [26]. In orange the values of g aγ for an axion that makes up 100% of the DM density following Eqs. Galactic center region such as J1745-2900 [80]. This pulsar would be immersed in an ALP DM medium with density ρ DM of 46.75 GeV/cm 3 , assuming a distance to Earth of 8.3 kpc as in [81] and using the same DM profile as above. One of the problems with this strategy is that the presence of hot, ionized gas in the central part of our Galaxy may lead to a decrease in the observed flux density from neutron stars residing in this region. Nevertheless, future surveys such as the Square Kilometre Array (SKA) [82] and Next Generation Very Large Array (ngVLA) [83] are expected to probe a sizable population of pulsars in the Galactic center. This should provide better constraints because of the larger DM density in those regions (the limit on g aγ scales linearly with the square root of the density), and also robust ones because of the large number of expected detections. In Fig. 7 we show a forecast of the potential reach of an analysis of 10 pulsars, each including 500 observations spanning 5 years with a dispersion and precision similar to the J0437-4715 data, and with ρ DM = 46.75 GeV/cm 3 . Another intriguing possibility is to consider polarization data from pulsars hosted by globular clusters. At the moment there is no consensus on the DM profiles of globular clusters, see e.g. [84][85][86]. Nevertheless, it is plausible that order 1% of their content is still made of DM today [87], after being drastically reduced because of tidal stripping and dark matter thermalization [88] 12 . Even if only O(1)% of the total amount of matter in globular clusters is dark matter, its density in the inner cores may be O(10 2 − 10 3 ) larger than the local one. This would imply again a sensitivity to an ALP-photon coupling 10 − 30 times smaller than the one probed in the present work, assuming the same precision on the measurements. Furthermore, globular clusters provide a natural environment for pulsars to form, see [91]. In both cases, it is very likely that future bounds will reduce the allowed parameter space by more than an order of magnitude in the coupling g aγ , an achievement hard to foreseen with any other probes for this range of masses.

Conclusions
We have studied the effects of ULDM on the polarized radio waves that we detect from pulsars. We focused on the coupling g aγ in (2.4), which makes the ULDM behave as a chiral medium harmonically affecting the polarization angle of photons while they are in transit from the astronomical source to Earth. Our work improves on previous analyses in several ways. First, we use techniques based on periodograms to analyze the time series of the polarization measurements. The calculation of the false-alarm probability rates, using bootstrap (Montecarlo sampling) techniques, provide us with robust upper limits on the strength of harmonic signals in the data. Second, we include systematically the environmental variables in the analysis, see Eq. (2.19). In particular, we carefully considered the intrinsic degree of stochasticity of this signal, arising from the virialized state of ULDM in the Milky Way's dark matter halo. At the practical level, this transforms the local energy density ρ DM into a random variable Rayleigh-distributed with scale parameter given by the standard value in WIMP DM models, with a different value in each coherent patch of size l c in Eq. (2.13). And, finally, we perform a global analysis of 20 pulsars measured by the Parkes PTA and of the Crab pulsar measured by the experiment QUIJOTE for ∼4.5 years. This allows us to circumvent statistically our lack of knowledge of the exact values of the environmental variables by sampling their known probability distributions. In addition, it allows us to exclude possible harmonic signals of a more prosaic origin, such as source-specific astrophysical phenomena. Including two different experiments in the analysis is also important, to search for effects appearing at different frequencies and for excluding possible signals of instrumental origin.
All in all, our analysis leads to the strongest combined constraints on the ALP ULDM scenario for the range of masses 10 −23 m a 2 × 10 −21 eV, as shown in Fig. 6. In sec. 2.3 we discussed the role of these bounds for models of ALP DM, where the relation between the coupling g aγ and the mass m a (connected to the ALP decay constant f a ), is encapsulated by κ, considered a free parameter.
Finally, by using an idealized model of a pulsar near the Galactic center in a region with a high energy density of DM (and also pulsars in globular clusters), we have highlighted some prospects for probing further the ALP parameter space.
The code to reproduce these analyses can be found at https://github.com/jorgetc16/ ALP_meter.

A Generalized Lomb Scargle Periodogram
Here we define the Lomb-Scargle periodogram P LS (ν). Consider a time series with N observations, where y = {y i }, σ = {σ i } and t = {t i } are the data, errors, and time "vectors" respectively. The normalized errors weights are where W = 1/σ 2 i , and with w i = 1. We denote for each frequency C = {c i } = {cos (2πνt i )} and S = {s i } = {sin (2πνt i )}. In addition, we define the products [33,72] The generalized LS periodogram with errors weights is then determined as where 0 ≤ P LS (ν) ≤ 1. Note that the dot product between two "vectors" (in bold face) expresses their scalar product and between two scalars it expresses their multiplication. An alternative way to understand the LS periodogram is as a statistical minimization analysis given by fitting an harmonic model [71,72] y = a ν cos(2πνt), where the amplitudes a ν are determined by a χ 2 minimization process for the different frequencies ν considered (one can also add an offset c ν to the fit [71]). With this procedure, one can recover the generalized LS periodogram, P LS (ν) as whereχ ν is the minimum value for that frequency ν andχ 0 is the minimum value for a non-relevant reference model.

B Effects of the ionospheric corrections
As an example of the impact of FARROT and IRI models of the ionospheric corrections, in Figure 10  We also show FAPs at 32% (red solid line) and 5% (red dashed line) false-positive rates estimated using a bootstrap method with a 1000 random resamplings of the data set at the same temporal coordinates.
peaks clearly disappear once the corrections are applied. Instead in the case of J1730-2304, the modelling of ionospheric corrections seems to introduce "spurious" peaks over FAPs of 5%, which are absent in the original data without ionospheric corrections. Although this affects the background model for the pulsar J1730-2304 and its individual bound on φ 95 , the combined bound exclude these "fake" peaks effects in the ALP signals searches. Since the IRI method gives better results in modelling ionospheric corrections (lower seasonal peaks and less spurious ones), we use the data implementing this model in our analyses.

C Comparison Stochastic vs. Deterministic analysis
In general, one might expect a degradation of the bounds on g aγ in the stochastic case compared to the deterministic one. In order to understand this effect quantitatively, we perform toy MC simulations in two limits that are relevant in our analysis: (i) homogeneous, where the DM density at source and observer are the same ρ s = ρ o , which is the configuration we approximately have in our dataset (see Tab. 1); and (ii) overdense, when the DM density at the source is much higher than the one at Earth, ρ s ρ o , which is relevant for our forecasts 13 . The amplitude of the ALP-induced oscillation in Eq. (2.19), in these two cases and implementing the stochasticity of the ALP field, simplifies as In the deterministic case (α i = 1) this simplifies further to We perform several pseudo-experiments where we sample a measurement of φ a , ∆ and of the α i and we derive the limits in the deterministic and stochastic cases to compare them. For this toy MC we assume observations φ i over i = 1, . . . , N s sources all featuring the same precision (that is, the same scale factor for the Rayleigh distribution) and with the same ρ DM .
Surprisingly, we find that the upper limit on g aγ improves in the stochastic analysis by a factor 2.6 in the homogeneous case, while it worsens (as expected) by a factor 0.55 for the overdense one. The improvement in the homogeneous limit seems to be due to a nontrivial interplay between the random ∆ and α i variables. The impact of stochasticity is intuitively more severe for the overdense case, as the sensitivity of a given pulsar to g aγ hinges only on the particular value of α s (compared to the former where both α s and α o are relevant). If we neglect the effect of ∆ by e.g. taking the root mean square value of Eq. (2.19) (as often done in the literature [34,43]) then we would obtain a limit that is worse by a factor 0.92 also in the homogeneous case.
Expanding on the toy example, we repeat the pseudo-experiments for N s sources, obtaining a "global" value of g aγ as the weighted average of those measured at each source. In the stochastic case, g aγ ∝ (α 2 o + α 2 s − 2α o α s cos ∆) −1/2 , and a larger α s or α o will translate into a smaller g aγ . When this happens, the variance of g aγ that determines the weight in the average also shrinks, giving a higher weight to that measurement. Moreover, in the case of ∆ → π, the interference with the cos ∆ term will be constructive and its effect become greater. It is not likely to get a high value of α s and a ∆ close to π for a single pulsar; however, as we add more sources the probability of getting at least one with a high value of α s increases. In the deterministic case, g aγ ∝ (2 − 2 cos ∆) −1/2 is bounded by 2, but as N s increases these larger denominators (and therefore smaller g aγ and variances) will dominate the weighted average, making the resulting PDF thinner and closer to the stochastic one.
In the overdense scenario, for the stochastic case, we have g aγ ∝ α −1 s (the interference term is subleading). In the deterministic case, the environmental variables have no effect and the PDF will be constant no matter how many pulsars one adds. This is why we see a stronger bound in the deterministic case for a small number of N s while, as N s increases, the probability of getting a high value of α s increases, and the stochastic bound will eventually overcome the deterministic one.
All these effects are seen in Table 2, where we show the size of the bounds at 95% C.L. obtained in the deterministic scenario relative to those obtained in the stochastic one as a function of N s in the two cases discussed above. Interestingly, in the homogeneous case the difference becomes small as more sources are added to the analysis. For the overdense case, which is relevant in the forecasts presented in Sec. 5, a boost in the sensitivity could be achieved by studying many sources with similar statistical power.  Table 2: Ratio of the bounds at 95% C.L. obtained in the deterministic scenario over the ones obtained in the stochastic scenario for different number of sources, N s .

D Comparison between periodograms and the frequentist method
The frequentist method used in [32] relies on the comparison between the χ 2 obtained for a fit to a function with some parameters, in this case to an harmonic oscillating function, and the χ 2 of the null hypothesis, i.e., assuming any fluctuation in the data is due to statistical noise. A value of the parameters of the fit will be ruled out at certain C.L. when ∆χ 2 ≡ χ 2 null −χ 2 fit > N , with N depending on the number of parameters of the fit.
As discussed in Appendix A, the value of the periodogram at a given frequency is equivalent to the result of a minimization of a χ 2 with also a harmonic function at that frequency. However, in the periodogram approach the C.L. regions are derived by using the FAP procedure, and not assuming a χ 2 distribution given by the statistical nominal errors of the data. This is certainly not the case as the χ 2 /d. o. f found for both the null hypothesis and the periodic function are extremely large (e.g. χ 2 /d. o. f is 488.6 for the polarization angle of J0437-4715 measurements using IRI corrections), indicating that none of them is a good fit to the data. Namely, for some combination of the parameters, one finds that the χ 2 fit is significantly smaller than χ 2 null , from which one could deduce that the presence of a periodic signal is favoured by the data with a very high frequentist C.L. Still, both fits are disfavoured by the large χ 2 /d. o. f, which should be taken into account in the final bounds. On the other hand, the periodogram approach, with FAPs estimated using bootstrap methods, provides a robust estimate of the C.L. for the presence or absence of harmonic signals in the time series of the polarization measurements.
In Fig. 9 we compare the deterministic limits reported for the J0437-4715 using the frequentist approach in [32] with the ones that we obtain using the periodograms. One can see that the frequentist interpretation of the C.L. can lead to limits an order of magnitude stronger than the more realistic ones given by the FAPs.
In terms of the range of frequencies considered, for J0437-4715 the range of masses studied in [32] goes down to m a = 5 × 10 −24 eV, which is equivalent to a frequency of ν = 10 −4 days −1 or a period of T = 26.6 years. As discussed in Sec. 4.1, periodic harmonic signals with such low frequencies can be confused with a linear signal. Since our method critically relies on the study of periodic signals in the data, it is not possible for us to explore the ultra low mass regime with periodograms, as we can see in Fig. 9.
Finally, in [32], only one of the available pulsars was taken into account. As discussed in the main text, using several sources is important, not only because the signal of ULDM at a given frequency should appear in all sources, but because it diminishes the stochastic effects that affect the constraint. This work Caputo et al. Figure 9: Bounds on the photon-ALP coupling at 95% C.L. derived with the pulsar J0437-4715. The orange area corresponds to the bounds derived with the method in [32], while the blue region is derived in this work (under the deterministic assumption).