# Searches for signatures of ultra-light axion dark matter in polarimetry data of the European Pulsar Timing Array

N.K Porayko, P Usynina, J Terol-Calvo, J Martin Camalich, G.M Shaifullah, A Castillo, D Blas, Lucas Guillemot, M Peel, C Tiburzi, et al.

► To cite this version:

N.K Porayko, P Usynina, J Terol-Calvo, J Martin Camalich, G.M Shaifullah, et al. Searches for signatures of ultra-light axion dark matter in polarimetry data of the European Pulsar Timing Array. Physical Review D, 2025, 111 (6), pp.062005. 10.1103/PhysRevD.111.062005. hal-04845970

HAL Id: hal-04845970

https://hal.science/hal-04845970v1

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Searches for signatures of ultra-light axion dark matter in polarimetry data of the
European Pulsar Timing Array

N. K. Porayko,1,2,* P. Usynina,3,† J. Terol-Calvo,4,5,6,‡ J. Martin Camalich,7,5 G. M. Shaifullah,1,7,8 A. Castillo,9 D. Blas,10,11 L. Guillemot,12,13 M. Peel,14 C. Tiburzi,8 K. Postnov,3,15 M. Kramer,2,16

J. Antoniadis,17,2 S. Babak,18 A.-S. Bak Nielsen,2,19 E. Barausse,20,21 C. G. Bassa,22 C. Blanchard,12,13

M. Bonetti,1,7,23 E. Bortolas,1,7,24 P. R. Brook,25 M. Burgay,8 R. N. Caballero,26 A. Chalumeau,22

D. J. Champion,2 S. Chanlaridis,17 S. Chen,27 I. Cognard,12,13 G. Desvignes,2 M. Falxa,12,18

R. D. Ferdman,28 A. Franchini,1,7 J. R. Gair,29 A. Golden,30 B. Goncharov,31,32 E. Graikou,2

J.-M. Grießmeier,12,13 Y. J. Guo,2 H. Hu,2 F. Iraci,8,33 D. Izquierdo-Villalba,1,7 J. Jang,2

J. Jawor,2 G. H. Janssen,22,34 A. Jessner,2 R. Karuppusamy,2 E. F. Keane,35 M. J. Keith,16

M. A. Krishnakumar,2,19 K. Lackeos,2 K. J. Lee,36,37,38,39 K. Liu,27,2 Y. Liu,19,37 A. G. Lyne,16

J. W. McKee,40,41 R. A. Main,2 M. B. Mickaliger,16 I. C. Nițu,16 A. Parthasarathy,2 B. B. P. Perera,42

D. Perrodin,8 A. Petiteau,43,18 A. Possenti,8 H. Quelquejay Leclere,18 A. Samajdar,44 S. A. Sanidas,16

A. Sesana,1,7,23 L. Speri,29 R. Spiewak,16 B. W. Stappers,16 S. C. Susarla,30 G. Theureau,12,13,45

E. van der Wateren,34,22 A. Vecchio,25 V. Venkatraman Krishnan,2 J. Wang,46 L. Wang,16 and Z. Wu,37

(EPTA Collaboration)

¹Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

²Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

³Sternberg Astronomical Institute, Moscow State University, Universitetsky pr., 13, Moscow 119234, Russia

⁴Instituto de Astrofísica de Canarias, C/ Vía Láctea, s/n E38205 - La Laguna, Tenerife, Spain

⁵Universidad de La Laguna, Dpto. Astrofísica, E38206 - La Laguna, Tenerife, Spain

⁶Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via P. Giuria 1, I-10125 Torino, Italy

⁷INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

⁸INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius (CA), Italy

⁹Departamento de Física, Universidad Antonio Nariño. Cra 3 Este # 45-03, Bogotá, Colombia.

¹⁰Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology,
Campus UAB, 08193 Bellaterra (Barcelona), Spain

¹¹Institució Catalana de Recerca i Estudis Avançats (ICREA),
Passeig Lluís Companys 23, 08010 Barcelona, Spain

¹²LPC2E, OSUC, Univ Orleans, CNRS, CNES, Observatoire de Paris, F-45071 Orleans, France

¹³Observatoire Radioastronomique de Nançay, Observatoire de Paris,
Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

¹⁴Imperial College London, Blackett Lab, Prince Consort Road, London SW7 2AZ, UK

¹⁵Institute for Nuclear Research, Moscow 117312, Russia

¹⁶Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Manchester M13 9PL, UK

¹⁷Institute of Astrophysics, FORTH, N. Plastira 100, 70013, Heraklion, Greece

¹⁸Université Paris Cité CNRS, Astroparticule et Cosmologie, 75013 Paris, France

¹⁹Fakultät für Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany

²⁰SISSA — International School for Advanced Studies,
Via Bonomea 265, 34136, Trieste, Italy and INFN, Sezione di Trieste

²¹IFPU — Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

²²ASTRON, Netherlands Institute for Radio Astronomy,
Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands

²³INAF - Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano, Italy

²⁴INAF - Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio, 5, I-35122 Padova (PD), Italy

²⁵Institute for Gravitational Wave Astronomy and School of Physics and Astronomy,
University of Birmingham, Edgbaston, Birmingham B15 2TT, UK

²⁶Hellenic Open University, School of Science and Technology, 26335 Patras, Greece

²⁷Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, P. R. China

²⁸School of Physics, Faculty of Science, University of East Anglia, Norwich NR4 7TJ, UK

²⁹Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Muühlenberg 1, 14476 Potsdam, Germany

³⁰Ollscoil na Gaillimhe — University of Galway, University Road, Galway, H91 TK33, Ireland

³¹Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy

³²INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy

³³Universitá di Cagliari, Dipartimento di Fisica,
S.P. Monserrato-Sestu Km 0,700 - 09042 Monserrato (CA), Italy

³⁴Department of Astrophysics/IMAPP, Radboud University Nijmegen,
P.O. Box 9010, 6500 GL Nijmegen, The Netherlands
Ultra-light axion-like particles (ALPs) can be a viable solution to the dark matter problem. The scalar field associated with ALPs, coupled to the electromagnetic field, acts as an active birefringent medium, altering the polarisation properties of light through which it propagates. In particular, oscillations of the axionic field induce monochromatic variations of the plane of linearly polarised radiation of astrophysical signals. The radio emission of millisecond pulsars provides an excellent tool to search for such manifestations, given their high fractional linear polarisation and negligible fluctuations of their polarisation properties. We have searched for the evidence of ALPs in the polarimetry measurements of pulsars collected and preprocessed for the European Pulsar Timing Array (EPTA) campaign. Focusing on the twelve brightest sources in linear polarisation, we searched for an astrophysical signal from axions using both frequentist and Bayesian statistical frameworks. For the frequentist analysis, which uses Lomb-Scargle periodograms at its core, no statistically significant signal has been found. The model used for the Bayesian analysis has been adjusted to accommodate multiple deterministic systematics that may be present in the data. A statistically significant signal has been found in the dataset of multiple pulsars with common frequency between $10^{-8}$ Hz and $2 \times 10^{-8}$ Hz, which can most likely be explained by the residual Faraday rotation in the terrestrial ionosphere. Strong bounds on the coupling constant $g_{a\gamma}$, in the same ballpark as other searches, have been obtained in the mass range between $6 \times 10^{-24}$ eV and $5 \times 10^{-21}$ eV. We conclude by discussing problems that can limit the sensitivity of our search for ultra-light axions in the polarimetry data of pulsars, and possible ways to resolve them.

# I. INTRODUCTION

The search for dark matter (DM) continues to be
among the most active research fields in fundamental
physics, astrophysics and cosmology [1, 2]. In recent
years, different models have been scrutinized and their
parameter space constrained by a combination of novel
phenomenological approaches and improvements in the
quantity and quality of data. Regarding the first aspect,
weakly interacting particles with very light masses (of
about $10^{-22}$ eV), which can be non-thermal cosmolog-
ical relics, are playing a special role as candidates for
ultra-light dark matter (ULDM) [3–6]. Irrespective of the
DM problem, such pseudo-scalar Goldstone bosons have
a possible connection to the strong-CP problem in QCD
(the QCD axion) and appear in various Standard Model
extensions of particle physics [7]. They have a distinct

impact on astrophysical observables as compared to mod-
els of DM particles with masses above $\sim$ 1 eV [8]. Very
low masses of ULDM particles imply that their de-Broglie
wave is much larger than their particle separation, which
enables a purely classical description (sometimes referred
to as ‘wave DM’, see [9] for a review and references). As
far as data is concerned, here we will focus on astrophys-
ical observations with an ever-growing capacity to study
new phenomena with increasingly high accuracy. Precise
measurements of pulsars play a privileged role in this
program.

A remarkable experiment probing ultra-light axion-like
particles (ALPs) using pulsars is based on the search
for periodic changes in the time-of-arrivals (TOAs) of
pulses propagating through the axionic field. The con-
densate of the ULDM particles with phenomenologi-
cally interesting masses $m_a \sim 10^{-23} - 10^{-22}$ eV, forms
clumps with a typical size of the de-Broglie wavelength
of $\sim 60$ pc $\times (10^{-22}\text{eV}/m_a)$, where $m_a$ is the axionic
mass. Within each clump, ULDM experiences fast co-
herent oscillation with a typical frequency defined by
$m_a$, inducing variations of the space-time metric with

<footnote n="35"><sup>35</sup>School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland</footnote>

<footnote n="36"><sup>36</sup>Department of Astronomy, School of Physics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="37"><sup>37</sup>National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, P. R. China</footnote>

<footnote n="38"><sup>38</sup>Yunnan Astronomical Observatories, Chinese Academy of Sciences, Kunming 650216, Yunnan, P. R. China</footnote>

<footnote n="39"><sup>39</sup>Beijing Laser Acceleration Innovation Center, Huairou, Beijing 101400, P. R. China</footnote>

<footnote n="40"><sup>40</sup>E. A. Milne Centre for Astrophysics, University of Hull,
Cuttingham Road, Kingston-upon-Hull, HU6 7RX, UK</footnote>

<footnote n="41"><sup>41</sup>Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM),
University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK</footnote>

<footnote n="42"><sup>42</sup>Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA</footnote>

<footnote n="43"><sup>43</sup>IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="44"><sup>44</sup>Institut für Physik und Astronomie, Universität Potsdam,
Haus 28, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany</footnote>

<footnote n="45"><sup>45</sup>Laboratoire Univers et Théories LUTh, Observatoire de Paris,
Université PSL, CNRS, Université de Paris, 92190 Meudon, France</footnote>

<footnote n="46"><sup>46</sup>Ruhr University Bochum, Faculty of Physics and Astronomy,
Astronomical Institute (AIRUB), 44780 Bochum, Germany</footnote>

<footnote>(Dated: December 4, 2024)</footnote>

<footnote>* nataliya.porayko@unimib.it</footnote>

<footnote>† usynina.pg19@physics.msu.ru</footnote>

<footnote>‡ jortecal@protonmail.com</footnote>
twice the frequency. This effect was originally described by [10] and was searched for in the timing data of various Pulsar Timing Arrays (PTAs) [11–16]. The tightest bounds to date using this method have been obtained in [17], constraining the ULDM density to be below a few tenths of the observed DM abundance in the mass range $m_a \sim [10^{-24} \text{ eV}, 10^{-23.3} \text{ eV}]$.

If ULDM couples non-gravitationally to Standard Model particles, it will affect the spin of pulsars as well as the frequency of the atomic clocks used to record pulsar timing measurements, leading to another monochromatic signal in the timing observations at the frequency $\omega = m_a$ [16, 18–20]. In addition to coherent oscillations, ULDM spatial density fluctuations in the Galactic halo source stochastic perturbations of the TOAs at frequencies $\omega \lesssim m_a \sigma^2$ enabling searches for ALPs in the mass range $m_a \sim [10^{-18} \text{ eV}, 10^{-14} \text{ eV}]$ [20].

In the present work, we go beyond the aforementioned searches for ULDM in the timing residuals of PTAs and look for the complementary effect of ALPs emerging in the polarimetry data. In particular, the effective interaction between photons and ALPs, characterised by the dimensionful coupling $g_{a\gamma}$,<fnref n="1" /> induces sinusoidal oscillations at frequency $\omega = m_a$ of the polarisation plane of linearly polarised emission of astrophysical signals. This effect, also known as cosmic birefringence [22, 23], allows to impose constraint on $g_{a\gamma}$ through multiple astrophysical probes, including polarisation measurements of parsec-scale jets in active galaxies [24], searches for birefringence effect in the polarisation of the cosmic microwave background (CMB) [25], the non-observation of a reduction of the net polarised fraction of the CMB light as compared to the $\Lambda$CDM prediction [26]. In addition, two other types of experiments that examine axion production in extreme environments, i.e. solar axion searches [27] and searching for ALPs-to-photon conversion in the supernova SN 1987A [28], enable to set constraints on the coupling constant which are independent of the relic axion abundance (local DM density), in contrast to the other probes mentioned above. Pulsars, which emit highly linearly polarised radiation with unprecedented temporal stability, offer an exciting opportunity to probe and directly detect cosmic birefringence due to ultra-light axions [29]. The latest limits from pulsar polarimetry were obtained using the data collected within the Parkes PTA campaign [21, 30]. Here, based on previous studies, for the first time we investigate the presence of ULDM signatures in the polarisation data of the European PTA (EPTA) covering a total period of 11 years. To this aim, we analyse the polarisation calibrated EPTA dataset,

clean it from known systematics, i.e. from the effects of the terrestrial plasma, and construct the time series of position angles (PAs), which characterise the orientation of the plane of linearly polarised pulsar radiation. Finally, we search for ULDM signals in the PA time series and set upper limits on the coupling constant applying both frequentist and Bayesian data analysis methods.

