Stochastic Wave Dark Matter with Fermi-LAT γ-Ray Pulsar Timing Array

Pulsar timing arrays (PTAs) can detect disturbances in the fabric of spacetime on a galactic scale by monitoring the arrival time of pulses from millisecond pulsars (MSPs). Recent advancements have enabled the use of γ-ray radiation emitted by MSPs, in addition to radio waves, for PTA experiments. Wave dark matter (DM), a prominent class of DM candidates, can be detected with PTAs due to its periodic perturbations of the spacetime metric. In response to this development, we perform in this Letter a first analysis of applying the γ-ray PTA to detect the ultralight axion-like wave DM, with the data of Fermi Large Area Telescope (Fermi-LAT). Despite its much smaller collecting area, the Fermi-LAT γ-ray PTA demonstrates a promising sensitivity potential. We show that the upper limits not far from those of the dedicated radio-PTA projects can be achieved. Moreover, we initiate a cross-correlation analysis using the data of two Fermi-LAT pulsars. The cross-correlation of phases, while carrying key information on the source of the spacetime perturbations, has been ignored in the existing data analyses for the wave DM detection with PTAs. Our analysis indicates that taking this information into account can improve the sensitivity to wave DM by ≳50% at masses below 10−23 eV.


INTRODUCTION
Wave dark matter (DM) represent a class of bosonic DM candidates with a mass ≲ 1 eV (for a review on wave DM, see, e.g., Hui 2021).Such a candidate has a large occupation number at a scale of de Broglie wavelength in the environment of Milky Way galaxy and hence demonstrates a coherent wave nature.Wave DM with a mass ≲ 10 −18 eV is usually dubbed "ultralight" (Marsh 2016), hnluu@connect.ust.hktaoliu@ust.hkrenjing@ihep.ac.cn for its astronomical-scale de Broglie wavelength.The "Fuzzy" DM (FDM, Hu et al. 2000;Schive et al. 2014;Hui et al. 2017), as one special class of ultralight wave DM with a mass ∼ 10 −22 eV, is of particular interest in astronomy.As its de Broglie wavelength is comparable to the size (∼ 0.5 kpc) of dwarf-galaxy dark cores, the FDM can be applied to address several small-scale problems on astronomical structures, such as core-cusp problem (de Blok 2010) (namely a contradiction between the predicted Navarro-Frenk-White "cusps" for the DM halo profile (Navarro et al. 1996) in N-body simulations of standard collisionless cold DM, and the observation of approximately constant DM density in the central parts of galaxies) which are challenging the standard cold DM scenario.
Ultralight wave DM, including FDM, could be detected by employing its impacts on either the high-z Universe (Marsh 2016), such as the formation of largescale structure, or the astronomy of local galaxy, e.g., the timing and polarization of pulsar pulses.The successful launching of James Webb Space Telescope provides a great tool for achieving the goal with the first method, while the active and upcoming pulsar timing programs enable us to explore this puzzle with the second one.
Given that pulsars can emit electromagnetic pulses with extraordinary regularity, a Pulsar Timing Array (PTA, Sazhin 1978;Detweiler 1979;Hellings & Downs 1983) consisting of broadly distributed and welltimed millisecond pulsars has been proposed as an astronomical interferometer to detect the gravitational waves (GWs) of 10 −9 − 10 −6 Hz.In 2013, Khmelnitsky and Rubakov pointed out that the wave DM can periodically perturb galactic metric, at a level of O(v 2 ), and hence be detected with the PTA (Khmelnitsky & Rubakov 2014).The concept of Pulsar Polarization Array (PPA) is much newer (Liu et al. 2023).This new methodology was developed recently to reveal new phenomena in astronomy with a common correlated polarization signal, such as the one caused by ultralight axion-like wave DM (Liu et al. 2023).Because its Chern-Simons coupling with electromagnetic field strength is parity-violating, the axion-like wave DM can modulate the position angle of linearly polarized pulsar light while the light travels across the DM halo.The PPA is thus highly suited for its detection (Liu et al. 2023).Notably, the PTA and PPA methods here are based on gravitational and nongravitational effects, respectively.So combining them can further strengthen the capability of arrayed pulsars for exploring the nature of DM.
Due to the good sensitivity of large radio telescopes, the PTAs have been mostly realized by monitoring radio millisecond pulsars.In addition to the ongoing programs such as the Parkes PTA (PPTA, Manchester et al. 2013), the European PTA (EPTA, Desvignes et al. 2016), the Chinese PTA (CPTA, Lee 2016) and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, NANOGrav Collaboration et al. 2015) PTA, in the near future the international Square Kilometer Array radio telescope (Dewdney et al. 2009) is also expected to underpin the construction of global PTA.These dedicated efforts were extended to a higher-frequency band recently.With the data accumulated by Fermi Large Area Telescope (Fermi-LAT), the first γ-ray PTA (γ-PTA) constraints for the stochastic GWs have been delivered (FERMI-LAT Collaboration et al. 2022).While its sensitivity is limited by statistics, such a high-frequency PTA benefits from a suppression of intrinsic red noise of dispersion measure variance and relative weak effects from white noises such as radio interference and jitter.The γ-PTA thus well-complements the radio one.Such a "radio + high-frequency" strategic development could be extended to the PPA measurements in the near future, with the launch of dedicated X-ray polarimetry programs such as the missions of Imaging X-ray Polarimetry Explorer (IXPE, Weisskopf et al. 2022) and Enhanced X-ray Timing and Polarimetry (eXTP, Zhang et al. 2019).For these telescopes, a precision ∼ O(1 • ) can be achieved for measuring polarization angle of bright sources such as bright pulsars (Doroshenko et al. 2023) and accretion flows around black holes (De Rosa et al. 2019).
Given its high scientific value, in this Letter we will demonstrate a first application of high-frequency PTA for detecting the ultralight axion-like wave DM, with the Fermi-LAT data.

