High-precision search for dark photon dark matter with the Parkes Pulsar Timing Array

The nature of dark matter remains obscure in spite of decades of experimental efforts. The mass of dark matter candidates can span a wide range, and its coupling with the Standard Model sector remains uncertain. All these unknowns make the etection of dark matter extremely challenging. Ultralight dark matter, with $m \sim10^{-22}$ eV, is proposed to reconcile the disagreements between observations and predictions from simulations of small-scale structures in the cold dark matter paradigm, while remaining consistent with other observations. Because of its large de Broglie wavelength and large local occupation number within galaxies, ultralight dark matter behaves like a coherently oscillating background field with an oscillating frequency dependent on its mass. If the dark matter particle is a spin-1 dark photon, such as the $U(1)_B$ or $U(1)_{B-L}$ gauge boson, it can induce an external oscillating force and lead to displacements of test masses. Such an effect would be observable in the form of periodic variations in the arrival times of radio pulses from highly stable millisecond pulsars. In this study, we search for evidence of ultralight dark photon dark matter (DPDM) using 14-year high-precision observations of 26 pulsars collected with the Parkes Pulsar Timing Array. While no statistically significant signal is found, we place constraints on coupling constants for the $U(1)_B$ and $U(1)_{B-L}$ DPDM. Compared with other experiments, the limits on the dimensionless coupling constant $\epsilon$ achieved in our study are improved by up to two orders of magnitude when the dark photon mass is smaller than $3\times10^{-22}$~eV ($10^{-22}$~eV) for the $U(1)_{B}$ ($U(1)_{B-L}$) scenario.

