Constraints on millicharged dark matter and axion-like particles from timing of radio waves

TL;DR

This paper uses pulsar timing and fast radio burst data to set new constraints on millicharged dark matter and axion-like particles, improving our understanding of their possible properties and interactions.

Contribution

It introduces novel astrophysical constraints on millicharged dark matter and axion-like particles using radio wave timing data, expanding the parameter space limits.

Findings

Constraints on millicharged dark matter: ε/mₘᵢₗₗᵢ ≲ 10⁻⁸ eV⁻¹ for mₘᵢₗₗᵢ ≳ 10⁻⁶ eV.

Constraints on axion-like particles: g/mₐ ≲ 10⁻¹³ GeV⁻¹ / (10⁻²² eV).

Current pulsar data constrain new regions of parameter space, with future data expected to improve these bounds.

Abstract

We derive novel constraints on millicharged dark matter and ultralight axion-like particles using pulsar timing and fast radio burst observations. Millicharged dark matter affects the dispersion measure of the time of arrival of radio pulses in a way analogous to free electrons. Light pseudo-scalar dark matter, on the other hand, causes the polarization angle of radio signals to oscillate. We show that current and future data can set strong constraints in both cases. For dark matter particles of charge $\epsilon e$, these constraints are ${\epsilon}/{m_{\rm milli}} \lesssim 10^{-8}{\rm eV}^{-1}$, for masses $m_{\rm milli}\gtrsim 10^{-6}\,$eV. For axion-like particles, the analysis of signals from pulsars yields constraints in the axial coupling of the order of $g/m_a\lesssim 10^{-13} {\rm GeV}^{-1}/(10^{-22}{\rm eV})$. Both bounds scale as $(\rho/\rho_{\rm dm})^{1/2}$ if the energy density $\rho$ of the components is a fraction of the total dark matter energy density $\rho_{\rm dm}$. We do a detailed study of both effects using data from two samples of pulsars in the galaxy and in globular clusters, as well as data from FRB 121102 and PSR J0437$-$4715. We show that in both cases actual pulsar data constrain a new region of the parameter space for these models, and will improve with future pulsar-timing observations.

Figures

Tables

Table 1. Table 1: List of local pulsars considered in this work.

Pulsar	Parallax $\left(\mathrm{mas}\right)$	DM $\left(\mathrm{pc}{\mathrm{cm}}^{-3}\right)$	$n_e\left({\mathrm{cm}}^{-3}\right)$
J1024$-$0719	0.770 $\pm$ 0.23	6.4778 $\pm 0.0006$	0.009036
J1012+5307	0.710 $\pm$ 0.17	9.02314 $\pm$ 0.00007	0.007827
J2010$-$1323	0.300 $\pm$ 0.10	22.177 $\pm$ 0.005	0.004931
J2234+0611	0.700 $\pm$ 0.20	10.7645 $\pm$ 0.0015	0.008292
J1909$-$3744	0.810 $\pm$ 0.03	10.3932 $\pm$ 0.01	0.016935
B2020+28	0.370 $\pm$ 0.12	24.63109 $\pm$ 0.00018	0.012689
B1508+55	0.470 $\pm$ 0.03	19.6191 $\pm$0.0003	0.004691
J2017+0603	0.400 $\pm$ 0.20	23.92344 $\pm$0.00009	0.011004
B1534+12	0.860 $\pm$ 0.18	11.61944 $\pm$ 0.00002	0.009608
J0108$-$1431	4.200 $\pm$ 1.40	2.38 $\pm$ 0.19	0.009246
B0031$-$07	0.930 $\pm$ 0.08	10.922 $\pm$ 0.006	0.008167
J1023+0038	0.731 $\pm$ 0.022	14.325 $\pm$ 0.01	0.008209
B1237+25	1.160 $\pm$ 0.08	9.25159 $\pm$ 0.00053	0.008940

Table 2. Table 2: List of pulsars in globular clusters considered in this work.