METHODOLOGY FOR DETECTING STOCHASTIC WAVE DM
The profile of stochastic wave DM can be described as an uncorrelated superposition of particle waves.For a proof of concept, let us consider ultralight axion-like wave DM.Its profile can be put as Here C v is equal to α v f (v)(∆v) 3/2 with v ∈ Ω denoting lattice sites in phase space and ∆v being their spacing.ρ(x) and f (v) represent energy density profile and velocity distribution.The stochastic nature of wave DM is then encoded as a random amplitude α v ∈ (0, +∞) and phase ψ v , which follow Rayleigh distribution p(α v ) = α v e −α 2 v /2 and uniform distribution respectively (for a proof, see Sec.II in Foster et al. 2018).
Because of its gradient energy, the wave-DM halo in a galaxy can induce an oscillating gravitational potential Ψ(x, t) in proportion to its local density.This further yields a shift to the time of arrival (TOA) of pulsar pulses (Khmelnitsky & Rubakov 2014): where L p is a distance from the p-th pulsar to the Earth, x i denotes the position of the Earth (i = e) and pulsar (i = p), and ρ i = ρ(x i ) is the DM energy density at x i .∆t p (t) is then determined as a difference between an "Earth" and a "pulsar" terms.These terms oscillate with a frequency ∼ 2m a , with the frequency being set by the mass of axion-like particles m a .For the oscillation phases, the spatial dependence on pulsars is encoded in φ (e) vv ′ and φ (p) vv ′ contain two random phases due to a quadratic dependence of ∆t p (t) on the wave DM profile a(x, t).This is different from the polarization signal of axion-like wave DM where such a dependence is linear (Liu et al. 2023).
The radio-PTA analysis is usually performed with a TOA-based method, which involves the comparison of a time-referenced model of the pulse profile with the observed data.For many of its pulsars, the Fermi-LAT can only reconstruct a few TOAs per year due to the low statistics of received γ-ray photons from these pulsars.Thus, its ability to detect short-time-scale (f ≳ 1 yr −1 ) signals gets highly constrained in this case.A photonby-photon (PBP)-based method, which we will take below, thus has been introduced in FERMI-LAT Collaboration et al. (2022) to address this limitation.
At the Fermi-LAT (Atwood et al. 2009), the γ-ray data are recorded in terms of energy E i and arrival time t i for the i-th photon.The relevant physical observable is pulsar spin phase (PSP) at certain photon arrival time, i.e., 0 ≤ ϕ(t i ) < 1.The evolution of PSP can be predicted with a timing model.For the p-th pulsar, it is defined as There are in total three classes of contributions to ϕ (p) .ϕ tm denotes the contributions from a deterministic timing model, encoding the effects of pulsar spindown, position, proper motion, etc.For example, the spindown effect yields a contribution ϕ ), with ν and ν being pulsar spin frequency and its spin down rate, respectively.δϕ tn denotes the contributions arising from potential timing-noise processes such as red spin noise.It can be written as in terms of Fourier modes (FERMI-LAT Collaboration et al. 2022): Here, T ≈ 12.5 yrs represents the total observation time, which is approximately equal for all pulsars.ν (p) 0 de-notes the intrinsic spin frequency, and k's are the indices of Fourier amplitudes.The red spin noise can be modeled as a Gaussian process, with a covariant matrix element for the Fourier mode where P tn (f ) = A 2 tn 12π 2 f /yr −1 −Γtn yr 3 denotes the power spectrum, with A tn and Γ tn being its amplitude and spectral index, respectively, and f k = k/T is the frequency of k-th Fourier mode.At last, δϕ s sin(2m a t).This results in where These coefficients follow a multivariate Gaussian distribution with zero mean.Assuming that the DM velocity is isotropic and peaks sharply at the virial velocity of Milky Way, i.e., f (v) ≈ (4πv 2 ) −1 δ(v − v 0 ) with v 0 ∼ 10 −3 , the covariance matrix elements are given by with L pq = L p − L q and y ij = |x i − x j |/l c .Here l c ≡ 1/(m a v 0 ) is the coherence length of wave DM. sin y ij /y ij is a sinc function, encoding the spatial correlation of signal.To derive Eq. ( 7), ⟨α 2 v ⟩ = 2 and ⟨sin(α have been applied.Compared to the stochastic GWs with a coherence length l c ∼ 1/ω, the wave DM of the same frequency has a much longer l c due to its non-relativistic character.The "pulsar" term thus plays a more important role in cross-correlating the signals (Liu et al. 2023).This significance can be further enhanced for the pulsars close to the galactic center where the DM halo is dense (De Martino et al. 2017;Liu et al. 2023).
For a PTA with N pulsars, the log-likelihood function in the PBP method is given by FERMI-LAT Collaboration et al. ( 2022) Here the terms of first line describe a weighted Poisson distribution of PSP with a summation over the photons from all pulsars, f (p) denotes a template known in prior fitting the data, and w (p) i measures the probability for a i-th photon to originate from the p-th pulsar.The rest is simply the Gaussian log-likelihood for random noise and signal processes.The terms of C tn arise from red spin noise, with N1 , ..., β The terms of C a account for the contributions of stochastic wave DM, with 0 A (1) s , ..., ν Here the elements of 2N × 2N covariant matrix (C a ) pq are calculated in Eq. ( 8).Cross-correlating the data of pulsars thus can help to distinguish the signals from the uncorrelated noises and identify the nature of signal, as the Hellings-Downs curve does for the GW detection (Hellings & Downs 1983) 1 .