traditional particle like cold dark matter models show that the central density profile of dark matter in dwarf galaxies is much steeper than that inferred from observed rotation curves, and the number of satellite galaxies of Milky Way-like hosts predicted from simulations is higher than that inferred from observations [2]. These so-called "core-cusp problem" and "missing-satellite problem" impose challenges to the the conventional cold dark matter paradigm. Baryonic feedback effects [3] may be a viable solution, but it is a nontrivial task. Meanwhile, various alternative dark matter models have been proposed to address such small-scale shortcomings, e.g., warm dark matter [4], self-interacting dark matter [5], and fuzzy dark matter [6].
A dark photon is a hypothetical particle from the theory beyond the Standard Model of particle physics. Just like the ordinary photon, a dark photon is the gauge boson of a U(1) interaction. The existence of the dark photon is well predicted in many string-inspired models, such as large-volume string compactifications [7][8][9]. The dark photon mass can be generated by either the Higgs mechanism or the Stueckelberg mechanism, and it is naturally light [10]. When this U(1) symmetry is the baryon or lepton number or their linear combination, protons, neutrons, and electrons are "charged" under this symmetry, and there will be an extra force (i.e., the "fifth force") between objects that consist of ordinary matter.
In recent years, the dark photon has been widely considered as a viable dark matter candidate. Interestingly, when the dark photon is ultralight (m A ∼ 10 −22 eV), which represents a realization of the fuzzy dark matter paradigm, its de Broglie wavelength is up to subgalactic scale, i.e., O(0.1 − 1) kpc, and the occupation number in dark matter halos is very large [11]. These aspects imply that the "wave nature" of dark matter particles is pronounced; hence the dark matter can be properly treated as an oscillating background field rather than individual particles. Because of the formation of a "soliton core" instead of a density cusp at the center of galaxies [12][13][14], the density of the dark matter halo at a galactic center becomes flat. This provides a better fit to observations than a cuspy Navarro-Frenk-White profile [15] predicted in the cold dark matter scenario. In addition, the substructures formed under the fuzzy dark matter hypothesis are fragile against tidal disruptions, since they usually have less concentrated mass distributions. This solves the "missing-satellite" problem [16] at the same time. To date, the ultralight dark photon dark matter (DPDM) leads to predictions consistent with most existing observations at various scales [16], making it one of the most compelling dark matter candidates.
Conventional direct detection experiments [17][18][19] are not sensitive to fuzzy dark matter particles because of the extremely small energy and momentum depositions in elastic scatterings. If the ultralight dark matter particle is a U(1) B or U(1) B−L gauge boson (i.e., dark photon), the fifth-force experiments can be used to put stringent constraints on the existence of such particle [20], as well as the energy loss of compact binary systems [21]. Note that, a more general U(1) interaction that includes both U(1) B and U(1) B−L is discussed in earlier work [22,23], while the fifth-force constraint is derived in Ref. [24]. Furthermore, the dark matter background can cause displacements on terrestrial/celestial objects, leading to observable effects [25]. For example, depending on the mass of dark photon particles which determines the dark matter oscillation frequency, such motion can be probed using high-precision astrometry [26], spectroscopic [27], timing observations [28,29], as well as gravitational-wave detectors [30][31][32].
Among all these relevant observations, the pulsar timing array (PTA) experiments are of particular interest to us. Millisecond pulsars with very high rotational stability are observed to emit periodic electromagnetic pulses with incredible accuracy. For the best pulsars, the measurement uncertainties of pulse arriving times are below 100 ns. The information of possible new physics (including signatures of DPDM) is in the irregularities of pulse arrival times, which are usually called timing residuals. The U(1) B or U(1) B−L DPDM with an ultralight mass ∼ 10 −22 eV, induces displacements of Earth and pulsar that cause a periodic signal in timing residuals with frequency [33] Hz. This frequency falls into the sensitive region of PTA experiments, which have collected measurements of pulse arrival times on weekly to monthly cadence over time scales of 10-20 years.
In this work, we make use of the second public data release of the Parkes PTA project [35,36], to search for evidence of the DPDM. Subsets of these data have been used to search for stochastic gravitational waves [37] as well as the gravitational effects induced by ultralight scalar dark matter [28]. Compared with Ref. [28], we discuss the direct coupling between DPDM and the ordinary matter in this work, which results in different signals from the gravitational effect studied in Ref. [28]. We also improve the analysis through considering the correlation of pulsars in the DPDM field, the interference of DPDM wave functions, and excluding fake signals that arise from the time-frequency method. (see Secs. I and II of the Supplemental Material [38] for details.) Observations. -The second data release (DR2 [41]) of the PPTA project [36] includes high-precision timing observations for 26 pulsars collected with the 64-m Parkes radio telescope in Australia. The data span is 14.2 years, from 2004 February 6 to 2018 April 25, except for PSR J0437−4715 where pre-PPTA observations from 2003 April 12 are also included. The end date of this data set corresponds to the installation and commissioning of a new ultrawide bandwidth receiver on the Parkes telescope. The observations were made typically once every 3 weeks in three radio bands (10, 20, and 40/50 cm). Details of the observing systems and data processing procedures are described in Refs. [35,36].
Pulsar times of arrival (ToAs) in PPTA DR2 are obtained using the Jet Propulsion Laboratory Solar-System ephemeris DE436 and the TT(BIPM18) reference time scale published by the International Bureau of Weights and Measures (BIPM). We fit these ToAs to a timing model with the standard TEMPO2 [39,40] software package. We perform the Bayesian noise analysis and search for DPDM signals using the enterprise [42] package. The PTMCMC [43] sampler is used for the stochastic sampling of parameter space and the calculation of Bayesian upper limits. We also use PyMultiNest [44] to calculate the Bayes Factor while performing model selection.
The names of pulsars, basic observing information, and noise properties of the PPTA DR2 pulsars are listed in Supplemental Table S1 [38].
Timing model. -The ToAs of a pulsar include quite a few terms. There are several astronomical effects that should be accounted for, such as the intrinsic pulsar spin-down, the proper motion of the pulsar, the time delay due to interstellar medium, and the Shapiro delay caused by Solar-System planets. In addition to these deterministic effects, both uncorrelated noise (white noise) and correlated noise (red noise) may be present in pulsar timing data [45][46][47][48][49][50]. We model the ToAs as follows: (1) where t TM is the "timing model" that accounts for the deterministic effects, ∆t Noise represents stochastic noise terms, and ϕ TM and ϑ Noise are parameters in the timing model and noise model, respectively. ∆t DPDM is used to describe the "signal" term caused by the DPDM with paramaters ϑ DPDM . We make use the TEMPO2 and enterprise software packages to determine the timing model and marginalize over model uncertainties in our Bayesian analysis. More details can be found in Refs. [28,51].
The noise term includes white noise and red noises from the rotational irregularities of the pulsar, the variations of the dispersion measure, and the band noise. The noise is assumed to follow a multivariate Gaussian distribution with a covariance matrix. For the red noises, power-law frequency dependence is assumed with parameters independently fitted for each noise term. We describe the modelings of noise in detail in the Supplemental Material [38].
The DPDM-induced ToA residuals. -Within coherence region where l < l c ∼ 2π/(m A v vir ), the DPDM can be approximately described as a plane wave, A(t, x) = A 0 sin m A t − k · x + α(x) . Here m A is the dark photon mass, the phase term α is a function of x [12,13], and k is the characteristic momentum. The direction of k is random and its magnitude is ∼ m A v vir where v vir is the virial velocity in our Galaxy. A 0 is the gauge potential of the DPDM background, whose direction is another random vector and is unrelated to that of k in the non-relativistic limit. The averaged value of the magnitude | A 0 | 2 is determined by the local dark matter energy density, 2ρ 0 /m 2 A , with ρ 0 = 0.4 GeV cm −3 near the Earth [52]. The interference among dark photons induces a random O(1) fluctuation to the magnitude. This plane-wave approximation breaks down when the measurements are performed at a time scale longer than the coherence time, t c ∼ 2π/(m A v 2 vir ), or at a length scale larger than the coherence length l c . Note that we only consider the vector components of the gauge potential for dark photon field because we adopt the Lorentz gauge where the contribution from the scalar component is negligible [30]. In the Supplemental Material [38] we describe the simulations of the DPDM distribution in the local vicinity of the Solar System.
The weak coupling between the dark photon background and the "dark charge" results in an acceleration of a test body [30] as where m is the mass of the test body, an characterizes the coupling strength of the new gauge interaction that is normalized to the electromagnetic coupling constant e. The dark charge q equals the baryon number for the U(1) B interaction or baryon-minus-lepton number for the U(1) B−L interaction. Such an acceleration allows the detection of the DPDM in several ways, e.g., using high-precision astrometry [26] or gravitational-wave detectors [30]. The acceleration causes a displacement to a test object, which is approximately Note that here we neglect the contribution to the displacement from the spatial part −k · x + α(x). The virial velocity of dark matter is v vir ∼ 10 −3 c. For the dark photon mass range of interest, 10 −23.5 eV≤ m A ≤ 10 −21 eV, the coherence length ranges from 0.04 to 13 kpc. For an observation of O(10) years, the proper motion of a celestial body (the Sun or the pulsar) is about 3 × 10 −3 pc which is much smaller than the coherence length of the DPDM field, assuming a proper motion velocity of ∼ 10 −3 c. The phase change for the Earth and pulsars, induced by the −k · x + α(x) term, can be safely ignored and the DPDM field can be treated as spatially uniform around the Earth or the pulsar. Therefore we can rewrite Eq. (3) as here A 0 e,p and α e,p represent the amplitude and phase of the dark photon field at the locations of the Earth and pulsars, respectively.
Most of the pulsars studied in this work have distances between 0.1 and several kpc, which are comparable to the coherence length. Therefore, we perform the analysis in two limits: (A) Completely uncorrelated: The phases and amplitudes of dark photon background for each pulsar are independent.
(B) Completely correlated: The phases are independent phases but with a common amplitude.
In both cases, we treat the phases of the DPDM field at each pulsar as independent free parameters because the locations to most of these pulsars are highly uncertain. Such an uncertainty results in a random phase to dark photon field characteristic to each pulsar. Note that for m A > 10 −22 eV we only perform the uncorrelated analysis (A), since in this high-mass range most of the pulsars should lie in the uncorrelated regime (see Supplemental Fig. S2 [38]). A hybrid of correlated and uncorrelated treatment can in principle be applied based on the comparison between the coherence length and the separation of each pulsar pair. We expect the results of such a hybrid analysis to fall between the two cases discussed here.
We further assume a flat spacetime and obtain the time residual, ∆t DPDM , induced by the DPDM where t and t are the times when the pulsar emits and the Earth receives a pulse respectively, and d is the position vector pointing from the Earth to the pulsar, and n = d/|d|. The timing residual caused by the DPDM is anisotropic. Specifically, the Earth term is a dipole contribution in terms of angular correlation, which is distinct from the monopole signal that derives from the gravitational effect of the general fuzzy dark matter model [28,53]. Now we write the timing residuals explicitly as follows where q (B) e,p , q (B−L) e,p , m e and m p are the B charge, B − L charge, the mass of the Earth and the pulsar, respectively. For U(1) B , q (B) /m is approximately 1/GeV for both the Earth and the pulsar. For U(1) B−L , q (B−L) /m is about 1/GeV for the pulsar and 0.5/GeV for the Earth. Note that we absorb the time difference between a pulsar and the Earth into the phase term.
Statistical analysis -We perform a likelihood-ratio test, which assesses the goodness of fit using the following statistic where H 0 is the null hypothesis that the timing residuals contain only the noise contributions, and H 1 is the signal hypothesis that the DPDM signal is present. In both hypotheses the noise parameters ϑ W (white noise parameters), ϑ SN (spin noise parameters), ϑ DM (dispersion-measure noise parameters), and ϑ BN (band noise parameters) are fixed at their best-fit values obtained from single pulsar analyses. The free parameters include ϑ BE (Bayes ephemeris parameters, vary in both H 0 and H 1 ; see Ref. [54] for details) and ϑ DPDM (DPDM parameters, vary only in H 1 ). The priors of all the parameters mentioned above are listed in Supplemental Table S2 [38]. We adopt the Bayesian approach to derive the constraints on the DPDM parameters. We set the upper limit on the coupling constant,¯ , at 95% credibility as where ϑ DPDM are the DPDM parameters excluding m A and , and P are the prior probabilities of those parameters (see Supplemental Table S2 [38]). Here the dark photon mass m A is singled out from ϑ DPDM as we search for possible signals for a given range of masses. In the computation of the upper limits, we fix both the white noise and the red noise parameters as their best-fit values obtained in single-pulsar analyses.
Results and discussion -Our main results about the 95% upper limits on the coupling constant 2 , where is the dimensionless coupling constant and e is the coupling strength of the U(1) B or U(1) B−L interaction with e being the electromagnetic coupling, are shown in Fig. 1. For comparison, the constraints on these parameters from experiments testing the weak equivalence principle (WEP), e.g., Ref. [20], are also presented. See the Supplemental Material [38] for a description of the WEP results. We find that our study imposes significantly improved limits in the low-mass region with m A 3 × 10 −22 eV for U(1) B and m A 10 −22 eV for U(1) B−L . Compared with the WEP results, the PPTA constraints on are stronger by one to two orders of magnitude.
We find that for a few particular frequencies, the signal hypothesis fits the PPTA DR2 data better than the null hypothesis. The strongest "signal" is at f 1.02 × 10 −7 Hz. The corresponding m A and are given in Supplemental Table S3 [38]. We note that the favored value of the coupling constant has already been ruled out by existing WEP experiments [20]. Therefore, it may be due to some unmodeled noise effect. A similar spurious signal also appeared in the analysis Ref. [28], which was speculated to be caused by the unmodeled perturbations of Solar System barycenter from Mercury.
In order to have a better understanding of the "signal" at f 1.02 × 10 −7 Hz, we reduce the number of pulsars and perform a consistency test. We choose two sets of pulsars: The first one includes five pulsars that contribute the most to this strongest "signal", the second one is the same except that we do not included PSR J0437-4715. Additionally, we allow red-noise parameters ϑ SN , ϑ DM , and ϑ BN to vary along with ϑ DPDM and ϑ BE . We calculate the Bayes factor between the null and the signal hypothesis. We find that the "signal" is mainly contributed by PSR J0437-4715. For the analysis with the other four pulsars excluding PSR J0437-4715 gives a logarithmic Bayes factor of 37 in favor of the H 1 hypothesis at f 6.58 × 10 −8 Hz, which is significantly smaller than the five-pulsar analysis. The results of the tests are presented in Supplemental Table S4 [38]. Note also that the best-fit frequency differs from that of the five-pulsar analysis, and the best-fit coupling constant lies above the WEP upper limit, indicating that this may be due to unmodelled noise. More careful studies need to be carried out in the future to understand the sources of these systematics.
We note that the mass range with the strongest constraint is below the typical fuzzy dark matter mass (10 −22 eV). Currently there is no consensus on the constraints on the mass of the fuzzy dark matter. The stellar kinematics of dwarf spheroidal galaxies tend to favor a lower mass [56,57] (a few times of 10 −23 eV), while the halo mass function derived from the Lyman-α forest suggests a higher lower limit [16,27] ( 10 21 eV). Most of these constraints rely on as- sumptions about the structure formation in the fuzzy dark matter scenario. More detailed studies of the interplay between the fuzzy dark matter and the baryonic matter, such as the dark matter-to-baryon mass ratio or various baryonic effects in fuzzy dark matter simulations, are required to reach a robust and consistent result. The gray shaded regions in Fig. 1 indicate the parameter space for which the "gravitational effect" due to the spacetime metric oscillation induced by the wavelike DPDM field dominates over the gauge interaction discussed in this work. See the Supplemental Material [38] for more details of the estimate of the gray regions. In such parameter regions, dedicated analysis with both effects being included in the signal model is required and will be studied in future work.
Summary -Ultralight fuzzy dark matter is proposed as an attractive candidate of dark matter in the universe which helps solve the small-scale crises of the classical cold dark matter scenario. Using the precise timing observations of 26 pulsars by the PPTA project, we study the possible couplings between ultralight dark matter and ordinary matter. Taking the DPDM as an example, we obtain by far the strongest constraints on the parameters of the DPDM. The upper limits on the dimensionless coupling constant derived in our study are improved by up to two orders of magnitude when the dark photon mass is smaller than 3×10 −22 eV (10 −22 eV) for the U(1) B (U(1) B−L ) scenario.
The search sensitivity for the DPDM is expected to improve significantly in the near future as more pulsars are monitored with continually extending data spans by the worldwide PTA campaigns including the PPTA, the North American Nanohertz Observatory for Gravitational Waves (NANOGrav [58]), and the European PTA (EPTA [59]), which have jointly formed the International PTA (IPTA [60]). The Five-hundredmeter Aperture Spherical Telescope (FAST [61]) and the Square Kilometer Array (SKA [62]) are also expected to join the IPTA collaboration. These efforts are likely to bring revolutionary progress in studying a wide range of dark matter models, and more generally in answering the related fundamental questions in physics.
Acknowledgements Following [28], we assume ∆t Noise follows a multivariate Gaussian distribution with a covariance matrix C tot . In the most general case, for each pulsar, C tot contains four terms: where N accounts for white noise, C SN , C DM and C BN arise from red-noise processes. C SN is the spin noise that describes the time-correlated spin noise which might be induced by rotational irregularities of the pulsar. C DM accounts for the time variations of the dispersion measures. C BN is the band noise which is used to account for unknown additional red noise that is only present in a specific radio-frequency band [48]. In our analysis, we only include band noise for two pulsars that are known to exhibit significant band noise [48], PSR J0437-4715 and J1939+2124. We apply the band noise for these two pulsars in four bands: 10 cm, 20 cm, 40 cm and 50 cm. All the terms in Eq. (10) are n × n matrices, with n being the number of ToAs for each pulsar. In practice, the white noise of the actual timing data is usually larger than the ToA measurement uncertainties. Empirically, for each receiver/backend system and each pulsar, a re-scaling factor EFAC (Error FACtor) is multiplied to the measurement uncertainty, and an extra white noise component EQUAD (Error added in QUADrature) is also included [S1, S2, S3]. The rescaled ToA uncertainty σ s is related to the original uncertainty σ We model the power spectra of the spin noise, band noise and dispersion-measure noise as power-law functions, where A and γ are the noise amplitude and noise spectral index at a reference frequency of yr −1 , respectively. The noise covariance matrices are, here t i and t j are i-th and j-th ToAs, ν obs i and ν obs j are their respective radio frequencies, f L is the low-frequency cutoff, ν ref = 1400 MHz is the reference radio frequency. Earlier work shows that setting f L = 1/T , with T being the observational time interval, is sufficient for pulsar timing analysis due to the fitting for the quadratic spin-down in the timing model [45,51]. Note that in Eq. (13), the elements of band-noise covariance matrics is non-zero only when radio frequencies of i-th and j-th ToAs belong to the same frequency band.
We adopt the "time-frequency" method to calculate the integral in Eq. (13,14). Namely, we approximate the integral as a discrete summation of different Fourier frequency modes. In this approximation, a uniform interval of Fourier frequency modes is usually adopted up to an upper limit. One may choose an interval of d f ≡ ∆ f = 1/T , which we call "T uniform". However, such a choice is sub-optimal [46]. Instead, we adopt the "2T uniform" sampling where ∆ f = 1/(2T ) for the approximation of Eq. (13,14). We found a better approximation is achieved with the "2T uniform" interval. On this basis, we find that additional periodic patterns may appear due to the approximation, when γ is small. We also compare the influence of the total number of Fourier frequency modes (timefrequency modes) M when M = 30 and M = 40, and find that the frequencies of their periodic patterns are different. This is used to remove the potential fake signal introduced during this calculation. The discussions above are demonstrated in Fig. S1.
In the "time-frequency" method, the covariance matrices in Eq. (13,14) can be rewritten as Here i is the index for ToAs, and α is used to label the discretized Fourier frequency. We follow Ref. [28] in terms of the definition and calculation of the likelihood function.
The parameter description and prior ranges of the noise model is given in Table S2. In Table S1 we also present the posterior credible intervals for spin noise and dispersionmeasure noise based on the results of single pulsar analyses, where we analyze the noise properties of each pulsar independently. In Bayesian analyses detailed below, we fix the noise parameters at their maximal likelihood values obtained from single pulsar analyses in some cases to reduce the computational costs. They are referred to as the "best fit" parameter values, which are available in our code repository.