Pulsar	Cluster	$d\left(\mathrm{pc}\right)$	DM $\left(\mathrm{pc}{\mathrm{cm}}^{-3}\right)$	$n_e\left({\mathrm{cm}}^{-3}\right)$
B1516+02B	M5	8500	29.47 $\pm 0.11$	0.000027
B1516+02A	M5	8500	30.08 $\pm$ 0.05	0.000027
J1518+0204D	M5	8000	29.3 $\pm$ 0.11	0.000042
J1518+0204E	M5	8000	29.3 $\pm$ 0.11	0.000042
J1518+0204C	M5	8000	29.3146 $\pm$ 0.006	0.000042
J2140-2310A	M30	9200	25.0640 $\pm$ 0.0041	0.000015
J2140-2310B	M30	9200	25.09 $\pm$ 0.12	0.000015
J0024-7204X	47Tuc	4690	24.539 $\pm$ 0.005	0.0083
J0024-7204Z	47Tuc	4690	24.47 $\pm$ 0.01	0.0083
J0024-7204Z	47Tuc	4690	24.29 $\pm$ 0.03	0.0083
B0021-72H	47Tuc	4690	24.37 $\pm$ 0.02	0.0083
B0021-72E	47Tuc	4690	24.236 $\pm$ 0.002	0.0083
J0024-7204R	47Tuc	4690	24.361 $\pm$ 0.007	0.0083

Equations

$$(\square -m_a^2)\varphi$$

$${\partial }_{\mu }{F}^{\mu \nu }$$

$${\cal M}:=\left(\begin{array}[]{ccc}\omega+\Delta_{xx}&0&0\\ 0&\omega+\Delta_{yy}&gB/2\\ 0&gB/2&\omega-m_{a}^{2}/(2\omega)\end{array}\right)\,.$$

$${\Delta }_{xx}^{\mathrm{plasma}}\sim {\Delta }_{yy}^{\mathrm{plasma}}\sim -\frac{{\omega }_p^2}{2\omega },$$

$$k_0$$

$$v_0^{-1}$$

$$v_{-}^{-1}$$

$$\Delta t_{\mathrm{DM}}^{\mathrm{astro}}$$

$$\sim 4.15\left(\frac{{\mathrm{DM}}_{\mathrm{astro}}}{\mathrm{pc}{\mathrm{cm}}^{-3}}\right){\left(\frac{\nu }{\mathrm{GHz}}\right)}^{-2}\mathrm{ms},$$

$$\Delta t_{\mathrm{DM}}^{\mathrm{milli}}=\frac{1}{2\pi }\frac{{ϵ}^2{e}^2}{m_{\mathrm{milli}}}\int \mathrm{d}l{n}_{\mathrm{milli}}\left({\nu }^{-2}-{\nu }_{\infty }^{-2}\right).$$

$${\rm DM{}}_{\rm milli}=\Big{(}\frac{\epsilon}{m_{\rm milli}}\Big{)}^{2}m_{e}\int{\rm d}l\,\rho_{\rm milli}\,,$$

$$\frac{ϵ}{m_{\mathrm{milli}}}\lesssim \frac{10^{-8}}{\mathrm{eV}}\sqrt{\frac{0.3\mathrm{GeV}/{\mathrm{cm}}^3}{{\rho }_{\mathrm{milli}}}}\sqrt{\frac{{\mathrm{DM}}_{\mathrm{obs}}}{20\mathrm{pc}/{\mathrm{cm}}^3}}\sqrt{\frac{400\mathrm{pc}}{d}},$$

$$\ln\mathcal{L}=-\frac{1}{2}\sum_{i=1}^{N}\frac{\Big{(}{\rm DM{}}^{i}_{\rm obs}-{\rm DM{}}^{i}_{\rm astro}-{\rm DM{}}^{i}_{\rm milli}\Big{)}^{2}}{\sigma_{i}^{2}}\,.$$

$$\frac{ϵ}{m_{\mathrm{milli}}}\lesssim \frac{4\times 10^{-9}}{\mathrm{eV}}\sqrt{\frac{0.3\mathrm{GeV}/{\mathrm{cm}}^3}{{\rho }_{\mathrm{milli}}}}\mathrm{at}95\%C.L.$$

$$\varphi (x,t)={\tilde{\varphi }}_0(x)\int {\mathrm{d}}^{3}v{e}^{-\frac{v^2}{{\sigma }_0^2}}{e}^{i({\omega }_{v}t-m_a\vec{v}\cdot \vec{x})+i{\phi }_{\mathbf{v}}}+c.c.,$$

$$\begin{array}{rl}\theta (t,T)\sim & 1.4\times 10^{-2}\sin (m_{a}t+\delta )\\ & \left(\frac{g}{10^{-12}{\mathrm{GeV}}^{-1}}\right)\frac{10^{-22}\mathrm{eV}}{m_a}\mathrm{rad},\end{array}$$