The paper is organised as follows. In Sec. II, we introduce the concept of ULDM and its effect on pulsar polarisation data. In Sec. III, we introduce the EPTA dataset and discuss the essential steps of data curation. Sec. IV is devoted to the description of the techniques developed to analyse the data, separately by frequentist and Bayesian methods. Our main results are presented in Sec. V. Finally, we summarise our study in Sec. VI. Throughout the paper, we work in the natural system of units $\hbar = c = 1$.

## II. BIREFRINGENCE OF PULSAR LIGHT DUE TO ULTRA-LIGHT AXIONS

A medium formed by an ALP relic<fnref n="2" /> can affect the propagation of electromagnetic (EM) waves by rotating their polarisation plane [32, 33]. This effect arises from the Lagrangian describing the ALP-photon interactions,

$$ \mathcal{L} = -\frac{1}{4} F_{\mu\nu} F^{\mu\nu} + \frac{g_{a\gamma}}{4} a F_{\mu\nu} \tilde{F}^{\mu\nu} + \frac{1}{2} (\partial_{\mu} a \partial^{\mu} a - m_a^2 a^2), \quad (1) $$

where $a$ is the ALP field, $F_{\mu\nu}$ is the EM Maxwell tensor and $\tilde{F}^{\mu\nu}$ is its dual. Note that the coupling constant $g_{a\gamma}$ has units of inverse mass.

The solution to the modified Maxwell equations in the adiabatic (WKB) approximation<fnref n="3" /> results in different angular frequencies for left and right circularly polarised waves

$$ \omega_{\pm} \approx k \pm \frac{1}{2} g_{a\gamma} (\partial_t a + \nabla a \cdot \hat{k}), \quad (2) $$

to lowest order in $g_{a\gamma}$. Therefore, the ALP medium is birefringent and leads to a phase shift between the polarisation components given by

$$ \Delta\phi = \frac{g_{a\gamma}}{2} \int_{t_s}^{t_o} \frac{da}{dt} dt = \frac{g_{a\gamma}}{2} (a_o - a_e) \equiv \frac{g_{a\gamma}}{2} \Delta\phi. \quad (3) $$

In the previous integral $t_o$ refers to the time of observation, while $t_s$ represents the time at emission (source term). This effect is independent of the frequency of the EM wave and depends only on the ALP field values at the source and the observer [32, 33].

<footnote n="1">¹ The photon-ALP interaction in Eq. (1) corresponds to a non-renormalizable operator with a coupling $g_{a\gamma}$ treated within the framework of effective field theories. In this context, the interaction represents the low-energy imprint of new interactions and degrees of freedom that become dynamical at energy scales roughly $\gtrsim 1/g_{a\gamma}$ (see Ref. [21] for further discussion and caveats).</footnote>

<footnote n="2">² Meaning that axions are produced in the early Universe through various mechanism, see, e.g. [31].</footnote>

<footnote n="3">³ Assuming that variations of the ALP field are much longer than the ALP oscillation period $T_a = 2\pi/m_a \sim 4 \times 10^7 \text{ s}(10^{-22} \text{ eV}/m_a)$.</footnote>
Ignoring the EM backreaction term in equations of motion arising from Eq. (1), the ALP field is described by plane waves that at any time $t$ and position $\vec{x}$ is

$$a(\vec{x}, t) = a_0(\vec{x}) \cos(m_a t + \delta(\vec{x})), \quad (4)$$

where $a_0(\vec{x})$ and $\delta(\vec{x})$ are the amplitude and phase of the field, respectively, at $\vec{x}$. This description holds for time intervals shorter than the coherence time $\tau_c = (m_a \sigma^2)^{-1} \simeq 10^{13}$s (10$^{-22}$eV/$m_a$)(200 km/s/$\sigma$)$^2$, where $\sigma$ is the dispersion velocity of dark matter, and for distances smaller than the coherence length $l_c = (m_a \sigma)^{-1}$. Within these limits, the distribution of ULDM behaves as a coherent field oscillating at the nominal Compton frequency

$$\nu = \frac{m_a}{2\pi}. \quad (5)$$

The amplitude of the ALP field is related to the local DM density but its nature is intrinsically stochastic. In a given coherence patch, labelled by $i$, the ALP amplitude is [21, 34, 35],

$$a_{0,i} = \frac{\sqrt{2\rho_i}}{m_a}\alpha_i, \quad (6)$$

where $\rho_i$ is the average DM density in the patch and with $a_{0,i}$ following a Rayleigh distribution.

We now consider the time-dependent change in the polarisation of light emitted by pulsars, induced by the ALP birefringent medium, as measured on Earth [21] (cf. Eq. (3)),

$$\Delta\phi(t) = \frac{g_{a\gamma}}{\sqrt{2}m_a} \left[ \sqrt{\rho_o}\alpha_o \cos(m_a t + \delta_o) - \sqrt{\rho_s}\alpha_s \cos(m_a(t-T) + \delta_s) \right], \quad (7)$$

where $t$ is the local time since the beginning of observations ($t=0$), $T$ is the light travel time from a pulsar to Earth, and $\delta_i$, $\rho_i$, and $\alpha_i$ are the phases, DM densities, and stochastic amplitudes at the source ($i=s$) or observer ($i=0$) patches, respectively.

The distribution of time-averaged inter-pulsar correlation coefficients is given by:

$$P_{ab} = \overline{\Delta\phi(t)_a \Delta\phi(t)_b} = \frac{1}{T_{\text{tot}}} \int_{-T_{\text{tot}}/2}^{T_{\text{tot}}/2} dt \Delta\phi(t)_a \Delta\phi(t)_b, \quad (8)$$

where $a$ and $b$ are pulsar labels and $T_{\text{tot}}$ is the total timespan of observations. In contrast to the case of gravitational waves (GWs) with inter-pulsar correlation pattern described by the Hellings-Downs curve [36], the ALP signal is expected to exhibit monopolar or fractional monopolar correlation, depending on the configuration of the Earth-pulsar systems. In particular, it is easy to show that if $\alpha_0 \gg \alpha_s$, signals for different pulsars are fully monopole-correlated. In Fig. 1 we construct the distribution of PTA-averaged inter-pulsar correlation coefficients $P_{ab}$ over multiple universe realisations. As one can

<figure><img src="image_2.png" /><figcaption>FIG. 1. Distribution of correlation coefficients averaged over multiple pulsar pairs. Different realizations of the Universe were generated using Eq. (7). The parameters $\alpha_0$ and $\alpha_s$ were sampled from a Rayleigh distribution, while $\delta_0$ and $\delta_s$ were assumed to be uniformly distributed from 0 to $2\pi$.</figcaption></figure>

see, there is a significant fraction of universes for which the PTA-averaged correlation coefficients are close to 1. However, there are also a number of universes, for which the signal between different pulsars is completely uncorrelated. In order to accelerate the analysis, we perform incoherent search and neglect any possible inter-pulsar correlation. In [37] authors suggested a method to take into account fractional monopolar correlations. The performance of these two approaches will be compared in future studies.

The right hand side of Eq. (7) can be rearranged in a more compact form:

$$\Delta\phi(t) = \phi_a \cos(m_a t + \varphi_a), \quad (9)$$

where

$$\phi_a = \frac{g_{a\gamma}}{\sqrt{2}m_a} \left( \rho_o \alpha_o^2 + \rho_s \alpha_s^2 - 2\sqrt{\rho_o \rho_s} \alpha_o \alpha_s \cos \chi \right)^{1/2}, \quad (10)$$

and $\chi = m_a T + \delta_o - \delta_s$. Given the large uncertainties on pulsar distances, $m_a T$ spans several cycles of the ALP phase. Hence, hereafter we assume that $\chi$ as well as $\varphi_a$ are distributed uniformly in $[0, 2\pi]$. In addition, we assume that the average value of the DM density at Earth and at the pulsar are the same and equal to the measured local DM density, $\rho_o = \rho_s = 0.4$ GeV cm$^{-3}$ [38–41]. Eq. (9) is the final expression that describes the form of the signal to be searched for in the polarisation measurements of pulsars.

### III. OBSERVATIONS AND DATA REDUCTION

The EPTA Collaboration has been continuously monitoring ~100 millisecond pulsars (MSPs) for over 25 years [42]. The observations are performed using the six largest
telescopes in Europe: the Effelsberg 100-m radio tele-
scope (EFF) in Germany, the 76-m Lovell telescope at
Jodrell Bank Observatory (JBO) in the United Kingdom,
the large radio telescope operated by the Nançay Radio
Observatory (NRT) in France, the 64-m Sardinia Radio
Telescope (SRT) operated by the Italian National Insti-
tute for Astrophysics (INAF) through the Astronomical
Observatory of Cagliari (OAC), and the Westerbork Syn-
thesis Radio Telescope (WSRT) operated by ASTRON,
the Netherlands Institute for Radio Astronomy.

The EPTA experiment is primarily designed to detect
GWs in the nHz frequency band using pulsars as natural
clocks, relying on the technique of high-precision pulsar
timing [43, 44]. Within this technique, TOAs of pulses
from a pulsar are obtained by comparing the observed
total intensity profile to a smoothed template profile, us-
ing a cross-correlation algorithm. Although pulsar tim-
ing traditionally utilises only the Stokes I or total inten-
sity profiles to obtain the TOAs of the pulses<fnref n="4" />, all six
EPTA instruments have recorded full polarimetry data,
such that all four Stokes parameters are available for fur-
ther analysis.

The first data release (DR1) of the EPTA [42] includes
the timing data from 42 MSPs observed between 1996
and 2014. In the summer of 2023, the EPTA collabora-
tion published the second data release (hereafter DR2,
[46]), covering the period from 13.6 to 24.7 years. Due
to the large number of telescopes and backend combina-
tions, as well as the diverse data acquisition and process-
ing schemes across the EPTA observatories, 25 MSPs
were optimally selected among those monitored by the
collaboration, based on the method described in [47].
Prior to TOA calculation, the raw data are folded mod-
ulo the pulsar spin-period with the dspsr software [48].
These folding-mode data cubes undergo a stage of stan-
dard initial processing using the PSRCHIVE software [49],
which includes the removal of radio frequency interfer-
ence (RFI) and polarisation calibration. The latter ac-
counts for the change of parallactic angle during the ob-
servation, the geometry of the feed and possible cross-
coupling between the two feeds of the receptor (see, e.g.
[50, 51]). Missing or inaccurate calibration will inevitably
introduce temporal variations in the polarisation proper-
ties of pulsar light, which will cause unwanted artefacts
and reduce the overall sensitivity of the dataset. For this
reason, an accurate polarimetric calibration is of vital
importance for this and similar studies.

For the current analysis, we make use of calibrated
full-Stokes data, which were pre-analysed as part of the
EPTA DR2 and subsequently utilised to produce TOA
time series. We have limited our analysis to pulsar obser-
vations collected by the radio telescopes with the largest

effective collective area, namely Effelsberg and Nançay.
Our final dataset comprises twelve MSPs with the high-
est signal-to-noise (S/N) in linear polarisation *L*, with
fractional linear polarisation greater than 20. The de-
tails of the observational set-up and the steps that have
been taken as part of the initial pre-analysis are briefly
described in the following paragraphs.

A. Effelsberg data

The 100-m Effelsberg radio telescope is operated by
the Max Planck Institut für Radioastronomie. The ob-
servations for the EPTA project are carried out with
monthly cadence at three widely distributed frequency
bands: 1.4 GHz with P200mm and P217mm receivers,
2.6 GHz with a S110mm receiver and 4.85 GHz with a
S60mm receiver (all summarised in [52]). The specified
integration time per source is 28 min followed by 2 min
of noise diode scan. Since 2011 the ROACH-based (Re-
configurable Open Architecture Computing Hardware)
backend PSRIX [52] digitises the electric signal from
the telescope and applies a polyphase filterbank. The
CoastGuard pipeline, which is part of the Toaster soft-
ware library, handles the subsequent necessary steps of
the raw data processing [52], i.e. the collected pulsar ob-
servations are coherently dedispersed and folded mod-
ulo pulsar period. The resultant data cubes are stored
in the psrfits data format with 10-sec and 1024-bin
resolution in time and pulse phase, respectively. The
RFI-contaminated profiles are identified and erased us-
ing CoastGuard's surgical algorithm. The data are
further calibrated for flux and polarisation. All of the
receivers used for the current study have circularly po-
larised feeds. Polarisation calibration of cleaned data is
performed using the 'SingleAxis' technique described in
[53]. Making use of noise diode scans, the relative gain
and phase difference of the two polarisation channels have
been adjusted. For more details on the initial stages
of raw data processing, see [46]. Polarisation-calibrated
data are available from early 2013 to October 2020 and
have been used in the current study.

B. Nançay data

The Nançay Radio Telescope (NRT) is a meridian
transit-type instrument with a collective area equivalent
to that of a 94-m parabolic dish, making its sensitiv-
ity comparable to that of the Effelsberg radio telescope.
It is equipped with cryogenically cooled receivers, cov-
ering the [1.1, 1.8] GHz and [1.7, 3.5] GHz frequency
ranges. Regular pulsar observations for the EPTA with
the NRT began in 2004. Since August 2011, pulsar ob-
servations with the NRT have been conducted with the
NUPPI backend (for further details, see [54, 55]), which
uses Graphics Processing Units (GPUs) to coherently
dedisperse and fold data from the receiver in real time

<footnote n="4">⁴ The matrix template matching (MTM) technique [45] has been implemented in the PSRCHIVE software to perform TOA extraction using all four Stokes parameters. For instance, MTM has been utilised for TOA extraction for the NUPPI observations.</footnote>
over a total bandwidth of 512 MHz. For the work presented in this article, we considered pulsar observations conducted with the low-frequency receiver at a central frequency of 1484 MHz, which represents the bulk of NRT pulsar observations. The duration of these observations ranges from 20 to 80 mins, and the average observing cadence is ~ 2 days. As with the Effelsberg observations, the frequency channels and time sub-integrations heavily contaminated by RFI have been removed.