B. Simulation of the DPDM background
Here we provide details on how to simulate the DPDM background. The interference among the wave-functions of these dark photon particles induces fluctuations on the dark matter field strength. Here we study this amplitude distribution, and it will be used to inform our priors in the statistical analysis. In addition, by studying the spatial correlations among various quantities of the dark photon background, we demonstrate the concept of coherence length. Within the coherence length, the DPDM can be approximated as a plane wave, whereas this approximation breaks down once the spatial separation becomes larger. The dark photon background A (tot) can be modeled as the summation of many, N 1, freely propagating planewaves, where n is the index for each dark photon particle and α (n) is the phase term in a plane wave. A (n) 0 is the gauge potential vector, whose magnitude is the same for all particles, but the direction is isotropically distributed. ω (n) and k (n) are the particle's energy and momentum, distributed according to Boltzmann velocity distribution with virial velocity v vir ∼ 220 km/s. The normalization of the dark photon field is determined by the local dark matter energy density ρ 0 , where V and T need to be sufficiently large, much larger than the coherence volume and time. Since v vir c, we have ω (n) m A . For the mass range that we are interested in, the coherence time is much larger than the observation time. Also as discussed previously, the phase variations, induced by the motions of pulsars and the Earth in the Galaxy during the period of observation, are negligible. Thus in our analysis, we can approximate the wave-function for the DPDM as where the index i labels the spatial component of the dark photon field, i.e. {x, y, z}. The averaged magnitude of the gauge potential component can written as For simplicity we introduce a normalized DPDM amplitudẽ To study the statistic behaviors of the DPDM, we set N = 2 15 and simulate the background field for 10 5 times. This ensures the results we obtain are statistically stable. We first focus on the DPDM amplitude. In the left panel of Fig. S2, we show the distribution of (Ã (tot) 0,x ) 2 at x = 0. This distribution can be well fit by an exponential function, which will be used as the prior probability of (Ã (tot) 0,i ) 2 in our Bayesian analysis. Next we study the coherence of the dark photon field as a function of the spatial separation between two points. The coherence between two quantities, X and Y, can be described by their two-point correlation function as where µ X , µ Y are expectation values, and σ X , σ Y are standard deviations. Let us consider X ∈ {Ã (tot) 0,i (0, 0, 0), α i (0, 0, 0)} and Y ∈ {Ã (tot) 0,i (d, 0, 0), α i (d, 0, 0)}, where d represents the separation between two points. In the right panel of Fig. S2, we show corr(X, Y) as a function of spatial separation. Here we clearly see that the correlation function reduces significantly once the spatial separation is longer than the coherence length, i.e., the de Broglie wavelength. We compare our result with the separations among the Earth and pulsars. We find that when m A > 10 −22 eV, it is safe to assume that the dark photon field at these locations are completely uncorrelated.