Full text

††preprint: KCL-PH-TH/2019-14

Unraveling the nature of dark matter (DMa) is among the most urgent issues in fundamental physics. Indirect searches aim at detecting the effects of DMa in astrophysical observations, beyond its pure gravitational interaction. Given the feeble interaction of DMa with standard model fields, precise measurements are particularly promising for these searches. When one requires precision, a particular measurement stands out in astrophysics: the time of arrival (TOA) of radio waves from pulsars and fast radio bursts (FRBs). The use of pulsar timing has already been suggested to study the effects of dark matter Khmelnitsky and Rubakov (2014); Porayko and Postnov (2014); Pani (2015); Clark et al. (2016); Blas et al. (2017); Schutz and Liu (2017); De Martino et al. (2017); Caputo et al. (2018); Dror et al. (2019). In this Letter we present new results for DMa models directly coupled to light from the propagation of radio pulses from pulsars and FRBs. A more comprehensive exploration will be presented elsewhere Blas et al. (tion).

If DMa is coupled to the electromagnetic field, one expects modifications in the emission, propagation, and detection of radio pulses. We focus here on the effects during the propagation, which are robust under astrophysical uncertainties. In particular, we derive stringent constraints on millicharged DMa and axion-like particles (ALPs) based on dispersion measurements (DM) of radio signals from pulsars and FRBs, and on the modulation of the light polarization angle due to axion-like DMa in the Milky Way.

We give a unified treatment, where the millicharged DMa and ALPs are considered as independent species. In the former case we consider that (a fraction of) the DMa is made of particles with mass $m_{\rm milli}$ and electric charge $q=\epsilon e$ ($\epsilon\ll 1$) De Rujula et al. (1990); Perl and Lee (1997); Holdom (1986); Sigurdson et al. (2004); Davidson et al. (2000); McDermott et al. (2011); Berlin et al. (2018); Ejlli (2017). As an example, this coupling arises in models where the DMa is charged under a dark photon, which is kinematically coupled to the visible photon McDermott et al. (2011); Berlin et al. (2018). In our analysis we remain agnostic to the origin of this term and other possible model-dependent signatures behind the charge of the DMa, and focus on constraining $\epsilon$. Regarding ALPs, we assume the existence of axion-like Peccei and Quinn (1977); Weinberg (1965); Preskill et al. (1983); Abbott and Sikivie (1983); Dine and Fischler (1983), pseudo-scalar DMa of mass $m_{a}$ (represented by the field $\phi$ below).

The relevant field equations read

[TABLE]

where $g$ is the ALP-photon axial coupling, $j^{\nu}$ is the ordinary electron current, whereas $j^{\nu}_{\rm milli}$ is the current from millicharged particles. The role of this term in the propagation of radio-waves will be studied in the next section, under the assumption of a cold distribution of the millicharge DMa  component.

Dispersion in the TOA. We consider the propagation of a light signal of frequency $\nu=\omega/(2\pi)$ along the $z$ direction in the presence of a homogeneous background magnetic field polarized along (say) the $y$ direction, $\vec{B}=(0,B,0)$. We neglect a possible $B_{z}$ component in this section since its role in dispersion of the light signal in a medium of particles of mass $m_{q}$ and charge $q$ is suppressed by $qB_{z}/(m_{q}\omega)$, always small for the cases we study. For the first part of this work, DMa is considered as a cold-medium with vanishing background values for the fields appearing in (1)-(2). When $\omega\gg m_{a}$, the propagation of the light signal in this medium is described by the first-order system $i\frac{\partial}{\partial z}\left|\psi(z)\right\rangle={\cal M}\left|\psi(z)\right\rangle$, where the $\left|\psi(z)\right\rangle$ is a linear combination of the two photon polarizations along the $x$ and $y$ directions and of the ALP state Raffelt and Stodolsky (1988). The $3\times 3$ mixing matrix reads Dupays et al. (2005)

[TABLE]