Until late 2019, NUPPI pulsar observations were polarisation calibrated with a simple calibration scheme: the calibration was done using observations of a noise diode conducted before each pulsar observation and using the SingleAxis method of PSRCHIVE, which assumes that the polarisation feeds are perfectly orthogonal and that the noise diode signal is 100% linearly polarised and illuminates the two feeds equally and in phase. Since late 2019, observations of bright, highly linearly polarised pulsars have been conducted regularly, approximately monthly, in a special mode where the horn is made to rotate across the observation. These observations, analyzed with the Reception model of PSRCHIVE enable the determination of the full polarimetric response of the NRT. For the work presented here, we used NRT pulsar archives recorded at L-band from October 2019 to April 2023, i.e., calibrated with this much-improved polarisation calibration scheme. Further details about the NRT, pulsar observations with the NRT and the polarisation calibration of the NRT data can be found in [56].

### C. Ionospheric correction

The effect of ULDM is not the only process that can alter the polarisation properties of pulsar emission. In fact, the dominant contribution comes from the Faraday rotation that takes place in the magneto-ionic plasma between the source and the observer, causing the rotation of the plane of polarisation of the linearly polarised part of the radiation. In astronomy, the orientation of the polarisation plane is characterised by the position angle (PA):

$$PA = \Psi_0 + RM \lambda^2, \quad (11)$$

where angle $Ψ_0$ characterises the initial orientation of the plane of polarisation (as emitted by a source) and $λ$ is the observing wavelength. The rate at which the plane of polarisation rotates is characterised by the rotation measure (RM):

$$RM = 0.81 \left[ \frac{\text{rad}}{\text{m}^2} \right] \int_{\text{LoS}} n_e(\mathbf{B} \cdot \text{dl}), \quad (12)$$

where $n_e$ (cm$^{-3}$) is the electron number density of the intervening plasma and the dot product in brackets represents the projection of the magnetic field $\mathbf{B}$ (measured in $\mu$G) onto the infinitesimal interval $dl$ of the distance along the line of sight (LoS) expressed in pc.

As can be seen from Eq. (11), Faraday rotation is strongly dependent on the observational frequency, whereas the effect of ultra-light axions is independent of the EM wave frequency (the signal is achromatic), allowing the two to be disentangled. However, that would only be possible for the observing systems with large fractional bandwidth$^5$, so that the plane of polarisation rotates considerably along the bandwidth. For instance, for the Nançay L-band receiver the relative change in PA between the upper and lower part of the band is only ~ 2 deg, which can not be probed with the current sensitivities. Therefore the two signals due to the ULDM and Faraday rotation in the intervening plasma, have similar manifestation in the data for the given observational setup.

As we are not sensitive to time-independent modifications of the PA (given that the initial polarisation configuration is unknown), only time-variable effects should be considered in the analysis. Generally, there are two known chromatic, i.e. which depend on observational radio frequency, processes that can alter the polarisation properties of pulsar radiation dynamically. The main contributor to variations in RM and hence PA is the terrestrial plasma, i.e. ionosphere and plasmasphere [57–59]. The second type of temporal RM fluctuations is driven by spatial turbulence in the ionised interstellar medium (ISM), over which the LoS scans as a pulsar moves in the tangential plane. However, the latter is estimated to be five orders of magnitude lower than the former and cannot be measured with current sensitivities (see more details in [59]). Therefore, only the atmospheric contribution is accounted for in the current analysis.

The ionosphere is modelled using the conventional single-layer model (SLM), which assumes that the terrestrial plasma is encapsulated in an infinitesimally thin slab fixed at a height of ~ 450 km. Within this approximation, Eq. (12) reduces to a simple product of the LoS projected geomagnetic field $B_E$ and integrated electron density, also known as slant total electron content (STEC). Both values are evaluated at the ionospheric pierce point (IPP), where the LoS intersects the ionospheric layer. Information about the Earth’s magnetic field is provided by the World Magnetic Model [60]. The STEC values are recalculated from the vertical TECs (VTECs), which are electron column densities in the direction of the zenith. VTECs are extracted from the interpolated global ionospheric maps (GIMs), which contain VTEC values defined on a regular temporal and spatial grid reconstructed from the Global Positioning System (GPS) observations. Specifically, for the current work, we have used UQRG GIMs produced by the UPC-IonSat collaboration [61, 62]. The UQRG GIMs were shown to surpass other publicly available GIMs when comparing different

<footnote>$^5$ Fractional bandwidth is the ratio of the total observed bandwidth to the central frequency of operation.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Lomb-Scargle periodogram of the PA residuals for PSR J1600−3053. The spectra of ΔPA, when the ionospheric contribution has been subtracted, before (pre-processing scheme) and after (post-processing scheme) PA extraction are shown in blue (dashed) and turquoise (solid), respectively. The Lomb-Scargle periodogram of ΔPA un-corrected for the ionospheric effect is shown in black (dot-dashed).</figcaption></figure>

ionospheric models using the observed Faraday rotation of pulsars at low frequencies of LOFAR [63]. However, the authors also showed that none of the existing models, including UQRG, have been able to capture the full complexity of ionospheric physics. In particular, it was found that two types of systematics emerge: quasi-sinusoidal variations with a one-year period, and long-term linear trends strongly correlated with the 11-yr Solar cycle. See more details in [63, 64]. The possible presence of this interference was taken into consideration within a Bayesian framework (see Section IV B).

The ionospheric correction was implemented at each fully frequency-resolved pulsar archive (before the stage of PA extraction) by modifying the relevant entry in the psrfits header using --aux_rms command. In this case, the possibly present (weak) effect of depolarisation due to the ionosphere, which can alter the shape of the polarisation profile and bias the PA measurements, is compensated. It is worth noting that we have also tested the performance of the ionospheric correction in the post-processing stage (after PA extraction), neglecting the frequency dependence of the Faraday rotation within an observing band. As expected, for instruments with small fractional bandwidths that we are using in the current project, the difference in the resultant standard deviation of the residuals for the two schemes is negligible (see Fig. 2).

## D. PA extraction

The astrophysical signal of interest manifests itself as a periodic variation of the PA (see Eq. (7)). The following steps are undertaken to capture these changes in PA relative to a baseline. We start by correcting the pulsar data, stored in psrfits format, for the ionospheric contribution, as described in Sec III C. Afterwards, we build a template by summing up a number of observations with a high S/N ratio, followed by the smoothing. We force-align averaged pulse profiles taken at different observational epochs in phase, using the Stokes I template constructed on a previous step as a benchmark. This procedure is necessary as multiple stochastic processes, primarily the intrinsic spin noise, can phase-shift the pulse profiles and bias the final measurements of the PA (as PAs of an individual archive and of a template generally may not match the same phase). Finally, we calculate the difference between the PA across an individual profile and a template. The phase-resolved PA residuals are constructed directly from Stokes Q and U observables, using the following expression:

$$ \Delta PA(\phi) = \frac{1}{2} \left( \arctan \frac{U_{\text{epoch}}}{Q_{\text{epoch}}} - \arctan \frac{U_{\text{tmp1}}}{Q_{\text{tmp1}}} \right) = \frac{1}{2} \arcsin\{U_{\text{epoch}}(\phi)Q_{\text{tmp1}}(\phi) - Q_{\text{epoch}}(\phi)U_{\text{tmp1}}(\phi)\}, \quad (13) $$

where $Q_{\text{epoch}}/U_{\text{epoch}}$ and $Q_{\text{epoch}}/U_{\text{tmp1}}$ are Stokes Q and U profiles as a function of a pulse phase $\phi$ of an individual observation and the template, respectively. The sought ΔPA is a weighted average of ΔPA($\phi$), where we use the inverse square of measurement uncertainties $1/\sigma_{\text{PA}}^2$ as weights:

$$ \sigma_{\text{PA}} = \frac{\sqrt{Q_{\text{epoch}}^2(\phi)\sigma_U^2 + U_{\text{epoch}}^2(\phi)\sigma_Q^2}}{2(Q_{\text{epoch}}^2(\phi) + U_{\text{epoch}}^2(\phi))}. \quad (14) $$

Here, $\sigma_Q$ and $\sigma_U$ are the off-pulse root-mean-squares. The average uncertainty on the ΔPA is calculated as the standard deviation of ΔPA($\phi$) normalised by the number of on-pulse measurements.

The obtained dataset is publicly available on zenodo<fnref n="6" /> and is displayed in Fig. 3. The reported $\chi^2$ values indicate that the uncertainties are considerably underestimated. A similar behaviour has been reported previously by the Parkes PTA Collaboration [PPTA 65]. This bias is caused by the presence of systematics (i.e., deviations from non-Gaussianity) in phase-resolved PA residuals, calculated using Eq. (13). An example of such residual systematics between the epoch-wise and template PA measurements for PSR B1937+21 is shown with the waterfall plot in Fig. 4. We expect that such systematics

<footnote n="6">⁶ https://doi.org/10.5281/zenodo.10.5281/zenodo.1400595</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 3. The PA residuals of the subset of twelve pre-selected EPTA DR2 MSPs. Multiple colours correspond to four different observing systems: P217-Asterix (red circles), P200-Asterix (orange stars), S110-Asterix (green filled points), S60-Asterix (turquoise triangles), 1484-NUPPI (blue points). The shadow area, which extends from MJD 59460 to 59560, represents the malfunction of one of the local oscillators at the Nançay observatory. One of the consequences of this problem are the jumps in PA residuals, which are particularly apparent in the data of PSR J1012+5307, J1643–1224 and J1600–3053. The χ² values indicated at the panels are calculated after eliminating the epochs affected by the malfunction issue.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

FIG. 4. Intensity and polarisation configurations of PSR B1937+21. On the left: template intensity profile (upper panel) and individual PA($\phi$) profiles in black overlapped with the template PA profile in yellow (lower panel). On the right: waterfall plots of the $\Delta$PA($\phi$) observed with NRT. Each line corresponds to a specific epoch of observation. The red line on the left panel is the PA swing of the smoothed template. For the sake of presentation, only forty epochs have been plotted. It is clear that $\Delta$PA($\phi$) exhibits non-Gaussian systematic behaviour. These observed features are particularly pronounced in the regions close to the orthogonal jump of the PA.

may be caused by the scattering in the IISM, or errors in polarisation calibration, or other instrument-related factors, and will be further explored in future work (see Sec. VI for discussion).

## IV. METHODS

In this section, we introduce the methods used to perform parameter estimation and search for the ULDM signature in the data using frequentist and Bayesian inferences. For both types of analyses, we assume that the Earth and pulsar terms fall in different coherent patches and that the monochromatic, i.e. sinusoidal with period unambiguously defined by the mass of the ALP, signal from ULDM is uncorrelated between different pulsars. For our sample of pulsars, this assumption holds for $\ge 5 \times 10^{-23}$ eV. In App. A we further discuss the justification for this assumption and investigate how upper limits are modified for different pulsar-Earth configurations. We also treat the noise in the PA residuals as white Gaussian. In the Bayesian framework, we additionally take measures to account for the underestimation in

the uncertainties and presence of long-term systematics with white noise modifiers, an extra error in quadrature (EQUAD), which represents an additional unaccounted noise, and an error scaling factor (EFAC) [66], and inclusion of linear and quadratic trends in the model, respectively.

### A. Frequentist analysis

The frequentist approach involves searching for periodic signals in unevenly distributed time series with the Lomb-Scargle (LS) periodogram [67, 68], based on a combination of signal analysis techniques [69]. In this work we employ the *generalised* LS periodogram (GLSP) [69, 70] $P_{LS}(\nu)$ and follow the procedure described in [21].

First, we look for peaks in $P_{LS}(\nu)$ and quantify their significance by comparing them with the so-called false alarm probability (FAP). This value quantifies the probability that a significant peak emerges from purely random fluctuations, and we estimate it using the *bootstrap* method, which resamples the time series while keeping
the temporal coordinates fixed using Monte Carlo simulations [69]. The highest peaks of these simulations form a probability distribution that is used directly to compute the FAP.

After verifying that there are no peaks compatible with the expected DM signal, we perform $10^4$ Monte Carlo simulations, injecting a harmonic signal into the time series,

$$ \Delta\phi_{\text{sim}} = \phi \cos(2\pi\nu t + \varphi), \quad (15) $$

with a given frequency $\nu$, amplitude $\phi$, and a random phase $\varphi$ uniformly sampled from [0, 2π]. Using the GLSP, we compute the probability density function (PDF) of $P_{LS}(\nu, \phi)$ and determine $\phi_{95}$ such that $P_{LS}(\nu, \phi_{95})$ matches the experimental $P_{LS}(\nu)$. This provides an upper bound on $\phi$ at 95% confidence level (C.L.). We use a rolling mean method to smooth out the fluctuations, averaging $\phi_{95}$ within a window centred at each frequency. This process is repeated for all frequencies and pulsars.

Finally, we perform a combined analysis, conducting a Monte Carlo computation similar to one in the previous section. Using the PDF of $\phi_a$ obtained earlier, we generate $10^5$ pseudo-experiments, drawing values of $\phi_{a,i}$ for each pulsar, along with values for $\alpha_o$, $\alpha_{s,i}$, and $\chi_i$ for the pulsars. Using Eq. (15), we calculate $g_{a\gamma,i}$ for each pulsar. The global value of $g_{a\gamma}$ is obtained as a weighted mean, with weights determined by propagating the root squares of the variances of $\phi_{a,i}$, directly related to the scale factor of the Rayleigh distribution. By running multiple pseudo-experiments, we construct its PDF and report the upper limit on the axion-photon coupling at the 95% confidence level from the resulting global $g_{a\gamma}$ distribution.

## B. Bayesian analysis