SENSITIVITY ANALYSIS
We analyze the Fermi-LAT γ-ray data for 28 out of 35 pulsars (Kerr & Parthasarathy 2022), excluding J0312-0921, J0613-0200, J1513-2550, J1543-5149, J1741+1351, J1858-2216, J2034+36322 , with the log-likelihood in Eq. ( 10).Since the distance of these pulsars to the Earth is small compared to the scale of DM density variations, we assume ρ p ≃ ρ e in Eq. ( 8) and then absorb ρ e into a dimensionless parameter Ψ c = πGρ e /m 2 a .Here Ψ c is an amplitude measure for the Ψ(x, t) oscillation induced by the axion-like wave DM.It can be viewed as a counterpart of h c /f 1/2 , namely the GW characteristic strain with a frequency factor.To derive the upper limit on Ψ c , we first analytically marginalize over the nuisance parameters ϕ 0 , ν, ν, β tn and β a to derive the marginal likelihood L m , following the procedure in FERMI-LAT Collaboration et al. (2022).Then, we numerically marginalize over the amplitude A tn and spectral index Γ tn for the red spin noise, and the pulsar distance {L p } for p = 1, ..., N in L m3 .The exclusion limits of Ψ c or ρ e at 95% C.L. are finally computed with where x = {A tn , Γ tn , L p } and P (x), P (Ψ c ) denote the relevant priors.For the dimensionless parameter Ψ c , we set a uniform prior with Ψ c ∈ U 10 −20 , 10 −10 .For the red spin noise, we set A tn ∈ U 10 −20 , 10 −10 and Γ tn ∈ U [1, 6].These prior ranges have been chosen such that the associated PSP for the signals does not exceed the timing model contribution and hence can be treated perturbatively with respect to the timing model.We employ the package PTMCMCSampler (Ellis & Van Haasteren 2017) to perform the numerical marginalization on nuisance parameters for the cases where the timing noises and/or the distance uncertainties are involved.
As a "fiducial" demonstration, let us consider the sensitivities with the terms independent of pulsar distance only in Eq. ( 8), for both auto-and cross-correlations.The covariance matrix (C a ) pq in Eq. ( 12) is reduced to In this case, L m is also independent of the pulsar distance, and the numerical marginalization of {L p } in Eq. ( 13) becomes trivial.This fiducial signal model will be referred to as simplified correlation below.
With this fiducial signal model, we display the upper limits of the Fermi-LAT γ-PTA on Ψ c and ρ e (95% C.L.) in Fig. 1 and 2 respectively, as a function of the wave DM mass m a .As the Fermi-LAT is not a dedicated PTA mission, its statistics for the PTA is relatively limited.The sensitivity potential demonstrated for the γ- PTA however is quite promising4 .For m a ≲ 1/T ∼ 10 −23 eV, these limits are comparable to the PPTA and NANOGrav ones (Afzal et al. 2023).For larger m a values, they become approximately one order of magnitude weaker than the PPTA limits (Porayko et al. 2018).Interestingly, as the TOA shift ∆t p is ∝ Ψ c for the axionlike wave DM and ∝ h c /f 1/2 for the GWs, these γ-PTA limits exhibit a scaling behavior with m a (or frequency) similar to those obtained with a TOA-based method on the characteristic strain (Hazboun et al. 2019) after factoring out f 1/2 .Approximately, we have Ψ 95 ∝ m a for m a > 1/T ∼ 10 −23 eV, and Ψ 95 ∝ m −2 a for m a < 1/T .For both PTA and PPA detections, cross-correlating the data from the arrayed pulsars plays a crucial role in recognizing the nature of signals, which might be manifested as a specific pattern for the signal spatial correlation, or a distinguishable way for the correlated signals to depend on the relative position of these pulsars.For the axion-like wave DM, such an effect arises from a full consideration of its wave phase, or more explicitly the ∼ m a v • x term in this phase (see Eq. ( 1)), in calculating the covariance matrix C a in Eq. ( 12), which however has been ignored in previous studies (Porayko et al. 2018;Afzal et al. 2023;Smarra et al. 2023) and our analysis so far.As indicated in Eq. ( 8), such a treatment trivializes the Earth and pulsar crossing terms also in the correlation functions.These terms depend on the pulsar spatial position and exist for both auto-correlation and cross-correlation.
To look into the pulsar cross correlation, next let us extend the sensitivity analysis to the case where all terms in Eq. ( 8) are turned on.We consider J1231-1411 and J0614-3329, namely the two pulsars yielding the strongest upper limits individually, instead of using the full data of 28 pulsars as above.The pulsar spatial coordinates and their distance to the Earth are taken from the ATNF Pulsar Catalogue (Manchester et al. 2005).For the marginalization in Eq. ( 13), we take the priors for pulsar distance introduced in Appendix A. We also assume no timing noise for simplifying the analysis, as its impacts on the exclusion limits tend to be tiny (see Fig. 1 and Fig. 2).Three signal models are then comparatively analyzed: simplified correlation where (C a ) pq is simplified as Eq. ( 14), auto-correlation only where (C a ) p̸ =q terms are turned off in Eq. ( 12) and full correlation where all terms in Eq. ( 12) are turned on.We demonstrate the upper limits on Ψ c in Fig. 3, as a function of m a .The cross correlation plays a particularly important role for m a ≲ 0.2 kpc/v 0 ∼ 10 −23 eV, i.e., while the de Broglie wavelength of wave DM becomes bigger than the spatial separation between the two considered pulsars.Because of its contributions, the upper limits of full correlation get visibly strengthened compared to those of auto-correlation.We find this feature to be generic by performing the same analysis with the Fermi-LAT data of various pulsar pairs.For larger m a values, all sinc function terms in Eq. ( 8) tend to be suppressed.Then, the contributions of auto-and cross-correlations to the covariance matrix and hence to the sensitivities become comparable.Notably, while the measurement of pulsar distance to the Earth is challenging, the information of cross-or full-correlations remains largely unaffected unless the distance uncertainties be-come larger than the de Broglie wavelength of wave DM 5 (see Appendix A for details).