The terms $\Delta_{xx}$ and $\Delta_{yy}$ contain both QED vacuum polarization effects and plasma effects Raffelt and Stodolsky (1988); Dupays et al. (2005). The first ones are of order $\Delta_{xx}^{\rm QED}\sim\Delta_{yy}^{\rm QED}\sim\omega\frac{e^{2}}{45\pi}\left(\frac{B}{B_{c}}\right)^{2}$, where $B_{c}\approx 4\times 10^{13}\,{\rm G}$ Adler (1971). We shall only consider interstellar magnetic fields, for which $B\ll B_{c}$ and $\Delta^{\rm QED}$ effects are negligible. Plasma effects arise from the presence of free charges. In the limit where the photon energy is much smaller than the mass of the charged, cold particles Gell-Mann et al. (1954); Goldberger and Watson (1964); Latimer (2013),

[TABLE]

where $\omega_{p}^{2}:=\sum_{i}\frac{4\pi n_{i}q_{i}^{2}}{m_{i}}$ is the plasma frequency for particles with charge $q_{i}$, mass $m_{i}$, and number density $n_{i}$. The normal modes corresponding to (3) satisfy

[TABLE]

with $\Delta_{\omega}=(m_{a}^{2}-\omega_{p}^{2})^{2}+4B^{2}g^{2}\omega^{2}$. The last term in $\Delta_{\omega}$ is always subdominant and we treat it perturbatively.

The TOA of a signal traveling at speed $v=\partial\omega/\partial k$ across a distance $d$ is $T=\int_{0}^{d}\frac{{\rm d}l}{v}=\int_{0}^{d}{\rm d}l\frac{\partial k}{\partial\omega}$ along the line of sight. From the previous expressions one finds for the relevant polarizations,

[TABLE]

In the absence of new physics $(\epsilon=g=0$), the previous modes propagate with velocity $v_{0}$. For a photon with frequency $\nu$, a background of cold free electrons yields a time delay

[TABLE]

relative to a photon with high enough energy ($\nu_{\infty}$ in the previous formula) Lorimer and Kramer (2005). Here ${\rm DM{}}_{\rm astro}:=\int n_{e}{\rm d}l$ is the standard dispersion measure (DM) from electrons with number density $n_{e}$ along the light of sight. The last line is also the observational definition of the dispersion measure, ${\rm DM}_{\rm obs}$. Comparing this number with the ALP-photon coupling term in Eq. (7) one sees that the modifications from the interstellar or intergalactic magnetic fields ($B\lesssim\mu$G) are only relevant for $g>\rm GeV^{-1}$, which is already excluded by other methods, e.g. Anastassopoulos et al. (2017). We ignore these terms in the following. We have checked that the high magnetic field of the pulsar magnetosphere is also not relevant for our studies and we ignore it. Finally, the local conditions of FRBs are not known. It is rather unlikely that they play a role in the DM and even more that they cancel the effects from the DMa plasma, Eqs. (6) and (7). We hence restrict our analysis of the TOA to the millicharged DMa.

TOA constraints on millicharged DMa. As we explained above, we now focus on the case of millicharged DMa, i.e. $g=0$. The contribution of the millicharged DMa to the time delay is given by an an expression analogue to (8), now considering the DMa  particles as the dispersive medium

[TABLE]

In this case the observed DM is dominated by the sum of the contributions from ordinary electrons and millicharged particles (see also Gardner and Latimer (2010)), ${\rm DM{}}_{\rm obs}={\rm DM{}}_{\rm astro}+{\rm DM{}}_{\rm milli}$, where the millicharged contribution is obtained by comparing (8) and (9),

[TABLE]

where $\rho_{\rm milli}$ is the density of millicharged particles, which is equal to or smaller than the full DMa density $\rho_{\rm dm}$. While the effect of ${\rm DM{}}_{\rm astro}$ and ${\rm DM{}}_{\rm milli}$ are completely degenerate, for a source at a distance $d$ any measurement of the DM can be translated into a conservative upper bound on ${\epsilon}/{m_{\rm milli}}$ by simply requiring that all the DM  is due to DMa, i.e. ${\rm DM{}}_{\rm milli}<{\rm DM{}}_{\rm obs}$. This yields

[TABLE]

where we normalized the quantities by typical values within the galaxy. This estimate gives already a rather stringent bound, which can be refined through a Bayesian analysis. In the following we closely follow Shao and Zhang (2017). Given our theoretical hypothesis (${\rm DM{}}_{\rm obs}={\rm DM{}}_{\rm astro}+{\rm DM{}}_{\rm milli}$), and the set of measurements of ${\rm DM{}}_{\rm obs}$ from $N$ pulsars, we construct the log-likelihood as

[TABLE]