In the Bayesian framework we assume that the $\Delta PA$ measured at observational epoch $t_i$, $\Delta PA^{\text{obs}}(t_i)$, for a single pulsar is a sum of an astrophysical signal of interest $\delta\phi^{\text{ULDM}}(t_i)$, given by Eq. (9), and a Gaussian white noise $n(t_i)$. Given the flexibility of the Bayesian scheme in comparison to the GLSP, we complement the frequentist approach by considering the possible presence of multiple deterministic processes $PA^{\text{det}}(t_i)$ which may emerge due to, e.g. polarisation calibration errors or poor ionospheric modelling:

$$ \Delta PA^{\text{obs}}(t_i) = \delta\phi^{\text{ULDM}}(\mathbf{\Theta}, t_i) + PA^{\text{det}}(\mathbf{\Psi}, t_i) + n(t_i), \quad (16) $$

where $\mathbf{\Theta}$ is a parameter vector of the ALP model and $\mathbf{\Psi}$ parameterise the deterministic contribution. We assume that $PA^{\text{det}}$ linearly depends on $\mathbf{\Psi}$, so that it can be expressed as:

$$ PA^{\text{det}} = M\Psi, \quad (17) $$

where $M$ is the $(N \times m)$ design matrix of partial derivatives of each parameter characterising deterministic trends ($N$ is the total number of measurements and $m$ is the dimensionality of $\mathbf{\Psi}$). In our work, we consider linear and quadratic trends, as well as sinusoidal variations with 1-yr and 11-yr periods responsible for commonly known biases in the ionospheric modelling [59, 71]. The inclusion of the deterministic term in the model removes possible residual trends unrelated to the ALPs signal, but also provides more conservative estimates of the target parameter $g_{a\gamma}$. Within the Bayesian framework, all systems, noises and signals of interest are considered simultaneously using the multivariate Gaussian likelihood:

$$ \mathcal{L}(\Delta\mathbf{PA}|\mathbf{\Theta}, \mathbf{\Psi}) = \frac{1}{\sqrt{|\det 2\pi\mathbf{C}|}} \times \exp \left\{ -\frac{1}{2}(\Delta\mathbf{PA}^{\text{obs}} - \delta\mathbf{\phi}^{\text{ULDM}} - \mathbf{PA}^{\text{det}})^T \mathbf{C}^{-1} (\Delta\mathbf{PA}^{\text{obs}} - \delta\mathbf{\phi}^{\text{ULDM}} - \mathbf{PA}^{\text{det}}) \right\}, \quad (18) $$

where $\mathbf{C} = \sigma_{\text{PA}}^2 \delta_{ij}$ is the covariance matrix of the stochastic processes and $\delta_{ij}$ is the Kronecker delta. Given the misestimation of the uncertainties, two additional parameters, EFAC and EQUAD, were added to the analysis, so that the diagonal elements of the covariance matrix are modified as follows: $C_{ii} = \text{EFAC}^2 \sigma_{\text{PA},ii}^2 + \text{EQUAD}^2$. This likelihood function can be marginalised analytically over the parameters $\mathbf{\Psi}$ and further reduced to a more compact form [72, 73]:

$$ \mathcal{L}(\Delta\mathbf{PA}|\mathbf{\Theta}) = \frac{\exp(-\frac{1}{2}\delta\mathbf{PA}^T\mathbf{G}(\mathbf{G}^T\mathbf{C}\mathbf{G})^{-1}\mathbf{G}^T\delta\mathbf{PA})}{\sqrt{(2\pi)^{N-m}\det(\mathbf{G}^T\mathbf{C}\mathbf{G})}}, \quad (19) $$

where $\delta\mathbf{PA} = \Delta\mathbf{PA} - \delta\mathbf{\phi}^{\text{ULDM}}$ and the $(N \times N - m)$ matrix $\mathbf{G}$ was obtained from the singular-value decomposition of the design matrix $M$. Specifically, $M = \text{UDV}^T$, where $\mathbf{U}$ and $\mathbf{V}$ are $(N \times N)$ and $(m \times m)$ orthogonal matrices, respectively and $\mathbf{D}$ is an $(N \times m)$ diagonal matrix with singular values of $M$. The $\mathbf{G}$ matrix is formed from the $\mathbf{U}$ matrix such that $\mathbf{U} = (\mathbf{G}_c \ \mathbf{G})$. See more details in, e.g. [74].

After constructing the likelihood function, the Bayesian analysis proceeds with the determination of the posterior distributions of the parameters $p^{\text{post}}(\mathbf{\Theta}|\Delta\mathbf{PA})$ by updating the prior distribution of those parameters $p^{\text{prior}}(\mathbf{\Theta})$ through the likelihood function $\mathcal{L}(\Delta\mathbf{PA}|\mathbf{\Theta})$:

$$ p^{\text{post}}(\mathbf{\Theta}|\Delta\mathbf{PA}) = \frac{\mathcal{L}(\Delta\mathbf{PA}|\mathbf{\Theta})p^{\text{prior}}(\mathbf{\Theta})}{Z}, \quad (20) $$

where the denominator $Z$ is called evidence and can be expressed as:

$$ Z = \int \mathcal{L}(\Delta\mathbf{PA}|\mathbf{\Theta}) p^{\text{prior}}(\mathbf{\Theta}) d^m(\mathbf{\Theta}). \quad (21) $$
Evidence is a key criterion for the hypothesis testing (or model selection) problem. One can use it to calculate the Bayes factor (BF) to choose between two hypotheses, which in our case are: i) signal present in the data; and ii) signal missing in the data. For the task of hypothesis testing, we do not evaluate the above integrals, but use the Savage-Dickey ratio instead, which is defined as the ratio of the posterior to the prior density at the critical value $\Theta_0$, corresponding to the zero amplitude $g_{a\gamma} = 0$:

$$B_{21} = \frac{p^{\text{post}}(\Theta_0 | \Delta \text{PA})}{p^{\text{prior}}(\Theta_0)} \quad (22)$$

Conventionally, the Bayes factors $B_{21} = 1$ and $B_{21} = 2$ are considered substantial and evidence, respectively. For the parameter estimation problem, evidence $Z$ plays the role of normalisation constant. In the case of flat uninformative priors, the Eq. (20) is reduced to $p^{\text{post}}(\Theta | \Delta \text{PA}) \sim \mathcal{L}(\Delta \text{PA} | \Theta)$, which is an expression used to construct the posterior distribution of the signal parameters in practice. The model parameters $\Theta$ are split into two categories: i) common to all the pulsars, which are the coupling constant $g_{a\gamma}$, and the axion mass $m_a$; ii) individual to a pulsar, which are the phase of oscillation $\phi_a$, the normalisation parameter $B = \sqrt{\alpha_o^2 + \alpha_s^2 - 2\alpha_o\alpha_s \cos \chi}$, the error in quadrature EQUADs (one per observing system), the error multiplier EFACs (one per observing system) and the OFFSETs (one per observing system), which are introduced to describe the jumps between different observing systems. The distribution of the parameter $B$ has been estimated numerically using the priors of the relevant input parameters ($\alpha_0, \alpha_s$ and $\chi$) and approximated using a kernel density estimator (KDE) (see Fig. 5). The whole set of parameters and the corresponding priors are summarised in Table I.

To reduce the complexity of the current analysis, we neglect any possible cross-correlation between datasets of different pulsars (see Fig. 1). In this case the likelihood in Eq. (19) can be factorised as:

$$\mathcal{L}(\Delta \text{PA} | g_{a\gamma}, m_a, \mathbf{A}_1, \mathbf{A}_2, \dots, \mathbf{A}_N) = \prod_{k=1}^{N} \mathcal{L}(\Delta \text{PA} | g_{a\gamma}, m_a, \mathbf{A}_k), \quad (23)$$

where $\mathbf{A}_k$ are groups of parameters individual to a pulsar. Using this expression, one can perform an efficient search for an astrophysical signal in a reduced parameter space for each pulsar dataset separately and construct the joint posterior probability in the after-processing. Eq. (23) is the final expression that we used to perform the Bayesian run.

## V. RESULTS

We now turn to the description of the main results of the analysis. We present the results of the frequentist

<figure><img src="image_6.png" /><figcaption>FIG. 5. Numerically reconstructed prior distribution of the normalisation parameter $B$ and its KDE approximation (in blue).</figcaption></figure>

and the Bayesian analysis separately, discussing first the detection and then upper limit estimations.

### A. Frequentist analysis

Following the procedure described in Sec. IV A, we first search for peaks in the GLSP. We only find a peak below the 1% FAP, in PSR J1600-3053, at a frequency of $\sim 3 \times 10^{-8}$ Hz corresponding to a periodicity of a year. The lack of similar peaks in the GLSP of all the other pulsars points to a different explanation for the ULDM case. In fact, the 1-year frequency aligns well with a seasonal ionospheric signal, still remaining after the corrections, see Sec. III C. After confirming that there is no signal compatible with the ULDM, we proceed to set constraints on periodic signals for each pulsar, finding the maximum 95% C.L. amplitude, $\phi_{95}$. We show them in Fig. 6.

We see that all pulsars set an upper limit of around a degree in amplitude, approximately constant over a wide range of frequencies, except for the lower frequencies, where the sensitivity drops. In order to translate these results to ULDM parameters, i.e. to the coupling $g_{a\gamma}$ in Eq. (1), we perform the global analysis. We first asses the local DM density of each source and, to take into account stochastic effects, the distance between all of them to check if there are sources in the same coherence patch. We also produce the results under the assumptions that the local DM density at each source and the Earth is the same, $\rho_{DM} \approx 0.4$ GeV/cm$^3$ and neglecting stochastic effects, i.e. setting the stochastic parameter in Eq. (6), $\alpha = 1$. Given that the DM local density is very similar in all cases, the effect is small, and in Fig. 7 we plot the case with $\rho = 0.4$ GeV/cm$^3$, taking into account stochastic effects. See Appendix A for more details on the impact of including stochastic effects. The best sensitivity is found
<table><caption>TABLE I. The set of parameters, used for the Bayesian analysis</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th><th>Comments</th></tr></thead><tbody><tr><td colspan="4"><em>Noise parameters</em></td></tr><tr><td>EFAC</td><td>White-noise modifier per backend</td><td>Uniform [0, 4]</td><td>One parameter per backend</td></tr><tr><td>EQUAD</td><td>Additive white noise per backend</td><td>log<sub>10</sub>-Uniform [-8, 3]</td><td>One parameter per backend</td></tr><tr><td>OFFSET</td><td>An offset between observing systems</td><td>Uniform [-2, 2]</td><td>One parameter per backend</td></tr><tr><td colspan="4"><em>Signal parameters</em></td></tr><tr><td>g<sub>aγ</sub></td><td>photon-axion coupling constant</td><td>log<sub>10</sub>-Uniform [-6, 5] (search)<br>Linear-Exp [-7, 6] (upper limits)</td><td>One parameter per PTA</td></tr><tr><td>f(Hz)</td><td>Oscillation frequency</td><td>log<sub>10</sub>-Uniform [1/T, 10<sup>-6</sup>]</td><td>Fixed, regular grid in log-scale</td></tr><tr><td>φ<sub>a</sub></td><td>Oscillation phase</td><td>Uniform[0, 2π]</td><td></td></tr><tr><td>B</td><td>B = √(α<sub>o</sub><sup>2</sup> + α<sub>s</sub><sup>2</sup> - 2α<sub>o</sub>α<sub>s</sub> cos χ)</td><td>Estimated numerically using KDE<br>(see Fig. 5)</td><td></td></tr><tr><td colspan="4"><em>Additional parameters</em></td></tr><tr><td>α<sub>0</sub></td><td>Stochastic amplitude on Earth</td><td>Rayleigh distribution</td><td>p(α) = αexp(-α<sub>o</sub><sup>2</sup>/2)</td></tr><tr><td>α<sub>s</sub></td><td>Stochastic amplitude on a pulsar</td><td>Rayleigh distribution</td><td>p(α) = αexp(-α<sub>s</sub><sup>2</sup>/2)</td></tr><tr><td>χ</td><td>Phase</td><td>Uniform[0, 2π]</td><td>χ = m<sub>a</sub>T + δ<sub>o</sub> - δ<sub>s</sub></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Frequentist limits on the amplitude at the 95% C.L. ($\phi_{95}$) for the 12 pre-selected EPTA pulsars. These constraints will be used to produce a combined bound directly on $g_{a\gamma}$, which should also include the effects of stochastic amplitudes (so, is not illustrated in the current Fig.).</figcaption></figure>

for the smallest mass considered, $m_a \approx 8 \times 10^{-24}$ eV, where we set an upper limit of $g_{a\gamma}^{95\%} \approx 4 \times 10^{-14}$ GeV$^{-1}$.

## B. Bayesian analysis

To tackle the problem in the Bayesian framework, we first split the frequency range of interest [$10^{-9}$ Hz, $10^{-6}$ Hz] into a regular grid of 48 bins that are evenly distributed in logarithmic scale. We perform hypothesis testing and parameter estimation in each bin separately, while fixing the oscillation frequency $f$. As a result, the only varying parameter that is common for all the pulsars

in Eq. (23) is the coupling constant $g_{a\gamma}$.

The Markov chain Monte Carlo sampler “emcee” [75] is used to estimate the posterior distribution of the model parameters for each pulsar separately. The joint posterior distribution of $g_{a\gamma}$ is reconstructed by taking a product of these individual pulsar posteriors (see Eq. (23)). The obtained distributions are used to both estimate the BF in the Savage-Dickey approximation (see Eq. (22)) and set upper limits on the coupling constant. We choose two different kinds of priors for $g_{a\gamma}$: i) uniform in linear scale for model comparison; and ii) log$_{10}$-uniform for parameter estimation. The obtained BFs are shown in Fig. 8. For most frequency bins log$_{10}$(BF) is less than 1, indicating that no common harmonic signal was found at these particular frequencies. There are six bins between $10^{-8}$ Hz and $2 \times 10^{-8}$ Hz with log$_{10}$(BF) larger than 3, meaning that the hypothesis that a signal with a period of roughly 2 yrs is present in the data is supported. This frequency corresponds to a boson mass of $m_a \sim 5 \times 10^{-23}$ eV. In Fig. 9, we show the factorised likelihood for one of the bins with $f = 1.4 \times 10^{-8}$ Hz, clearly showing non-flat posteriors for multiple pulsars. Notably, a similar bump, but at half the frequency has been detected in [76–78], when searching for continuous GWs and ULDM signatures in the latest EPTA DR2 and North American Nanohertz Observatory for Gravitational Waves datasets, respectively. Although such a signal could arise for different reasons, the most viable is the residual Faraday rotation in the terrestrial plasma which is not sufficiently well-modelled by the ionospheric model, or the effect of the Solar wind.