DISCUSSION AND CONCLUSION
Before concluding, let us take a review on the results demonstrated in last section.Let us start with a comparison of constraints between different radio PTA experiments.As shown in Fig. 1 and Fig. 2, the EPTA constraints are stronger than the NANOGrav and PPTA ones in the low mass regime (m a ≲ 1/T ∼ 10 −23 eV), and become comparable to them for the high mass regime (m a ≳ 1/T ).The gap in the low mass regime could be partly attributed to the longer data span of EPTA, which is 24.7 years, in comparison to NANOGrav's 15 years and PPTA's 12 years.As discussed in Smarra et al. (2023), such a difference may result in an enhancement of constraints on Ψ c and ρ e by several times.At large masses, the longer data span of EPTA becomes less helpful in improving the constraints.The cadence instead becomes more relevant to the detection and recognization of fast-oscillating signals from the dominant pulsar white noises (Smarra et al. 2023).
Unlike these radio PTAs, the Fermi-LAT is not a dedicated PTA mission.But, the time span of 12.5 years for its pulsar observation (FERMI-LAT Collaboration et al. 2022) is only slightly shorter than that of NANOGrav.Particularly, the Fermi-LAT γ-PTA has collected the data uninterruptedly during its entire operation due to the constant experimental setup.Such a stability in data taking (and also an overall gain from the suppression of the intrinsic red noise of dispersion measure variance) makes the Fermi-LAT γ-PTA especially promising at masses below 10 −23 eV (FERMI-LAT Collaboration et al. 2022).The constraints comparable to the NANOGrav ones have thus been achieved in this mass regime, despite the telescope's much smaller collecting area.At high masses, however, the limitation stemming from the small number of photons becomes evident, resulting in a sensitivity gap between the Fermi-LAT γ-PTA and the radio PTAs which is enlarged gradually as m a increases.
Following the approach introduced in Hazboun et al. (2019), we can simply estimate the sensitivity prospects 5 Note, as the GWs are relativistic, their de Broglie wavelength are about three orders of magnitude shorter than that of wave DM, for a given frequency.So, different from what occurs to the wave DM, the "pulsar terms" contributions to the PTA detection of stochastic GWs are subject to a strong suppression caused by the distance uncertainty of pulsars, for the frequency range of interest.For example, at nano-Hertz the "pulsar terms" can contribute coherently for the stochastic GWs only if the distance uncertainty is ≲ 1 pc, in comparison to ≲ 1 kpc for the wave DM.
of detecting ultralight wave DM with the γ-PTA.In the low mass regime, the sensitivities are roughly scaled with ρ e,95 ∝ T −3 .This implies that, by extending the observation time by another 12.5 years, the Fermi-LAT constraints could get close to the EPTA ones (Smarra et al. 2023) and begin ruling out a major contribution from wave DM.As for the high mass regime, the Fermi-LAT sensitivity will not improve significantly due to the dominance of white noise.To further improve the sensitivity, more dedicated γ-ray telescopes with significantly larger collecting areas and instant luminosities, either groundbased or spaceborne (Aharonian et al. 2001;Fan et al. 2022;Xie et al. 2023), are needed.
Finally, let us comment on the prospects of generalizing the full correlation γ-PTA analysis from two to more pulsars.Unlike the simplified correlation model, the calculation of likelihood function for the full correlation model involves inverting the covariance matrix C a with off-diagonal components.As the pulsar number increases, the complexity of this process grows quickly, as ∼ 8N 3 .Furthermore, the dimension of parameter space becomes higher with more pulsars.More sampling points are thus required to fully explore the posterior distributions.Additionally, a more accurate analysis with red noise turned on would request more computing resources.It is therefore important to develop new likelihood calculation techniques for the future full correlation data analysis with more pulsars.
In summary, we have explored the potential of detecting the stochastic axion-like wave DM with the γ-PTA, using the recently released Fermi-LAT data.As it benefits from a suppression of the intrinsic red noise of dispersion measure variance, such a high-frequency PTA wellcomplements the traditional radio ones.Inspiringly, despite the much smaller collecting area of telescope, the upper limits generated for the energy density of axionlike wave DM are not far from those obtained in the more dedicated radio-PTA programs such as PPTA and NANOGrav, particularly for a mass range below the inverse of observation period.Such a strategic development could be extended to the PPA measurements in the near future.By then we would expect the hypothesis of wave DM to be examined gravitationally and nongravitationally, in both radio and high-frequency bands.This will tremendously advance our exploration on the nature of DM.
[Note added] While this Letter was being finalized, Xia et al. (2023) appeared on arXiv.The application of γ-PTA for detecting ultralight wave DM was explored in both papers, with the Fermi-LAT data.However, this Letter distinguishes itself from Xia et al. (2023) in terms of both signal nature (stochastic versus deterministic) and analysis method (PBP-based versus TOA-based).Particularly, a methodology to cross-correlate the data of arrayed pulsars, which is crucial for identifying the signal nature, has been developed and the relevant data analysis was fulfilled.
denotes the PSP deviation induced by a signal such as the stochastic wave DM.It has the same format as that of δϕ a (t) is related to the signal ∆t p (t) in Eq. (2) as δϕ