Here $\sigma_{i}$ is the dispersion for each pulsar, obtained adding in quadrature statistical uncertainties on ${\rm DM{}}_{\rm obs}^{i}$ and the astrophysical ones on ${\rm DM}_{\rm astro}^{i}$. We used a uniform prior on $\epsilon/m_{\rm milli}>0$ and verified that our results do not depend on this choice.

We shall consider two datasets of pulsars extracted from the ATNF Pulsar Catalogue Manchester et al. (2005) as explained in the Supplement Material. In both cases we assume a Navarro-Frenk-White profile for the DMa density, normalized to a local value of $\rho_{\rm dm}\approx 0.3\,{\rm GeV}/{\rm cm}^{3}$. The first dataset comprises $N=13$ local pulsars with the smallest values of ${\rm DM{}_{\rm obs}}/d$ and for which parallax measurements of the distance $d$ are available. We only choose pulsars located away from the galactic plane. This is to minimize the effect of the evacuation of DMa from the galactic plane for millicharged DMa. While early studies argue that this effect is relevant for \epsilon\gtrsim 5.4\times 10^{-22}\Big{(}\frac{m_{\rm milli}}{\rm eV}\Big{)} McDermott et al. (2011); Chuzhoy and Kolb (2009), a recent study Dunsky et al. (2018) suggests that this bound may be too restrictive. We also consider a second dataset of $N_{\rm cluster}=13$ pulsars located in globular clusters within $8\,{\rm kpc}$ from the galactic center and off the disk, again with the smallest ${\rm DM{}_{\rm obs}}/d$. Distances of clusters can be determined by different methods Krauss and Chaboyer (2003) not relying on the DM, and their uncertainty is usually of a few percent. We therefore assign a conservative error of $10\%$ to the value of $d$ for the pulsars in this second dataset. Even if the effect of the galactic magnetic field on the density of millicharged DMa away from the galactic disk is uncertain, we do not expect DMa to be evacuated at high galactic latitudes, and our analysis should provide realistic constraints.

For each pulsar we compute ${\rm DM{}}_{\rm astro}^{i}\approx\langle n_{e}\rangle_{i}d_{i}$, where $\langle n_{e}\rangle_{i}$ is an average electron density along the line of sight obtained using the YMW16 model M. Yao et al. (2016), while $d_{i}$ is the pulsar distance obtained from parallax (for the first dataset) or from the location of the globular cluster (for the second dataset). In the former case, we assign $\langle n_{e}\rangle_{i}$ a $20\%$ error to take into account potential systematics in the electron density model. This is a conservative approach given the uncertainties in M. Yao et al. (2016). We perform a Monte-Carlo Markov chain analysis using the Python ensemble sampler Emcee Foreman-Mackey et al. (2013) to explore the posterior distribution. For our datasets, $10^{5}$ samples are accumulated with $20$ chains. The chains show good acceptance rate and convergence. The results are similar for the two datasets:

[TABLE]

which we compare to other existing bounds in Fig. 1. In particular these results are compatible with $\epsilon=0$. For completeness, we also show a similar (weaker) bound estimated from the dispersion of the fast radio burst FRB 121102 Chatterjee et al. (2017). This line falls in the ballpark of the estimate (11). A more comprehensive analysis for FRBs will be presented elsewhere Blas et al. (tion).

The mass range in Fig. 1 is limited on the left because the expression (4) is valid as long as the energy of the photon is smaller than $m_{\rm milli}$. For radio waves from pulsars, $m_{\rm milli}\gtrsim\omega\sim{\rm GHz}\sim 10^{-6}~{}{\rm eV}$. Since the bound is more stringent for small masses, these constraints could improve as $1/m_{\rm milli}$ for sub-GHz pulsar measurements in systems with properties similar to the ones used in our analysis. Low-frequency measurements are indeed possible, see e.g. Ref. Pilia et al. (2016), though we leave a more systematic study of the sources for the future. Figure 1 shows that our bounds are competitive for masses below the Tremaine-Gunn bound on fermionic DMa, $m_{TG}\gtrsim\rm keV$ Tremaine and Gunn (1979). Hence, they apply to scalar charged DMa or to models with a fraction of millicharged fermionic DMa (see Eq. (13) for the scaling of the bound with $\rho_{\rm milli}$).