The upper limits on $g_{a\gamma}^{95\%}$, which are defined as 95%-quantile of the posterior distribution, are shown in the left panel of Fig. 10. For conservative analysis, within which a fit for nuance deterministic parameters is taken into account, the most stringent bounds are obtained at
<figure><img src="image_8.png" /><figcaption>FIG. 7. Upper limits on the photon-ULDM coupling, $g_{a\gamma}$, obtained with the frequentist (dashed gray) and Bayesian (solid gray) for the conservative approach and dotted gray for the optimistic analysis) analysis over a range of frequencies and masses. For the optimistic bayesian and frequentist analyses we cut below a minimum frequency $\nu_{min} = 1/T \approx 2.4 \times 10^{-9}$ Hz, which corresponds to demand at least a full oscillation cycle in our data. Other constraints obtained in previous studies are also shown: the horizontal blue line marks the start of the upper limit set by the helioscope experiment CAST [27], likewise the red horizontal line shows the limit set by the non-observation of a gamma-ray signal from SN1987A produced by axions [28]. The area in cyan (marked with long dashes) shows the parameter space excluded by the search of the birefringence effect in AGN polarisation observations by MOJAVE VLBA [24], in purple (dotted) we show the area excluded by the BICEP-Keck search for the birefringence effect in CMB polarisation [25], the orange area (dotted dashed) is excluded by the same search in polarisation data from PPTA and QUIJOTE [21] and lastly the green area (dashed) is excluded by the non-observation of a 'washout effect' in the polarisation of the CMB produced at last scattering using Planck [26].</figcaption></figure>

$m_a \simeq 3 \times 10^{-23}$ eV, constraining $g_{a\gamma}$ to be below $2 \times 10^{-13}$ GeV$^{-1}$ with 95% confidence. In addition, we compute optimistic upper bounds, where we exclude a fit for deterministic systematics $\text{PA}^{\text{det}}$ from the model. The latter can be directly compared with the frequentist results. These upper limits are shown in the right panel of Fig. 10 and, as expected, provide more stringent constraints at lower frequencies.

The validation of the Bayesian pipeline as well as the method of factorised likelihood is demonstrated in App. B, where a simulated signal is injected in the real dataset and successfully recovered.

In Fig. 7 we show the final results of our analysis for both frequentist and Bayesian analysis for the complete twelve-pulsar dataset. Our bounds are overplotted with the results from other campaigns aimed at searching for ULDM signals in other astrophysical observables.

## VI. CONCLUSIONS

In this paper, we analysed the polarisation data of the latest EPTA DR2 with the purpose of searching for signatures of pseudo-scalar ULDM (axion-like particle ALP) coupled to the EM field. In order to reduce the complexity of our analysis, we focused on full Stokes observations of twelve pulsars gathered with the two largest radio telescopes in Europe, Effelsberg and Nançay. The presence of ULDM, interacting with the EM field with the coupling described by Eq. (1), results in sinusoidal variations of the PA residuals, that we searched for with both frequentist and Bayesian techniques. Another special feature of the effect is its achromaticity, i.e. it does not depend on observing frequency so that it can be effectively decoupled in the analysis from the Faraday rotation caused by the magnetised plasma.

In the frequentist framework, the search and parameter estimation have been performed using GLSP. No ev-
<figure><img src="image_9.png" /><figcaption>FIG. 8. Log₁₀(BF) computed using Savage-Dickey ratio from the joint posterior probabilities reconstructed from the individual posteriors distributions using the method of factorised likelihood. The red line indicates a threshold value of log₁₀(BF)=1.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. Factorised likelihood analysis on the log₁₀(gₐγ) from the EPTA dataset at *f* = 1.4 × 10⁻⁸ Hz. Individual posterior probabilities for different pulsars are highlighted with colored (dashed and dotted dashed) histograms. The joint posterior probability is in black.</figcaption></figure>

idence for ULDM was found in the data, and we have
obtained stringent upper limits across a wide range of
ALP masses. Our best sensitivity is achieved at m_a =
8 × 10⁻²⁴ eV, where we get the 95%-upper limit on the
coupling constant g_aγ⁹⁵% = 5 × 10⁻¹³ GeV⁻¹.

For the Bayesian analysis, the model was extended to
include additional fits for linear trends, and 11-yr and
1-yr periodic signals that may be present in the data due
to polarisation calibration issues and residual ionospheric
systematics. We have performed a frequency-resolved
analysis by splitting the frequency band of interest into
multiple segments and performing search and parameter
estimation independently in each of them. We find strong

evidence for a monochromatic signal in 6 bins in the fre-
quency range between 10⁻⁸ Hz and 2 × 10⁻⁸ Hz. The
corresponding coupling constant of the detected signal is
log₁₀(gₐγ/GeV⁻¹) = -13.0 ± 0.86, which is in marginal
contradiction with Planck and PPTA-QUIJOTE con-
straints (see Fig. 7). The nature of this signal is still
unclear, but it is most likely related to poorly modelled
Faraday rotation in the terrestrial plasma.

The most stringent constraints on the coupling con-
stant in the Bayesian framework are obtained at m_a =
3 × 10⁻²³ eV, where g_aγ⁹⁵% = 2 × 10⁻¹³ GeV⁻¹. Addi-
tionally, we have calculated Bayesian upper limits for a
simplified model (excluding deterministic systematics),
which are in good agreement (within a factor of 2) with
the frequentist results.

There are a number of reasons that can cause the differ-
ence between the two types of analysis. As was discussed
earlier, in contrast to frequentist analysis, Bayesian mod-
els include white noise modifiers, EFAC and EQUAD,
that alter the uncertainties of the obtained model param-
eters, in particular, upper limits on g<sub>aγ</sub>. Moreover, the
GLSP analysis is constructed under the assumption that
the stochastic processes are white. As it was shown in
[79], the presence of coloured noises can introduce severe
errors in the GLSP parameter estimation and detection
statistics. In case of the Bayesian framework, correlated
systematics have been partially accounted for by intro-
ducing long-term deterministic trends in the model (as a
proxy for the ”red noise”).

Some interesting technical issues relating to the search
methods and the nature of the dataset have also emerged
from this analysis. They can be summarised as follows:

* The subtraction of an averaged PA template from individual PA profiles does not result in white Gaussian residuals, leading to a significant underestimation of the uncertainties on the resulting epoch-averaged PAs. A non-negligible change in the PA profile shapes is the main reason for this behaviour. The nature of these variations is uncertain but is most likely related to temporal variations in the IISM scattering, which can cause significant distortion of the PA swing (see [80]). Imperfect calibration, e.g. caused by the residual leakages between two orthogonal receivers, could provide another possible explanation. Similar underestimated uncertainties of the PA residuals and PA profile temporal variations were obtained with other radio instruments reported in [81–83]. The characteristics of the systematics that appear in the phase-resolved PA and possible ways to reduce them will be investigated in future studies.

* The Faraday rotation in the terrestrial ionospheric plasma causes a dynamical shift of the PA and can also be responsible for the distortion of the PA swing described above. Despite of being chromatic, the Faraday rotation can still be responsible for the emergence of ULDM-like signals and subse-
<figure><img src="image_11.png" /><figcaption>FIG. 10. Bayesian upper limits on the coupling constant $g_{a\gamma}$ as a function of frequency (boson mass). The different colours show the upper limits obtained for each pulsar separately, while the joint analysis is shown in black. Conservative upper limits for which the full deterministic model is taken into account are displayed on the left, while more optimistic upper bounds, for which the deterministic terms are excluded from the analysis, are shown in the right panel.</figcaption></figure>

quent false detection. For instance, we suspect that the excess of power that was detected in the current analysis may have been the result of the ionosphere not being modelled correctly. The ionospheric contribution can be accounted for with the external models of the ionosphere, as done in the current analysis. The shortcoming of this scheme is that none of the existing models can accurately describe the complex physics of the magneto-ionic content of the Earth's atmosphere. This problem can be overcome by performing a similar study at lower observing frequencies, in order to break the degeneracy between the chromatic Faraday rotation and the achromatic ULDM signal. In this case, the data can be *self-calibrated* for varying Faraday rotation (both, in the ionospheric and IISM) without the use of third-party models.

* For both the frequentist and Bayesian analysis, we have neglected any possible inter-pulsar correlations caused by the ULDM signal and effectively looked for an incoherent monochromatic signal in different pulsars. It is important to note that the astrophysical signal of interest may be partially monopolar-correlated. One possible way of including cross-correlations in the analysis was considered in [37]. Determining whether the detected signal exhibits any angular correlation can be performed in post-processing (similar to the optimal statistics analysis of the GW background [84]) and is going to be addressed in subsequent work.

* Our upper limits on the coupling constant can be ameliorated by selecting pulsars located in regions of high DM density. Millisecond pulsars that are located close to the Galactic centre and are known to have stable profiles will provide a sig-

nificant increase in sensitivity towards ULDM signals. In particular, long-term observations of the highly linearly polarised PSR J1744−2946, which was recently discovered in a Galactic centre filament, would be an excellent candidate to test the ULDM hypothesis [85].

The EPTA collaboration is an ongoing effort to detect GWs through performing regular observations of MSPs over several decades. The same data can be used to study a wide range of astrophysical problems beyond the detection of GWs. As it was shown in the current paper, searches for ultra-light DM candidates is one of the good examples. Here, we restrict our analysis to the data from Effelsberg and Nançay radio telescopes only. In the future, we plan to extend our analysis to include the data from three more observatories, namely JBO, SRT, and WSRT, as well as the LOFAR interferometer. LOFAR has recently been incorporated into the EPTA collaboration, and will provide low-frequency support to the campaign. Upcoming data releases with longer baselines and updated instrumentation are expected to further increase our sensitivity to ALPs, maximising the potential for a detection of ULDM.

## VII. DATA AVAILABILITY STATEMENT

The data that support the findings of this article are openly available [86], embargo periods may apply. Additional data products used in the current manuscript can be provided upon reasonable request.
ACKNOWLEDGMENTS

The European Pulsar Timing Array (EPTA) is a collaboration between European and partner institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), the University of Manchester (UK), the University of Birmingham (UK), the University of East Anglia (UK), the University of Bielefeld (GER), the University of Paris (FRA), the University of Milan-Bicocca (IT) and Peking University (CHN), with the aim to provide high precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. An Advanced Grant of the European Research Council to implement the Large European Array for Pulsars (LEAP) has also provided funding. The EPTA is part of the International Pulsar Timing Array (IPTA); we would like to thank our IPTA colleagues for their help with this paper.

Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK's Science and Technology Facilities Council (STFC). ICN is also supported by the STFC doctoral training grant ST/T506291/1. The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from "Programme National de Cosmologie and Galaxies" (PNCG), and "Programme National Hautes Energies" (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as a National Facility by the National Institute for Astrophysics (INAF).

The work is supported by the National SKA programme of China (2020SKA0120100), Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. This work is also supported as part of the "LEGACY" MPG-CAS collaboration on low-frequency gravitational wave astronomy. JA acknowledges support from the European Commission (Grant Agreement number: 101094354), the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of "Science and Society – Action Always strive for excellence – "Theodoros Papazoglou" (Project

Number: 01431). AF, AS, DIV, GS, MB acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant "Binary Massive Black Hole Astrophysics" (B Massive, Grant Agreement: 818691). AC and AP acknowledge financial support from the European Research Council (ERC) starting grant 'GIGA' (grant agreement number: 101116134). GD, RK and MK acknowledge support from European Research Council (ERC) Synergy Grant "BlackHoleCam" Grant Agreement Number 610058 and ERC Advanced Grant "LEAP" Grant Agreement Number 337062. JWMK is a CITA Postdoctoral Fellow: This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), (funding reference CITA 490888-16). KP acknowledges support from the State project "Science" by the Ministry of Science and Higher Education of Russia under the contract 075-15-2024-541. AV acknowledges the support of the Royal Society and Wolfson Foundation, and the UK Science and Technology Facilities Council (STFC). JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039) and by the NSF through AccelNet award #2114721. NKP funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projectnummer PO 2758/1–1, through the Walter–Benjamin programme. ASamajdar thanks the Alexander von Humboldt foundation in Germany for a Humboldt fellowship for postdoctoral researchers. AP, DP and MB acknowledge support from the research grant "iPeska" (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016 (Italy). RNC acknowledges financial support from the Special Account for Research Funds of the Hellenic Open University (ELKE-HOU) under the research programme "GRAVPUL" (grant agreement 319/10-10-2022). EvdW, CGB and GHJ acknowledge support from the Dutch National Science Agenda, NWA Startimpuls – 400.17.608.

The work of JTC is supported by the Ministerio de Ciencia e Innovación under FPI contract PRE2019-089992 of the SEV-2015-0548 grant, by the project "Theoretical Astroparticle Physics (TAsP)" funded by the INFN and from the research grant "Addressing systematic uncertainties in searches for dark matter No. 2022F2843L" funded by MIUR. JMC and JTC acknowledge support from the MICINN through the grant "DarkMaps" PID2022-142142NB-I00. D. Blas acknowledges the support from the Departament de Recerca i Universitats from Generalitat de Catalunya to the Grup de Recerca 00649 (Codi: 2021 SGR 00649). EB acknowledges support from the European Union's Horizon Europe programme under the Marie Sklodowska Curie grant agreement No 101105915 (TESIFA), the European Consortium for Astroparticle Theory in the form of an Exchange Travel Grant, the European Union's Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158), and the European Research Coun-
cil (ERC) under the European Union's Horizon 2020 research and innovation program ERC2018-CoG under-grant agreement N. 818691 (B Massive). The research leading to these results has received funding from the Spanish Ministry of Science and Innovation (PID2020-115845GB-100/AEI/10.13039/501100011033). IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Finally JMC and DB acknowledge funding received from the European Union through the grant "UNDARK" of the Widening participation and spreading excellence programme (project number 101159929).