Figure 1 .
Figure1.Upper limits of the Fermi-LAT γ-PTA (95% C.L.) on the dimensionless measure Ψc, as a function of ma."SS" and "DS" in the legend are abbreviations of "stochastic signal" and "deterministic signal".The limits are computed in the signal model with simplified correlation.The red dashed and solid curves denote the limits obtained with the Fermi-LAT data of 28 pulsars, with and without the timing noises, respectively.The best single-pulsar limits, which arise from pulsar J1231-1411, are shown as the green curves.As a reference, we present the most recent radio-PTA limits obtained with the PPTA data (black-dotted,Porayko et al. 2018),  the EPTA data (black-dashed, Smarra et al. 2023), and the NANOGrav data (black-solid,Afzal et al. 2023).Additionally, we show the theoretical prediction for Ψc as a light-blue band by setting ρe to a value of local DM density (de Salas & Widmark 2021).The lower sub-panel shows the same upper limits but with an extended mass range to highlight the scaling of Ψ95 with ma.

Figure 2 .
Figure 2. Upper limits of the Fermi-LAT γ-PTA (95% C.L.) on the local DM density ρe, as a function of ma.The description of curves is the same as that in Fig. 1.

Figure 3 .
Figure 3. Upper limits of the Fermi-LAT γ-PTA (95% C.L.) on Ψc, as a function of ma.The limits are computed in three signal models: simplified correlation, auto-correlation only and full-correlation, with the Fermi-LAT data from the two pulsars yielding the strongest upper limits individually.The lower sub-panel shows the relative change of these limits with respect to the full-correlation analysis: δΨc = Ψc/Ψ c,full -1.