Finally, the existence of milli-charge DMa  also impacts the cosmological 21-cm line and distortions of the CMB Ali-Ha moud et al. (2015); Munoz and Loeb (2018); Slatyer and Wu (2018). It seems possible that these observations also constrain the very light case considered here, though previous studies focus on much heavier DMa candidates, and it seems cautious not to extrapolate their conclusions at much lower masses. Instead, it would be interesting to extend these analyses to smaller masses in the future.

Polarization constraints on ALPs. We now consider the case where the millicharged particles are absent, $j^{\nu}_{\rm milli}=0$. As discussed before, the modification of the TOA from the terms depending on $g$ in Eq. (7) is negligible and we ignore it. Nevertheless, due to their pseudo-scalar nature, ALPs also induce an oscillating variation of light polarization Harari and Sikivie (1992); Ivanov et al. (2018); Sigl and Trivedi (2018); Fujita et al. (2018); Plascencia and Urbano (2018); Obata et al. (2018); Ejlli (2018). Parity-symmetry breaking leads to birefringence, i.e. different phase velocities for left- and right-handed modes, which in turn induces rotation of the linear polarization plane. At first approximation, we assume the ALP-DMa background in the Milky Way rest frame to be described by the field configuration Arvanitaki et al. (2018)

[TABLE]

where $\sigma_{0}\approx 10^{-3}$ corresponds to the virialized velocity of the Milky Way and $\varphi_{\bf v}$ are arbitrary phases. The value $\tilde{\phi}_{0}$ changes smoothly with $x$ to reproduce the DMa energy density. Finally, for this non-relativistic configuration one can assume that $\omega_{v}\approx m_{a}(1+v^{2}/2)$. For low DMa masses, this field configuration has only long modes as compared to the wavelength of radio signals and an eikonal approximation can be used to study the propagations of waves in this continous background Weinberg (1962). The leading result of this calculation yields an effect for the polarization angle of a photon propagating from time $t$ to $t+T$ Ivanov et al. (2018); Harari and Sikivie (1992)

[TABLE]

where $\delta$ is a phase over which we will marginalize. The characteristic time scale for the axion background oscillation is $T_{\rm ALP}\sim\frac{10^{-22}{\rm eV}}{m_{a}}\,{\rm yr}$; if one continuously observes the polarized light from the source during a time $t_{\rm obs}\gtrsim T_{\rm ALP}$, the observed variation of the polarization angle (15) may constrain the amplitude of the axion oscillations111Notice that after exploring a quarter of a period of oscillation, the original value of $\delta$ is not relevant. Hence, even if a system lives in a region with $\delta\ll 1$, the previous analysis is valid for masses satisfying $t_{\rm obs}\gtrsim T_{\rm ALP}$. , i.e. the coupling $g$ for a given mass $m_{a}$. Pulsars are observed for long periods and the polarization angle is measured to be almost constant with a precision of roughly one degree, that can be compared with Eq. (15). We use the polarization data from Ref. Yan et al. (2011) and in particular PSR J0437$-$4715, which is the pulsar with the highest number of observations of the polarization angle, spanning a period of roughly four years. The ionospheric contribution to the polarization angle was subtracted using the program Getrm-Iono Han et al. (2006). Similar results are obtained when the ionospheric contribution is subtracted with the alternative FARROT method developed at the Dominion Radio Astrophysical Observatory (DRAO), Penticton, Canada. We performed a likelihood estimation of the coupling $g$ for a set of fixed masses $m_{a}$. For each value of the mass, we marginalize over the unknown phase $\delta$ in Eq. (15) in the interval $[-\pi,\pi]$ and then obtained the $95\%$ C.L. exclusion value for $g$, which is our reported constraint. There is a caveat in using the bound from a single system: it may be that the pulsar of interest lives in a region where the amplitude of the field (14) is lower than expected from the NFW profile. This situation may happen, for instance, in certain ULDM models where the field $\phi$ interpolates between different domains of condensation. The chances for this to happen are slim. Still, it is important to take this caveat into consideration. The use of more pulsars in the future will likely reduce this possibility even more.