C. Limits from the WEP experiment
The experiments testing the weak equivalence principle (WEP) can also impose stringent constraints on the coupling constant of the U(1) B or U(1) B−L gauge force. In the WEP experiments, the dark photon is not required to be the dark matter. Here we compare our results with those obtained by WEP experiments. For the mass regime that we are interested in this study, the best limit is given by the satellite experiment, MICROSCOPE [20]. The effect is characterized by the Eötvös parameter as η E−s e 2 2 4πG where q 1,2,E and m 1,2,E are the "dark charges" and masses of the two test bodies in the MICROSCOPE satellite and the Earth, respectively. The best sensitivity achieved by the MI-CROSCOPE experiment [20] is |η| < 2.7 × 10 −14 at 2σ level. This can be translated to constraints on coupling constants as > 10 −23.28 for U(1) B and > 10 −24.07 for U(1) B−L . The results are consistent with Ref. [24]. In Fig. 1, we find our results improves upon that of MICROSCOPE when m A 3 × 10 −22 eV for U(1) B and m A 10 −22 eV for U(1) B−L , by up to two orders of magnitudes.

D. Gravitational effect of DPDM on the ToA
The quantum pressure of the ultralight bosonic dark matter field can induce a periodic oscillation of the metric. This will affect the propagation of photons and lead to detectable oscillatory effects at pulsar timing experiments. Such an effect has been studied in Refs. [28,53] in the case of ultralight scalar field. Ultralight vector field, such as the dark photon, can also be searched in this way. Thus it is instructive to compare this pure gravitational effect with the effect studied in this paper.
Here Ψ c ≡ Gρ 0 /(π f 2 ) with f = m DM /π, and m DM is the mass of the dark matter particle. The root-mean-square value of ∆t is where we take the maximal value of vector dark matter's timing residual. Meanwhile, the effect induced by the gauge coupling of the DPDM has the root-mean-square value as (∆t (B) dp ) 2 0.577 e 2ρ 0 where the prefactors are obtained through our numerical simulation. By comparing the root-mean-square value of these two effects, we find that the gauge interaction of the dark photon induces a stronger effect when × m A /eV > 10 −48.35 for U(1) B and 10 −48.24 for U(1) B−L . This comparison is shown as the black dashed line in Fig. 1.