We made extensive use of publicly available ` numpy` [89], ` scipy` [90], ` matplotlib` [91], ` emcee` [75] python libraries as well as ` PSRCHIVE` pulsar software [49].

## Appendix A: Influence of coherence patch size

As stated in the main text, the connection between the axionic field amplitude and the dark matter density takes a stochastic nature, see Eq. (6), that follows a Rayleigh distribution in each coherence domain [34, 35]. With the coherence length

$$l_c = (m_a \sigma)^{-1} \approx 65 \left( \frac{m_a}{10^{-22} \text{ eV}} \right)^{-1} \left( \frac{\sigma}{10^{-3}} \right)^{-1} \text{pc, (A1)}$$

we can determine whether two or more pulsars are in the same coherence patch for a given ULDM mass. For masses below ~5 × 10⁻²³ eV, there are several pulsars that lie in the same coherence patch. Then, to calculate the DM local density in each source and the Earth we have used a Navarro-Frenk-White (NFW) [98, 99] profile with the parameters $r_s = 14.48$ kpc and $\rho_s = 0.566$ GeV/cm³ [100]. To see the effect of assuming a common DM density, $\rho_{DM} \approx 0.4$ GeV/cm³, and a common coherence patch we performed the global frequentist analysis in four possible cases, as shown in Fig. 11. We labelled the analysis as deterministic when we neglect stochastic effects, i.e. $\alpha_0 = \alpha_p = 1$, and DM when the local DM density assuming NFW profile has been implemented for each pulsar. We can see that the result does not change significantly. For the main result presented in Fig. 7, we choose the stochastic case with equal DM density.

## Appendix B: Validation of the Bayesian search pipeline

To validate our methods and check the robustness of our Bayesian schemes, we have performed an injection-recovery test. The simulated ULDM signal was injected

into the real pulsar polarimetry data used in the current study. In order to recover the signal, we run a fully Bayesian analysis as described in Sec IV B. The posterior distribution of astrophysical parameters for

<figure><img src="image_12.png" /><figcaption>FIG. 11. Frequentist upper limits on the coupling constant $g_{a\gamma}$ produced in four different regimes of 'coherence'. See text on details of these modes.</figcaption></figure>

PSR J1744−1134 are shown in Fig. 12. The joint posterior distribution was constructed by multiplying individual 2D-posterior distributions ($g_{a\gamma}, f$) of 12 pulsars. The factorised distribution of $g_{a\gamma}$ marginalised over $f$ is shown in Fig. 13, demonstrating the successful recovery of the injected value.

## Appendix C: Observation of the residual systematics in the current analysis

After ionospheric corrections, see Sec. III C, some spurious signals from the modelling might still be left. In this section, we highlight some of these signatures that were detected in our dataset within both frequentist and Bayesian frameworks. For such a demonstration, we removed the $PA^{det}$ from the Bayesian model. For example, in PSR J1600−3053 we found a signal compatible with one year periodicity (see the left side of Figs. 14 and 15). For the case of PSR J1643−1224, in the Bayesian analysis a signal with a frequency of (11 yr)⁻¹ emerged. However, in the frequentist analysis, while there is a peak with relatively high power at that frequency, it is not statistically significant, i.e. it is not above the FAP level (see the right side of Figs. 14 and 15). Both signals are well-known residual quasi-periodicities after the subtraction of the SLM ionospheric model corresponding to a seasonal effect (~1 year) and to the solar cycle (~11 years).

<footnote>[1] G. Bertone and D. Hooper, History of dark matter, Reviews of Modern Physics 90, 045002 (2018), arXiv:1605.04909 [astro-ph.CO].</footnote>
<figure><img src="image_13.png" /><figcaption>FIG. 12. Posterior distribution of astrophysical model parameters which are summarised in Table I. The simulated signal has been injected into the real data of PSR J1744-1134. The injected values are shown with blue solid lines and the contours are 1-, 2- and 3-$\sigma$ intervals.</figcaption></figure>

[2] A. Arbey and F. Mahmoudi, Dark matter and the early Universe: A review, Progress in Particle and Nuclear Physics 119, 103865 (2021), arXiv:2104.11488 [hep-ph].

[3] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the invisible axion, Physics Letters B 120, 127 (1983).

[4] L. F. Abbott and P. Sikivie, A cosmological bound on the invisible axion, Physics Letters B 120, 133 (1983).

[5] M. Dine and W. Fischler, The not-so-harmless axion, Physics Letters B 120, 137 (1983).

[6] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO].

[7] J. Jaeckel and A. Ringwald, The Low-Energy Frontier of Particle Physics, Annual Review of Nuclear and Particle Science 60, 405 (2010), arXiv:1002.0329 [hep-ph].

[8] D. J. Marsh, Axion cosmology, Physics Reports 643, 1 (2016), axion cosmology.

[9] L. Hui, Wave Dark Matter, ARA&A 59, 247 (2021), arXiv:2101.11735 [astro-ph.CO].

[10] A. Khmelnitsky and V. Rubakov, Pulsar timing signal from ultralight scalar dark matter, J. Cosmology Astropart. Phys. 2014, 019 (2014), arXiv:1309.5888 [astro-ph.CO].

[11] N. K. Porayko and K. A. Postnov, Constraints on ultralight scalar dark matter from pulsar timing, Phys. Rev. D 90, 062008 (2014), arXiv:1408.4670 [astro-ph.CO].

[12] N. K. Porayko, X. Zhu, Y. Levin, L. Hui, G. Hobbs, A. Grudskaya, K. Postnov, M. Bailes, N. D. R. Bhat, W. Coles, S. Dai, J. Dempsey, M. J. Keith, M. Kerr, M. Kramer, P. D. Lasky, R. N. Manchester, S. Oslowski, A. Parthasarathy, V. Ravi, D. J. Reardon, P. A. Rosado, C. J. Russell, R. M. Shannon, R. Spiewak, W. van Straten, L. Toomey, J. Wang, L. Wen, X. You, and PPTA Collaboration, Parkes Pulsar Timing Ar-
<figure><img src="image_14.png" /><figcaption>FIG. 13. Factorised likelihood analysis on $g_{\alpha\gamma}$. The posterior distributions of individual pulsars are shown in colour (dashed and dotted dashed), while the joint posterior is in black. The injected value is highlighted in blue. One should note that the plot shows the marginalised distribution of a 2-D ($g_{\alpha\gamma}, m_\alpha$) posterior.</figcaption></figure>

ray constraints on ultralight scalar-field dark matter, Phys. Rev. D **98**, 102002 (2018), arXiv:1810.03227 [astro-ph.CO].

[13] R. Kato and J. Soda, Search for ultralight scalar dark matter with NANOGrav pulsar timing arrays, J. Cosmology Astropart. Phys. **2020**, 036 (2020), arXiv:1904.09143 [astro-ph.HE].

[14] K. Nomura, A. Ito, and J. Soda, Pulsar timing residual induced by ultralight vector dark matter, European Physical Journal C **80**, 419 (2020), arXiv:1912.10210 [gr-qc].

[15] X. Xue, Z.-Q. Xia, X. Zhu, Y. Zhao, J. Shu, Q. Yuan, N. D. R. Bhat, A. D. Cameron, S. Dai, Y. Feng, B. Goncharov, G. Hobbs, E. Howard, R. N. Manchester, A. Parthasarathy, D. J. Reardon, C. J. Russell, R. M. Shannon, R. Spiewak, N. Thyagarajan, J. Wang, L. Zhang, S. Zhang, and PPTA Collaboration, High-precision search for dark photon dark matter with the Parkes Pulsar Timing Array, Physical Review Research 4, L012022 (2022), arXiv:2112.07687 [hep-ph].

[16] A. Afzal, G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, P. T. Baker, B. Bécsy, J. J. Blanco-Pillado, L. Blecha, K. K. Boddy, A. Brazier, P. R. Brook, S. Burke-Spolaor, R. Burnette, R. Case, M. Charisi, S. Chatterjee, K. Chatziioannou, B. D. Cheeseboro, S. Chen, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, C. J. Cutler, M. E. Decesar, D. Degan, P. B. Demorest, H. Deng, T. Dolch, B. Drachler, R. von Eckardstein, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, N. Garver-Daniels, P. A. Gentile, K. A. Gersbach, J. Glaser, D. C. Good, L. Guertin, K. Gültekin, J. S. Hazboun, S. Hourihane, K. Islo, R. J. Jennings, A. D. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, V. S. H. Lee, N. Lewandowska, R. R. Lino Dos Santos, T. B. Littenberg, T. Liu, D. R. Lorimer, J. Luo, R. S. Lynch, C.-P.

Ma, D. R. Madison, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, P. M. Meyers, C. M. F. Mingarelli, A. Mitridate, J. Nay, P. Natarajan, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, T. T. Pennucci, B. B. P. Perera, P. Petrov, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, S. C. Sardesai, A. Schmiedekamp, C. Schmiedekamp, K. Schmitz, T. Schröder, L. Schult, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, P. Stratmann, J. P. Sun, A. Susobhanan, J. K. Swiggum, J. Taylor, S. R. Taylor, T. Trickle, J. E. Turner, C. Unal, M. Vallisneri, S. Verma, S. J. Vigeland, H. M. Wahl, Q. Wang, C. A. Witt, D. Wright, O. Young, K. M. Zurek, and Nanograv Collaboration, The NANOGrav 15 yr Data Set: Search for Signals from New Physics, ApJ **951**, L11 (2023), arXiv:2306.16219 [astro-ph.HE].

[17] C. Smarra, B. Goncharov, E. Barausse, J. Antoniadis, S. Babak, A. S. B. Nielsen, C. G. Bassa, A. Berthereau, M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay, R. N. Caballero, A. Chalumeau, D. J. Champion, S. Chanlaridis, S. Chen, I. Cognard, G. Desvignes, M. Falxa, R. D. Ferdman, A. Franchini, J. R. Gair, E. Graikou, J. M. Grießmeier, L. Guillemot, Y. J. Guo, H. Hu, F. Iraci, D. Izquierdo-Villalba, J. Jang, J. Jawor, G. H. Janssen, A. Jessner, R. Karuppusamy, E. F. Keane, M. J. Keith, M. Kramer, M. A. Krishnakumar, K. Lackeos, K. J. Lee, K. Liu, Y. Liu, A. G. Lyne, J. W. McKee, R. A. Main, M. B. Mickaliger, I. C. NiÅŁu, A. Parthasarathy, B. B. P. Perera, D. Perrodin, A. Peiteau, N. K. Porayko, A. Possenti, H. Q. Leclere, A. Samajdar, S. A. Sanidas, A. Sesana, G. Shaifullah, L. Speri, R. Spiewak, B. W. Stappers, S. C. Susarla, G. Theureau, C. Tiburzi, E. van der Wateren, A. Vecchio, V. V. Krishnan, J. Wang, L. Wang, Z. Wu, and European Pulsar Timing Array, Second Data Release from the European Pulsar Timing Array: Challenging the Ultralight Dark Matter Paradigm, Phys. Rev. Lett. **131**, 171001 (2023), arXiv:2306.16228 [astro-ph.HE].

[18] P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, and W. A. Terrano, Dark matter direct detection with accelerometers, Phys. Rev. D **93**, 075029 (2016), arXiv:1512.06165 [hep-ph].

[19] D. E. Kaplan, A. Mitridate, and T. Trickle, Constraining fundamental constant variations from ultralight dark matter with pulsar timing arrays, Phys. Rev. D **106**, 035032 (2022), arXiv:2205.06817 [hep-ph].

[20] C. Smarra et al., Constraints on conformal ultralight dark matter couplings from the European Pulsar Timing Array, Phys. Rev. D **110**, 043033 (2024), arXiv:2405.01633 [astro-ph.HE].

[21] A. Castillo, J. Martin-Camalich, J. Terol-Calvo, D. Blas, A. Caputo, R. T. G. Santos, L. Sberna, M. Peel, and J. A. Rubiño Martín, Searching for dark-matter waves with PPTA and QUIJOTE pulsar polarimetry, JCAP **06** (06), 014, arXiv:2201.03422 [astro-ph.CO].

[22] S. M. Carroll, G. B. Field, and R. Jackiw, Limits on a Lorentz- and parity-violating modification of electrodynamics, Phys. Rev. D **41**, 1231 (1990).

[23] S. M. Carroll and G. B. Field, Einstein equivalence principle and the polarization of radio galaxies, Phys. Rev. D **43**, 3789 (1991).

[24] M. M. Ivanov, Y. Y. Kovalev, M. L. Lister, A. G. Panin, A. B. Pushkarev, T. Savolainen, and S. V. Troitsky,
<figure><img src="image_15.png" /></figure>

FIG. 14. Log₁₀(BF) computed using the Savage-Dickey ratio from the individual posterior distributions of PSR J1600−3053 (on the left) and PSR J1643−1224 (on the right). The red line indicates a threshold value of log₁₀(BF)=1. The calculated BFs exhibit peaks at 1-yr and 11-yr.

<figure><img src="image_16.png" /></figure>

FIG. 15. Generalised Lomb-Scargle periodogram and the corresponding FAP at 1%, 5% and 32% of PSR J1600−3053 (on the left) and PSR J1643−1224 (on the right). We observe a peak at 1-yr in PSR J1600−3053 above the FAP level.