The excluded region in Fig. 2 spans roughly four orders of magnitude in the mass range, from $m_{a}\sim 10^{-19}\,{\rm eV}$ to $m_{a}\sim 10^{-23}\,{\rm eV}$. The lower limit is set by the total observation time ($\sim 4\,{\rm yr}$), whereas the upper limit is set by the resolution time in the data-set during each observation run (‘folding time’), that is roughly 1 hour for J0437$-$4715. The derived lower bounds scale as $1/m_{a}$ — with some modulation due to the fact that observations of the polarization angle for J0437$-$4715 are not homogeneous in time — and are stronger for smaller masses, i.e longer observation time. The bound scales as $\sim\sqrt{\rho_{\rm dm}}$, so it can be competitive even if ALPs form only a small fraction of the DMa. This is particularly important at low masses, where other astrophysical constraints require the mass of the ALP to be  $m_{a}\gtrsim 10^{-21}\,$eV if it constitutes all the DMa. These bounds are based on the clustering properties of the DMa candidate at small scales  Kobayashi et al. (2017), the modifications of rotation curves in the inner regions of galaxies Bar et al. (2018), and the mere existence of galaxies with very small gravitational binding energies Marsh and Niemeyer (2018). These constraints are subject to independent astrophysical uncertainties, though together they indicate that masses below $m_{a}\sim 10^{-21}\,$eV are in tension with current data. In Fig. 2 we represent the previous limit by a conservative line at $m_{a}=4\times 10^{-21}\,$eV. This limitation relaxes for fractional components.

Discussion. Several DMa models introduce dispersion effects in the photon propagation. Although small, these effects accumulate for photons coming from astrophysical sources and can be constrained through precision measurements. The effect of millicharged DMa is degenerate with that of ordinary plasma and improving models for the local plasma distribution will help strengthening the constraints from DM. On the other hand, the effect of ALP-photon coupling is more striking and requires a careful analysis of the TOA as a function of the frequency. In addition, in the upcoming era of the Square Kilometre Array, we will benefit from a much larger pulsar sample (possibly comprising sources near the galactic center, where the DMa density is higher than what assumed here), combined with a significantly improved timing precision Kramer and Stappers (2015); Shao et al. (2015); Bull et al. (2018). The prospects of using radio waves in probing DMa are very promising in the near future. For ALPs, their coupling to photons generates an oscillation of the polarization angle of photons in the ultra-light DMa case. Our results in Fig. 2 show that, for the mass range $10^{-23}-10^{-20}$eV, the constraints derived here are the best available and will greatly improve in the future with more data.

We have considered propagation in a weak magnetic field for which dispersion due to the ALP-photon coupling and QED vacuum polarization effects are negligible. However, our formalism can be easily extended to include such effects, which might be relevant for propagation in strongly magnetized regions. A discussion of this effect will appear elsewhere Blas et al. (tion).

Note: While this work was close to completion, Ref. Liu et al. (2019) appeared on the arXiv, estimating constraints on ALPs using the polarization angle of radio waves from pulsars similar to those derived in the second part of our work. Even though the idea is similar, our analysis, based on real data, is distinct and the results differ from the ones in Liu et al. (2019) by roughly a factor $\sqrt{\frac{10^{-22}{\rm eV}}{m_{a}}}\sqrt{\frac{400{\rm pc}}{d}}$ originating from a different assumption about the $\phi$ configuration.

Acknowledgments. We are grateful to Nikita Blinov, Richard Brito, Anson Hook, Georg Raffelt, Günter Sigl for interesting discussions. We thank Davide Racco and Mikhail M. Ivanov for pointing out a mistake in the first arXiv version of this work. AC acknowledge support from national grants FPA2014-57816-P, FPA2017-85985-P and the European projects H2020-MSCAITN-2015//674896-ELUSIVES and H2020-MSCA-RISE2015. PP acknowledge financial support provided under the European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480, and support from the Amaldi Research Center funded by the MIUR program “Dipartimento di Eccellenza” (CUP: B81I18001170001) and by the GWverse COST Action CA16104, “Black holes, gravitational waves and fundamental physics.” Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. LS was partially supported by the National Science Foundation of China (11721303), and XDB23010200.

Appendix A Supplemental material

We provide here additional details of the datasets analyzed in this work. For the millicharged DMa, we analyzed a first set of galactic pulsars selected for their minimal ${\rm DM}/d$, where $d$ is derived from parallax, and for their good agreement with the electron density model (Table 1). A second set of pulsars is selected in galactic clusters (Table 2). In this case, in addition to the aforementioned criteria, we also require that the pulsars are not further from the galactic center than the Solar System, $\sqrt{X^{2}+Y^{2}+Z^{2}}<8.3\,{\rm kpc}$, and are also located far from the galactic disk, $|Z|>1\,{\rm kpc}$.