Constraining the photon coupling of ultra-light dark-matter axion-like particles by polarization variations of parsec-scale jets in active galaxies, JCAP 02, 059, arXiv:1811.10997 [astro-ph.CO].

[25] P. A. R. Ade et al. (BICEP/Keck), BICEP/Keck XIV: Improved constraints on axionlike polarization oscillations in the cosmic microwave background, Phys. Rev. D 105, 022006 (2022), arXiv:2108.03316 [astro-ph.CO].

[26] M. A. Fedderke, P. W. Graham, and S. Rajendran, Axion Dark Matter Detection with CMB Polarization, Phys. Rev. D 100, 015040 (2019), arXiv:1903.02666 [astro-ph.CO].

[27] C. Collaboration, New cast limit on the axion–photon interaction, Nature Physics 13, 584–590 (2017).

[28] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, and A. Ringwald, Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles, JCAP 02, 006, arXiv:1410.3747 [astro-ph.HE].

[29] T. Liu, G. Smoot, and Y. Zhao, Detecting axionlike dark matter with linearly polarized pulsar light, Phys. Rev. D 101, 063012 (2020), arXiv:1901.10981 [astro-ph.CO].

[30] A. Caputo, L. Sberna, M. Frías, D. Blas, P. Pani, L. Shao, and W. Yan, Constraints on millicharged dark matter and axionlike particles from timing of radio waves, Phys. Rev. D 100, 063515 (2019), arXiv:1902.02695 [astro-ph.CO].

[31] P. Arias, N. Bernal, J. K. Osiński, and L. Roszkowski, Dark matter axions in the early universe with a period of increasing temperature, J. Cosmology Astropart. Phys. 2023, 028 (2023), arXiv:2207.07677 [hep-ph].

[32] S. M. Carroll, G. B. Field, and R. Jackiw, Limits on a Lorentz and Parity Violating Modification of Electrodynamics, Phys. Rev. D 41, 1231 (1990).

[33] D. Harari and P. Sikivie, Effects of a Nambu-Goldstone boson on the polarization of radio galaxies and the cosmic microwave background, Phys. Lett. B 289, 67
(1992).

[34] J. W. Foster, N. L. Rodd, and B. R. Safdi, Reveal-
ing the Dark Matter Halo with Axion Direct Detec-
tion, Phys. Rev. D 97, 123006 (2018), arXiv:1711.10489
[astro-ph.CO].

[35] G. P. Centers et al., Stochastic fluctuations of bosonic dark matter, *Nature Commun.* **12**, 7321 (2021), arXiv:1905.13650 [astro-ph.CO].

[36] R. W. Hellings and G. S. Downs, Upper limits on the isotropic gravitational radiation background from pulsar timing analysis., ApJ 265, L39 (1983).

[37] T. Liu, X. Lou, and J. Ren, Pulsar polarization arrays, Phys. Rev. Lett. 130, 121401 (2023).

[38] J. Bovy and S. Tremaine, On the Local Dark Matter Density, ApJ 756, 89 (2012), arXiv:1205.4033 [astro-ph.GA].

[39] J. I. Read, The local dark matter density, *Journal of Physics G Nuclear Physics* **41**, 063101 (2014), arXiv:1404.1938 [astro-ph.GA].

[40] S. Sivertsson, H. Silverwood, J. I. Read, G. Bertone, and P. Steger, The local dark matter density from SDSS-SEGUE G-dwarfs, MNRAS 478, 1677 (2018), arXiv:1708.07836 [astro-ph.GA].

[41] P. F. de Salas, Dark matter local density determina-
tion based on recent observations, in *Journal of Physics
Conference Series*, Journal of Physics Conference Se-
ries, Vol. 1468 (IOP, 2020) p. 012020, arXiv:1910.14366
[astro-ph.GA].

[42] G. Desvignes, R. N. Caballero, L. Lentati, J. P. W. Verbiest, D. J. Champion, B. W. Stappers, G. H. Janssen, P. Lazarus, S. Oślowski, S. Babak, C. G. Bassa, P. Brem, M. Burgay, I. Cognard, J. R. Gair, E. Graikou, L. Guillemot, J. W. T. Hessels, A. Jessner, C. Jordan, R. Karuppusamy, M. Kramer, A. Lassus, K. Lazaridis, K. J. Lee, K. Liu, A. G. Lyne, J. McKee, C. M. F. Mingarelli, D. Perrodin, A. Petiteau, A. Possenti, M. B. Purver, P. A. Rosado, S. Sanidas, A. Sesana, G. Shaifullah, R. Smits, S. R. Taylor, G. Theureau, C. Tiburzi, R. van Haasteren, and A. Vecchio, High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array, MNRAS 458, 3341 (2016), arXiv:1602.08511 [astro-ph.HE].

[43] M. V. Sazhin, Opportunities for detecting ultralong gravitational waves, Soviet Ast. **22**, 36 (1978).

[44] S. Detweiler, Pulsar timing measurements and the search for gravitational waves, ApJ 234, 1100 (1979).

[45] W. van Straten, Radio Astronomical Polarimetry and High-Precision Pulsar Timing, ApJ 642, 1004 (2006), arXiv:astro-ph/0510334 [astro-ph].

[46] EPTA Collaboration, J. Antoniadis, S. Babak, A. S. Bak Nielsen, C. G. Bassa, A. Berthereau, M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay, R. N. Caballero, A. Chalumeau, D. J. Champion, S. Chanlaridis, S. Chen, I. Cognard, G. Desvignes, M. Falxa, R. D. Ferdman, A. Franchini, J. R. Gair, B. Goncharov, E. Graikou, J. M. Grießmeier, L. Guillemot, Y. J. Guo, H. Hu, F. Iraci, D. Izquierdo-Villalba, J. Jang, J. Jawor, G. H. Janssen, A. Jessner, R. Karuppusamy, E. F. Keane, M. J. Keith, M. Kramer, M. A. Krishnakumar, K. Lackeos, K. J. Lee, K. Liu, Y. Liu, A. G. Lyne, J. W. McKee, R. A. Main, M. B. Mickaliger, I. C. Nițu, A. Parthasarathy, B. B. P. Perera, D. Perrodin, A. Petiteau, N. K. Porayko, A. Possenti, H. Quelquejay Leclere, A. Samajdar, S. A. Sanidas, A. Sesana,

G. Shaifullah, L. Speri, R. Spiewak, B. W. Stappers,
S. C. Susarla, G. Theureau, C. Tiburzi, E. van der Wa-
teren, A. Vecchio, V. Venkatraman Krishnan, J. P. W.
Verbiest, J. Wang, L. Wang, and Z. Wu, The second
data release from the European Pulsar Timing Array. I.
The dataset and timing analysis, A&A 678, A48 (2023),
arXiv:2306.16224 [astro-ph.HE].

[47] L. Speri, N. K. Porayko, M. Falxa, S. Chen, J. R. Gair, A. Sesana, and S. R. Taylor, Quality over quantity: Optimizing pulsar timing array analysis for stochastic and continuous gravitational wave signals, MNRAS 518, 1802 (2023), arXiv:2211.03201 [astro-ph.HE].

[48] W. van Straten and M. Bailes, DSPSR: Digital Signal Processing Software for Pulsar Astronomy, PASA 28, 1 (2011), arXiv:1008.3973 [astro-ph.IM].

[49] A. W. Hotan, W. van Straten, and R. N. Manchester, PSRCHIVE and PSRFITS: An Open Approach to Radio Pulsar Data Storage and Analysis, PASA 21, 302 (2004), arXiv:astro-ph/0404549 [astro-ph].

[50] C. Heiles, P. Perillat, M. Nolan, D. Lorimer, R. Bhat, T. Ghosh, M. Lewis, K. O'Neil, C. Salter, and S. Stanimirovic, Mueller Matrix Parameters for Radio Telescopes and Their Observational Determination, PASP 113, 1274 (2001), arXiv:astro-ph/0107352 [astro-ph].

[51] D. R. Lorimer and M. Kramer, *Handbook of Pulsar Astronomy* (Cambridge University Press, 2012).

[52] P. Lazarus, R. Karuppusamy, E. Graikou, R. N. Caballero, D. J. Champion, K. J. Lee, J. P. W. Verbiest, and M. Kramer, Prospects for high-precision pulsar timing with the new Effelsberg PSRIX backend, MNRAS 458, 868 (2016), arXiv:1601.06194 [astro-ph.IM].

[53] M. C. Britton, Radio Astronomical Polarimetry and the Lorentz Group, ApJ 532, 1240 (2000), arXiv:astro-ph/9911101 [astro-ph].

[54] G. Desvignes, I. Cognard, D. A. Smith, D. Champion, L. Guillemot, M. Kramer, P. Lespagnol, F. Octau, and G. Theureau, The SPAN512 mid-latitude pulsar survey at the Nançay Radio Telescope, A&A 667, A79 (2022), arXiv:2209.01806 [astro-ph.HE].

[55] I. Cognard, G. Theureau, L. Guillemot, K. Liu, A. Lassus, and G. Desvignes, Nançay contribution to the worldwide pulsar programs, in *SF2A-2013: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics*, edited by L. Cambresy, F. Martins, E. Nuss, and A. Palacios (2013) pp. 327–330.

[56] L. Guillemot, I. Cognard, W. van Straten, G. Theureau, and E. Gérard, Improving pulsar polarization and timing measurements with the Nançay Radio Telescope, A&A 678, A79 (2023), arXiv:2308.01599 [astro-ph.HE].

[57] D. Oberoi and C. J. Lonsdale, Media responsible for Faraday rotation: A review, Radio Science 47, RS0K08 (2012).

[58] C. Sotomayor-Beltran, C. Sobey, J. W. T. Hessels, G. de Bruyn, A. Noutsos, A. Alexov, J. Anderson, A. Asgekar, I. M. Avruch, R. Beck, M. E. Bell, M. R. Bell, M. J. Bentum, G. Bernardi, P. Best, L. Birzan, A. Bonafede, F. Breitling, J. Broderick, W. N. Brouw, M. Brüggen, B. Ciardi, F. de Gasperin, R. J. Dettmar, A. van Duin, S. Duscha, J. Eislöffel, H. Falcke, R. A. Fallows, R. Fender, C. Ferrari, W. Frieswijk, M. A. Garrett, J. Grießmeier, T. Grit, A. W. Gunst, T. E. Hassall, G. Heald, M. Hoeft, A. Horneffer, M. Iacobelli, E. Juette, A. Karastergiou, E. Keane, J. Kohler, M. Kramer, V. I. Kondratiev, L. V. E.
Koopmans, M. Kuniyoshi, G. Kuper, J. van Leeuwen, P. Maat, G. Macario, S. Markoff, J. P. McKean, D. D. Mulcahy, H. Munk, E. Orru, H. Paas, M. Pandey-Pommier, M. Pilia, R. Pizzo, A. G. Polatidis, W. Reich, H. Röttgering, M. Serylak, J. Sluman, B. W. Stappers, M. Tagger, Y. Tang, C. Tasse, S. ter Veen, R. Vermeulen, R. J. van Weeren, R. A. M. J. Wijers, S. J. Wijnholds, M. W. Wise, O. Wucknitz, S. Yatawatta, and P. Zarka, Calibrating high-precision Faraday rotation measurements for LOFAR and the next generation of low-frequency radio telescopes, A&A 552, A58 (2013), arXiv:1303.6230 [astro-ph.IM].

[59] N. K. Porayko, A. Noutsos, C. Tiburzi, J. P. W. Verbiest, A. Horneffer, J. Künsemöller, S. Osłowski, M. Kramer, D. H. F. M. Schnitzeler, J. M. Anderson, M. Brüggen, J. M. Grießmeier, M. Hoeft, D. J. Schwarz, M. Serylak, and O. Wucknitz, Testing the accuracy of the ionospheric Faraday rotation corrections through LOFAR observations of bright northern pulsars, MN-RAS 483, 4100 (2019), arXiv:1812.01463 [astro-ph.IM].

[60] NCEI Geomagnetic Modeling Team and British Geological Survey, World Magnetic Model 2020, NOAA National Centers for Environmental Information. doi: 10.25921/11v3-da71 (2020), doi: 10.25921/11v3-da71.

[61] A. J. Mannucci, B. D. Wilson, D. N. Yuan, C. H. Ho, U. J. Lindqwister, and T. F. Runge, A global mapping technique for GPS-derived ionospheric total electron content measurements, Radio Science 33, 565 (1998).

[62] M. Hernández-Pajares, J. M. Juan, J. Sanz, R. Orus, A. Garcia-Rigo, J. Feltens, A. Komjathy, S. C. Schaer, and A. Krankowski, The IGS VTEC maps: a reliable source of ionospheric information since 1998, Journal of Geodesy 83, 263 (2009).

[63] N. K. Porayko, A. Noutsos, C. Tiburzi, J. P. W. Verbiest, A. Horneffer, J. Künsemöller, S. Osłowski, M. Kramer, D. H. F. M. Schnitzeler, J. M. Anderson, M. Brüggen, J. M. Grießmeier, M. Hoeft, D. J. Schwarz, M. Serylak, and O. Wucknitz, Testing the accuracy of the ionospheric Faraday rotation corrections through LOFAR observations of bright northern pulsars, MN-RAS 483, 4100 (2019), arXiv:1812.01463 [astro-ph.IM].

[64] N. K. Porayko, M. Mevius, M. Hernández-Pajares, C. Tiburzi, G. Olivares Pulido, Q. Liu, J. P. W. Verbiest, J. Künsemöller, M. A. Krishnakumar, A.-S. Bak Nielsen, M. Brüggen, V. Graffigna, R.-J. Dettmar, M. Kramer, S. Osłowski, D. J. Schwarz, G. M. Shaifullah, and O. Wucknitz, Validation of global ionospheric models using long-term observations of pulsar Faraday rotation with the LOFAR radio telescope, Journal of Geodesy 97, 116 (2023).

[65] W. M. Yan, R. N. Manchester, W. van Straten, J. E. Reynolds, G. Hobbs, N. Wang, M. Bailes, N. D. R. Bhat, S. Burke-Spolaor, D. J. Champion, W. A. Coles, A. W. Hotan, J. Khoo, S. Osłowski, J. M. Sarkissian, J. P. W. Verbiest, and D. R. B. Yardley, Polarization observations of 20 millisecond pulsars, MNRAS 414, 2087 (2011), arXiv:1102.2274 [astro-ph.SR].

[66] L. Lentati, P. Alexander, M. P. Hobson, F. Feroz, R. van Haasteren, K. J. Lee, and R. M. Shannon, TEM-PONEST: a Bayesian approach to pulsar timing analysis, MNRAS 437, 3004 (2014), arXiv:1310.2120 [astro-ph.IM].

[67] N. R. Lomb, Least-squares frequency analysis of unequally spaced data, Astrophys. Space Sci. 39, 447 (1976).

[68] J. D. Scargle, Studies in astronomical time series analysis. 2. Statistical aspects of spectral analysis of unevenly spaced data, Astrophys. J. 263, 835 (1982).

[69] J. T. VanderPlas, Understanding the Lomb–Scargle Periodogram, Astrophys. J. Suppl. 236, 16 (2018).

[70] M. Zechmeister and M. Kurster, The generalised Lomb-Scargle periodogram. A new formalism for the floating-mean and Keplerian periodograms, Astron. Astrophys. 496, 577 (2009), arXiv:0901.2573 [astro-ph.IM].

[71] N. K. Porayko, M. Mevius, M. Hernández-Pajares, C. Tiburzi, G. Olivares Pulido, Q. Liu, J. P. W. Verbiest, J. Künsemöller, M. A. Krishnakumar, A.-S. Bak Nielsen, M. Brüggen, V. Graffigna, R.-J. Dettmar, M. Kramer, S. Osłowski, D. J. Schwarz, G. M. Shaifullah, and O. Wucknitz, Validation of global ionospheric models using long-term observations of pulsar Faraday rotation with the LOFAR radio telescope, Journal of Geodesy 97, 116 (2023).

[72] R. van Haasteren, Y. Levin, P. McDonald, and T. Lu, On measuring the gravitational-wave background using Pulsar Timing Arrays, MNRAS 395, 1005 (2009), arXiv:0809.0791 [astro-ph].

[73] R. van Haasteren and Y. Levin, Understanding and analysing time-correlated stochastic signals in pulsar timing, MNRAS 428, 1147 (2013), arXiv:1202.5932 [astro-ph.IM].

[74] S. R. Taylor, The Nanohertz Gravitational Wave Astronomer, arXiv e-prints , arXiv:2105.13270 (2021), arXiv:2105.13270 [astro-ph.HE].

[75] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, PASP 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM].

[76] J. Antoniadis, P. Arumugam, S. Arumugam, S. Babak, M. Bagchi, A. S. Bak Nielsen, C. G. Bassa, A. Bathula, A. Berthereau, M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay, R. N. Caballero, A. Chalumeau, D. J. Champion, S. Chanlaridis, S. Chen, I. Cognard, S. Dandapat, D. Deb, S. Desai, G. Desvignes, N. Dhanda-Batra, C. Dwivedi, M. Falxa, I. Ferranti, R. D. Ferdman, A. Franchini, J. R. Gair, B. Goncharov, A. Gopakumar, E. Graikou, J. M. Grießmeier, L. Guillemot, Y. J. Guo, Y. Gupta, S. Hisano, H. Hu, F. Iraci, D. Izquierdo-Villalba, J. Jang, J. Jawor, G. H. Janssen, A. Jessner, B. C. Joshi, F. Kareem, R. Karuppusamy, E. F. Keane, M. J. Keith, D. Kharbanda, T. Kikunaga, N. Kolhe, M. Kramer, M. A. Krishnakumar, K. Lackeos, K. J. Lee, K. Liu, Y. Liu, A. G. Lyne, J. W. McKee, Y. Maan, R. A. Main, S. Manzini, M. B. Mickaliger, I. C. Nitu, K. Nobleson, A. K. Paladi, A. Parthasarathy, B. B. P. Perera, D. Perrodin, A. Petiteau, N. K. Porayko, A. Possenti, T. Prabu, H. Quelquejay Leclere, P. Rana, A. Samajdar, S. A. Sanidas, A. Sesana, G. Shaifullah, J. Singha, L. Speri, R. Spiewak, A. Srivastava, B. W. Stappers, M. Surnis, S. C. Susarla, A. Sobhanan, K. Takahashi, P. Tarafdar, G. Theureau, C. Tiburzi, E. van der Wateren, A. Vecchio, V. Venkatraman Krishnan, J. P. W. Verbiest, J. Wang, L. Wang, and Z. Wu, The second data release from the European Pulsar Timing Array V. Search for continuous gravitational wave signals, arXiv e-prints , arXiv:2306.16226 (2023), arXiv:2306.16226 [astro-ph.HE].

[77] C. Smarra, B. Goncharov, E. Barausse, J. Antoniadis, S. Babak, A. S. B. Nielsen, C. G. Bassa, A. Berthereau,
M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay, R. N. Caballero, A. Chalumeau, D. J. Champion, S. Chanlaridis, S. Chen, I. Cognard, G. Desvignes, M. Falxa, R. D. Ferdman, A. Franchini, J. R. Gair, E. Graikou, J. M. Grießmeier, L. Guillemot, Y. J. Guo, H. Hu, F. Iraci, D. Izquierdo-Villalba, J. Jang, J. Jawor, G. H. Janssen, A. Jessner, R. Karuppusamy, E. F. Keane, M. J. Keith, M. Kramer, M. A. Krishnakumar, K. Lackeos, K. J. Lee, K. Liu, Y. Liu, A. G. Lyne, J. W. McKee, R. A. Main, M. B. Mickaliger, I. C. NiÅŁu, A. Parthasarathy, B. B. P. Perera, D. Perrodin, A. Peiteau, N. K. Porayko, A. Possenti, H. Q. Leclere, A. Samajdar, S. A. Sanidas, A. Sesana, G. Shaifullah, L. Speri, R. Spiewak, B. W. Stappers, S. C. Susarla, G. Theureau, C. Tiburzi, E. van der Wateren, A. Vecchio, V. V. Krishnan, J. Wang, L. Wang, Z. Wu, and European Pulsar Timing Array, Second Data Release from the European Pulsar Timing Array: Challenging the Ultralight Dark Matter Paradigm, Phys. Rev. Lett. 131, 171001 (2023), arXiv:2306.16228 [astro-ph.HE].

[78] G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, P. T. Baker, B. Bécsy, L. Blecha, A. Bräzier, P. R. Brook, S. Burke-Spolaor, R. Case, J. A. Casey-Clyde, M. Charisi, S. Chatterjee, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, M. E. Decesar, P. B. Demorest, M. C. Digman, T. Dolch, B. Drachler, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, N. Garver-Daniels, P. A. Gentile, J. Glaser, D. C. Good, K. Gültekin, J. S. Hazboun, S. Hourihane, R. J. Jennings, A. D. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, N. Lewandowska, T. Liu, D. R. Lorimer, J. Luo, R. S. Lynch, C.-P. Ma, D. R. Madi son, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, P. M. Meyers, C. M. F. Mingarelli, A. Mitridate, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, T. T. Pennucci, B. B. P. Perera, P. Petrov, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, S. C. Sardesai, A. Schmiedekamp, C. Schmiedekamp, K. Schmitz, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, A. Susobhanan, J. K. Swiggum, J. Taylor, S. R. Taylor, J. E. Turner, C. Unal, M. Vallisneri, R. van Haasteren, S. J. Vigeland, H. M. Wahl, C. A. Witt, O. Young, and Nanograv Collaboration, The NANOGrav 15 yr Data Set: Bayesian Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries, ApJ 951, L50 (2023), arXiv:2306.16222 [astro-ph.HE].

[79] W. Coles, G. Hobbs, D. J. Champion, R. N. Manchester, and J. P. W. Verbiest, Pulsar timing analysis in the presence of correlated noise, MNRAS 418, 561 (2011), arXiv:1107.5366 [astro-ph.IM].

[80] A. Karastergiou, The complex polarization angles of radio pulsars: orthogonal jumps and interstellar scattering, MNRAS 392, L60 (2009), arXiv:0901.1826 [astro-ph.SR].

[81] W. M. Yan, R. N. Manchester, G. Hobbs, W. van Straten, J. E. Reynolds, N. Wang, M. Bailes, N. D. R. Bhat, S. Burke-Spolaor, D. J. Champion, A. Chaudhary, W. A. Coles, A. W. Hotan, J. Khoo, S. Osłowski, J. M. Sarkissian, and D. R. B. Yardley, Rotation measure variations for 20 millisecond pulsars, Ap&SS 335, 485 (2011), arXiv:1105.4213 [astro-ph.SR].

[82] S. Osłowski, W. van Straten, P. Demorest, and M. Bailes, Improving the precision of pulsar timing through polarization statistics, MNRAS 430, 416 (2013), arXiv:1301.2374 [astro-ph.IM].

[83] L. Dey, M. A. McLaughlin, H. M. Wahl, P. B. Demorest, Z. Arzoumanian, H. Blumer, P. R. Brook, S. Burke-Spolaor, H. T. Cromartie, M. E. DeCesar, T. Dolch, J. A. Ellis, R. D. Ferdman, E. C. Ferrara, W. Fiore, E. Fonseca, N. Garver-Daniels, P. A. Gentile, J. Glaser, D. C. Good, R. J. Jennings, M. L. Jones, M. T. Lam, D. R. Lorimer, J. Luo, R. S. Lynch, C. Ng, D. J. Nice, T. T. Pennucci, N. S. Pol, S. M. Ransom, R. Spiewak, I. H. Stairs, K. Stovall, and J. K. Swiggum, Exploring pulsar timing precision: A comparative study of polarization calibration methods for NANOGrav data from the Green Bank Telescope, arXiv e-prints, arXiv:2406.13463 (2024), arXiv:2406.13463 [astro-ph.HE].

[84] S. J. Chamberlin, J. D. E. Creighton, X. Siemens, P. Demorest, J. Ellis, L. R. Price, and J. D. Romano, Time-domain implementation of the optimal cross-correlation statistic for stochastic gravitational-wave background searches in pulsar timing data, Phys. Rev. D 91, 044048 (2015), arXiv:1410.8256 [astro-ph.IM].

[85] M. E. Lower, S. Dai, S. Johnston, and E. D. Barr, A Milisecond Pulsar Binary Embedded in a Galactic Center Radio Filament, ApJ 967, L16 (2024), arXiv:2404.09098 [astro-ph.HE].

[86] N. K. Porayko, P. Usynina, J. Terol-Calvo, J. Martin-Camalich, and G. Shaifullah, Searches for signatures of ultra-light axion dark matter in polarimetry data of the European Pulsar Timing Array, 10.5281/zenodo.14005957 (2024).

[87] F. Wilczek, Quantum field theory, Reviews of Modern Physics Supplement 71, S85 (1999), arXiv:hep-th/9803075 [hep-th].

[88] E. G. M. Ferreira, Ultra-light dark matter, A&A Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO].

[89] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Shepard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS].

[90] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: fundamental algorithms for scientific computing in python, Nature Methods 17, 261 (2020), arXiv:1907.10121 [cs.MS].

[91] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Computing in Science and Engineering 9, 90 (2007).

[92] H. Kim and A. Mitridate, Stochastic ultralight dark matter fluctuations in pulsar timing arrays, Phys. Rev. D 109, 055017 (2024), arXiv:2312.12225 [astro-ph.D].
[hep-ph].

[93] M. M. Ivanov, Y. Y. Kovalev, M. L. Lister, A. G. Panin, A. B. Pushkarev, T. Savolainen, and S. V. Troitsky, Constraining the photon coupling of ultra-light dark-matter axion-like particles by polarization variations of parsec-scale jets in active galaxies, J. Cosmology Astropart. Phys. 2019, 059 (2019), arXiv:1811.10997 [astro-ph.CO].

[94] T. Fujita, R. Tazaki, and K. Toma, Hunting Axion Dark Matter with Protoplanetary Disk Polarimetry, Phys. Rev. Lett. 122, 191101 (2019), arXiv:1811.03525 [astro-ph.CO].

[95] M. A. Fedderke, P. W. Graham, and S. Rajendran, Axion dark matter detection with CMB polarization, Phys. Rev. D 100, 015040 (2019), arXiv:1903.02666 [astro-ph.CO].

[96] G. Sigl and P. Trivedi, Axion-like Dark Matter Constraints from CMB Birefringence, arXiv e-prints , arXiv:1811.07873 (2018), arXiv:1811.07873 [astro-ph.CO].

[97] A. Castillo, J. Martin-Camalich, J. Terol-Calvo, D. Blas, A. Caputo, R. T. Génova Santos, L. Sberna, M. Peel, and J. A. Rubiño-Martín, Searching for dark-matter waves with PPTA and QUIJOTE pulsar polarimetry, J. Cosmology Astropart. Phys. 2022, 014 (2022), arXiv:2201.03422 [astro-ph.CO].

[98] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of Cold Dark Matter Halos, ApJ 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph].

[99] J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal Density Profile from Hierarchical Clustering, ApJ 490, 493 (1997), arXiv:astro-ph/9611107 [astro-ph].

[100] M. Cirelli, A. Strumia, and J. Zupan, Dark Matter, arXiv e-prints , arXiv:2406.01705 (2024), arXiv:2406.01705 [hep-ph].

[101] P. Lazarus, R. Karuppusamy, E. Graikou, R. N. Caballero, D. J. Champion, K. J. Lee, J. P. W. Verbiest, and M. Kramer, Prospects for high-precision pulsar timing with the new Effelsberg PSRIX backend, MNRAS 458, 868 (2016), arXiv:1601.06194 [astro-ph.IM].