Introduction

Although there is little doubt that dark matter (DM) exists in our Universe, its nature remains mysterious, including its component mass(es). It may be an ultralight elementary particle, such as fuzzy DM with mass \(\sim\! 1{0}^{-22}\) eV1,2,3,4, or it may arise from stellar-mass objects, such as primordial black holes5.

One promising DM candidate in the ultralight mass regime is the dark photon (DP), which is the gauge boson of a \(U(1)\) gauge group. The DP can acquire its mass through the Higgs or Stueckelberg mechanism. As ultralight DM, the DP must be produced non-thermally, e.g., production from the misalignment mechanism6,7,8, parametric resonance production or tachyonic instability of a scalar field9,10,11,12, or from the decay of a cosmic string network13.

It was recently proposed in refs. 14,15 that a gravitational wave (GW) detector may be sensitive to dark photon dark matter (DPDM). The Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) consists of two 4-km dual-recycled Michelson Fabry–Perot interferometers in Livingston Louisiana (L1) and Hanford, Washington (H1). From the first two observing runs (coincident with the Virgo detector for several weeks of the O2 run), detections of ten binary black hole mergers and one binary neutron star merger have been reported16. These measurements require a differential strain measurement sensitivity better than 10\({}^{-21}\) for broadband transients with central frequencies of \(O\)(100 Hz), based on detecting minute changes in distance between the mirror pairs forming the Fabry–Perot interferometer arms.

Relevant to this search, the mirror separations can also change in response to a gradient in a DPDM field due to non-zero photon velocity. More explicitly, we consider a DP with mass \({m}_{{\rm{A}}}\) between \(1{0}^{-13} \sim 1{0}^{-11}\) eV. The DPDM is an oscillating background field, for which the rest-frame oscillation frequency satisfies: \({f}_{0}\approx \left[\frac{{m}_{{\rm{A}}}}{1{0}^{-12}\ {\rm{ev}}}\right]\)(241 Hz). We assume the DP is the gauge boson of the gauged \(U{(1)}_{{\rm{B}}}\) group so that any object, including a LIGO mirror, that carries baryon number will feel its oscillatory force, similar to that experienced by a macroscopic, electrically charged object in an oscillating electric field.

Using LIGO to look for DPDM bridges GW science and particle physics. In this paper, we present a \(U{(1)}_{{\rm{B}}}\) DPDM search using data from Advanced LIGO’s first observing run, O1. We confirm that the data from LIGO’s first observation run yields results already more sensitive than limits from prior experiments in a narrow DPDM mass range. The sensitivity will be improved significantly with more LIGO data, as well as with the growth of the global network of GW detectors. Meanwhile, the same search strategy can be directly applied to search for many other ultralight DM scenarios, with excellent sensitivities achievable.

Results

Estimating DPDM-induced effects

Through virialization, DPDM particles in our galaxy have a typical velocity around \({v}_{0}\equiv 1{0}^{-3}\) of the speed of light, and thus they are highly non-relativistic. The total energy of a DM particle is then the sum of its mass energy and kinetic energy, i.e., \({m}_{{\rm{A}}}(1+{v}_{0}^{2}/2)\). Here and in the following, we use natural units, i.e.,  \(c=\hslash =1\). Therefore, the oscillation frequency of the DP field is approximately a constant, \(\omega \simeq {m}_{{\rm{A}}}\), with \(O(1{0}^{-6})\) deviations.

Therefore, within a small period of time and spatial separation, the DP field can be treated approximately as a plane wave, i.e.,

$${A}_{\mu }\simeq {A}_{\mu ,0}\cos [{m}_{{\rm{A}}}t-{\bf{k}}\cdot {\bf{x}}+\theta ].$$
(1)

Here \({A}_{\mu ,0}\) is the amplitude of the DP field and \(\theta\) is a random phase. The DP field strength can be simply written as \({F}_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }\). We choose the Lorenz gauge, \({\partial }^{\mu }{A}_{\mu }=0\), in what follows. In the non-relativistic limit, the dark electric field is much stronger than the dark magnetic field, and \({A}_{t}\) is negligible relative to \({\bf{A}}\). The magnitude of the DP field can be determined by the DM energy density, i.e., \(| {{\bf{A}}}_{0}| \simeq \sqrt{2{\rho }_{{\rm{DM}}}}/{m}_{{\rm{A}}}\).

In Eq. (1), we neglect the kinetic energy contribution to the oscillation frequency. We also set the polarization and propagation vectors, i.e., \({{\bf{A}}}_{0}\) and \({\bf{k}}\), to be constant vectors. This approximation is valid only when the observation is taken within the coherence region, i.e., \({t}_{{\rm{obs}}}\ <\ {t}_{{\rm{coh}}}\simeq \frac{4\pi }{{m}_{{\rm{A}}}{v}_{{\rm{vir}}}^{2}}\) and \({l}_{{\rm{obs}}}\ <\ {l}_{{\rm{coh}}}\simeq \frac{2\pi }{{m}_{{\rm{A}}}{v}_{{\rm{vir}}}}\). For example, if the DP field oscillates at 100 Hz, the coherence time is only \(1{0}^{4}\) s, much shorter than the total observation time. In order to model the DPDM field for a time much longer than the coherence time, we simulate the DPDM field by linearly adding up many plane waves propagating in randomly sampled directions. More details are given in the “Methods” section below.

From the DPDM background field \({\bf{A}}(t,{{\bf{x}}}_{i})\), one can derive the acceleration induced by the DPDM on each test object, labeled by index \(i\). This can be written as:

$${{\bf{a}}}_{i}(t,{{\bf{x}}}_{i})\simeq \epsilon e\frac{{q}_{{\rm{D}},i}}{{M}_{i}}{\partial }_{t}{\bf{A}}(t,{{\bf{x}}}_{i}).$$
(2)

Here we use the approximation \({\bf{E}}\simeq {\partial }_{t}{\bf{A}}(t,{{\bf{x}}}_{i})\) for the dark electric field. \({q}_{{\rm{D}},i}/{M}_{i}\) is the charge–mass ratio of the test object in LIGO. Treating a DP as the gauge boson of \(U{(1)}_{{\rm{B}}}\), and given that the LIGO mirrors (test masses) are primarily silica, \({q}_{{\rm{D}},i}/{M}_{i}=1/{\rm{GeV}}\). We note that results from ref. 17 impose very strong constraints on the gauged \(U{(1)}_{{\rm{B}}}\) scenario, due to gauge anomaly. However, the results derived in these papers rely on an assumption of how to embed the model into a complete theory at high energy in order to cancel \(U{(1)}_{{\rm{B}}}\) anomalies. Extending the electroweak symmetry breaking sector or other anomaly cancellation mechanism can avoid such severe constraints. If one takes the model independent constraint on an anomalous gauge symmetry, new particles need to be added at energy scale \(O(\frac{4\pi {m}_{{\rm{A}}}}{\epsilon e})\), which gives \(O({\rm{TeV}})\) for the parameter space we are interested in. Since the required new particles carry only electroweak charges, they are safe from various collider searches. We label the DP–baryon coupling as \(\epsilon e\) where \(e\) is the \(U{(1)}_{{\rm{EM}}}\) coupling constant. We emphasize that we choose DP to be a \(U{(1)}_{{\rm{B}}}\) gauge boson as a benchmark model. The same analysis strategy presented in this study can be directly applied to many other scenarios, such as a \(U{(1)}_{{\rm{B}}-{\rm{L}}}\) gauge boson or scalar field, which couples through Yukawa interactions. More details on various models, as well as subtleties when observation time is longer than coherence time, will be described in the future work.

Signal-to-noise ratio (SNR) estimation

We approximate the DPDM field as a plane wave within a coherence region. For a DP field oscillating at frequency \(O(100)\) Hz, the coherence length is \(O(1{0}^{9}\ {\rm{m}})\), much larger than the separation between the two LIGO GW detectors at Hanford and Livingston. Thus these two GW detectors experience a nearly identical DPDM field, inducing strongly correlated responses. Exploiting the correlation dramatically reduces the background in the analysis.

The DPDM signal is exceedingly narrowband, making Fourier analysis natural. We first compute discrete Fourier transforms (DFT) from the time-domain data. The total observation time is broken into smaller, contiguous segments, each of duration \({T}_{{\rm{DFT}}}\), with a total observing time \({T}_{{\rm{obs}}}={N}_{{\rm{DFT}}}{T}_{{\rm{DFT}}}\). Denote the value of the complex DFT coefficient for two interferometers 1 and 2, DFT \(i\), and frequency bin \(j\) to be \({z}_{1(2),ij}\). The one-sided power spectral densities (PSDs) for two interferometers are related to the raw powers as \({{\rm{PSD}}}_{1(2),ij}=2{P}_{1(2),ij}/{T}_{{\rm{DFT}}}\). \({P}_{1(2),ij}\) are taken to be the expectation values for \(| {z}_{1(2),ij}{| }^{2}\), as estimated from neighboring, non-signal frequency bins, assuming locally flat noise (using a 50-bin running median estimate).

To an excellent approximation, the noise in the two LIGO interferometers is statistically independent, with the exception of particular very narrow bands with electronic line disturbances18, which are excluded from the analysis. Detailed descriptions of broadband LIGO noise contributions may be found in ref. 19, including discussion of potential environmental contaminations that could be correlated between the two LIGO detectors, but none of which would mimic a DPDM signal. The normalized signal strength using cross-correlation of all simultaneous DFTs in the observation time can be written as

$${S}_{j}=\frac{1}{{N}_{{\rm{DFT}}}}\sum _{i=1}^{{N}_{{\rm{DFT}}}}\frac{{z}_{1,ij}{z}_{2,ij}^{* }}{{P}_{1,ij}{P}_{2,ij}}.$$
(3)

In the absence of a signal, the expectation value is zero and the variance of the real and imaginary parts is

$${\sigma }_{j}^{2}=\frac{1}{{N}_{{\rm{DFT}}}}{\left\langle \frac{1}{2{P}_{1,j}{P}_{2,j}}\right\rangle }_{{N}_{{\rm{DFT}}}},$$
(4)

where \({\langle \rangle }_{{N}_{{\rm{DFT}}}}\) denotes an average over the \({N}_{{\rm{DFT}}}\) DFTs, which may have slowly varying non-stationarity. SNR can be defined by

$${{\rm{SNR}}}_{j}\equiv \frac{{S}_{j}}{{\sigma }_{j}}.$$
(5)

Taking into the account the small separation between the interferometers relative to the DP coherence length and their relative orientation (approximate 90-deg rotation of one interferometer’s arms projected onto the other interferometer’s plane), we expect the SNR\({}_{j}\) for a strong DPDM field to be primarily real and negative.

Efficiency factor

In order to use the observed real(SNR) values to set limits on DPDM coupling as a function of frequency, one must correct for the signal power lost from binning. The suggested nominal binning proposed in ref. 15 is \(\Delta f/f=1{0}^{-6}\), based on a Maxwell velocity distribution. The bin size in frequency space is set by \({T}_{{\rm{DFT}}}\), i.e., \(\Delta f=1/{T}_{{\rm{DFT}}}\), which is optimal at only \({f}_{{\rm{opt}}}\simeq 1{0}^{6}/{T}_{{\rm{DFT}}}\). For a frequency higher than \({f}_{{\rm{opt}}}\), the relative frequency binning is finer, implying loss of signal power in single-bin measurements. At frequencies lower than \({f}_{{\rm{opt}}}\), the relative frequency binning is coarser, implying full capture of the signal power but at the cost of unnecessarily increased noise. We note that it is possible the DM velocity distribution deviates from Maxwell distribution by an \(O(1)\) factor, e.g., refs. 20,21. However, the impact is small, as the single-bin search used here depends on the integrated power within a frequency bin and not so much on its shape.

In Fig. 1, we show the DPDM signal power spectrum as a function of frequency, where \({f}_{0}={m}_{{\rm{A}}}/2\pi\). We choose to normalize the x-axis by the intrinsic signal width, determined by the typical kinetic energy of DM particles. In this calculation, we include the Earth rotation effect. Without it, the signal PSD is proportional to \(vf(v)\) where \(f(v)\) is the Maxwell distribution. The Earth’s rotation broadens our signal by \(\Delta f\approx 2{f}_{{\rm{E}}}\). Different choices of \({f}_{0}\) result in slightly different deformations after including the rotation, but the changes are negligible in the frequency regime of interest. An analytical understanding of the PSD will be presented in the future work.

Fig. 1: An example of dark photon dark matter signal power spectrum and corresponding detector sensitivity.
figure 1

The dark photon dark matter (DPDM) signal power spectrum is shown in terms of characteristic strains \({h}_{{\rm{c}}}\) (red), with \(U{(1)}_{{\rm{B}}}\) coupling parameter \({\epsilon }^{2}=1{0}^{-41}\), DPDM oscillation frequency \({f}_{0}=500\) Hz, and typical velocity of DPDM \({v}_{0}=1{0}^{-3}\) of the speed of light. The Advanced LIGO design sensitivity in a small frequency window is approximately flat, which is shown as the black dashed line.

The power spectrum from numerical simulation is used to determine empirically the fractions of power falling into a single fixed \(\Delta f/f\) bin, where bin boundaries are systematically varied over the allowed range. Figure 2 shows the resulting efficiencies (power fractions) for \({T}_{{\rm{DFT}}}\) set to be 1800 s. The red dotted curve shows the best case, for which the bin boundary is optimal. The blue dashed curve shows the worst case, which necessarily approaches \(50 \%\) for coarse binning (low frequency), while the green solid curve shows the average maximum efficiency over all bin boundary choices. A fit to the green solid curve is used for deriving upper limits on DPDM coupling.

Fig. 2: Signal power single-bin detection efficiency as a function of relative frequency resolution for a fixed coherence time of 1800 s.
figure 2

The upper (red) curve is for an optimal bin boundary choice (a priori unknown) for a given signal. The lower (blue) curve shows the worst-case efficiency for the least optimal boundary choice. The middle (green) curve shows an average over randomly chosen boundary choices.

Data selection and analysis

The strain data used in this analysis were downloaded from the Gravitational Wave Open Science Center (GWOSC) web site22 and transformed to create 1786 1800-s coincident DFTs from the L1 and H1 interferometers. The GWOSC data sets exclude short periods during which overall data quality is poor. The choice of coherence time in this first DPDM search is somewhat arbitrary but allowed convenient comparison of spectral line artifacts observed with those reported from 1800-s DFTs in LIGO continuous GW searches, for which 1800 s is a common DFT duration chosen. A shorter coherence time would be more optimal at frequencies above \(\sim\! 500\) Hz for this single-bin detection analysis. In principle, a longer time would be more optimal for lower frequencies, but in practice, sporadic interruptions of interferometer operations during data taking lead to significant livetime loss for DFTs requiring very long contiguous periods of coincident Hanford–Livingston operations.

The search for detections and the setting of upper limits in the absence of detection is based on “loud” values of the detection statistic (Eq. (5)). Specifically, we look for large negative real values of the SNR. Since we search over \(\sim\! 4\) million DFT bins in the band 10–2000 Hz, we must correct for a large statistical trial factor in assessing what SNR value is deemed “significant.” We choose a nominal signal candidate selection of SNR \(< -\! 5.8\), corresponding to a \(\sim \!1\)% false alarm probability, assuming Gaussian noise. In practice, the noise in some frequency bands is not truly Gaussian, leading to excess counts at large SNR. To assess the severity of this effect, we also define and examine control bands (“frequency lags”) in which a DFT frequency bin in one interferometer is compared to a set of offset bins from the other interferometer such that a true DPDM signal would not contribute to a non-zero cross-correlation but for which single-interferometer artifacts or broadband correlated artifacts lead to non-zero correlation. This frequency lag method is analogous to the time lag method used in transient GW analysis. Specifically, we choose 10 lags of (\(-50\), \(-40\), ..., \(-10\), \(+10\), ..., \(+50\)) frequency bin offsets to assess the non-Gaussian background from these instrumental artifacts. To avoid contamination of both signal and control bands from known artifacts, we exclude from the analysis any band within \(\sim\! 0.056\) Hz of a narrow disturbance listed in ref. 18, where the extra veto margin is to reduce susceptibility to spectral leakage from strong lines. We also exclude the band 331.3–331.9 Hz, for which extremely loud narrow calibration excitations in the two interferometers lead to significant overlapping spectral leakage and hence non-random correlation.

Figure 3 shows the distributions of the real and imaginary parts of the SNR (Eq. (5)) for both the signal bins (“zero lag”) in magenta and the lagged bins in black. The distributions follow quite closely the ideal Gaussian curve shown, except for a slight excess visible in the tails beyond \(| {\rm{SNR}}| \ > \ 5\) (note there are \(\sim\! 10\) times as many lagged bins as signal bins in the graphs). The only signal bins with \(| {\rm{SNR}}| \ > \ 5.8\) arise from known continuous wave “hardware injections” used in detector response validation, for which the complex SNR can have an arbitrary phase in the cross-correlation that depends on the simulated source frequency and direction. An investigation was carried out of all other SNR outliers (10) with real or imaginary values having magnitudes >5. In all but three cases, lagged bins in neighboring bins within 0.2 Hz of the signal bin showed elevated noise, defined by an SNR magnitude >4, suggesting non-Gaussian contamination. The Gaussian noise expectation for this range [5.0–5.8] of subthreshold outlier magnitude (real or imaginary) is 4.1 events, consistent with observation in clean bands.

Fig. 3: Distributions of the real and imaginary parts of the signal-to-noise ratio.
figure 3

The signal-to-noise ratio (SNR) for the signal bins (“zero lag”) are labeled in magenta and the lagged (control) bins in black, along with the ideal Gaussian expectation in green.

Since no significant candidates were found, upper limits were set. In future searches, should significant candidates appear, it will be critical to assess their consistency with instrumental artifacts. A simple approach is to increase the number of control bins examined per candidate to assess better potential non-Gaussian single-interferometer contamination and broadband correlated artifacts. A greater concern would be a highly narrowband correlated disturbance, such as from identical electronic instruments at each observatory creating a sharp spectral line through electrical current draws in power supplies affecting interferometer controls. Detailed investigation using auxiliary instrumental and environmental channels would be warranted, to exclude such interference.

New constraints from LIGO O1 data

Our main results are presented in Fig. 4. We show the derived \(95 \%\) confidence level upper limits on the parameter \({\epsilon }^{2}\) for DP–baryon coupling, as a function of DPDM oscillating frequency. The broad red band shows the range of upper limits obtained with \(1/1800\) Hz binning, using the measured real part of the SNR detection statistic defined below and the Feldman–Cousins (FC) formalism23 and after applying an efficiency correction discussed below. The yellow curve shows the expected upper limit for an average measured real(SNR) = 0, applying the same FC formalism and efficiency correction. The dark blue curve shows a more optimal upper limit expected when the DFT binning adjusts with frequency to maintain \(\Delta f/f=1{0}^{-6}\) for the same 893-h observation time, for the same efficiency correction, and for an averaged detector sensitivity equal to that in the analysis. The yellow and dark blue curves agree well with each other at around 500 Hz, where \(1/1800\) Hz is the optimal choice of the bin size. The mean achieved upper limit is generally worse than the optimal sensitivity, because with fixed bin size at \(1/1800\) Hz, excess noise is included at low frequency and some signal power is lost at high frequency. The dashed curve shows upper limits derived from the Eöt-Wash group based on Equivalence Principle tests using a torsion balance24,25. Given the LIGO O1 data, under the assumption that DP constitutes all DM, we have already improved upon existing limits in a mass window around \({m}_{{\rm{A}}} \sim 4\;\times\,1{0}^{-13}\) eV.

Fig. 4: Derived \(95 \%\) confidence level upper limits on the coupling parameter \({\epsilon }^{2}\) for dark photon-baryon coupling.
figure 4

The broad red band shows the actual upper limits with \(1/1800\) Hz binning. The yellow curve shows the expected upper limit for an average measured real (SNR) = 0. The dark blue curve shows the “optimal” upper limit expected when the discrete Fourier transform (DFT) binning adjusts with frequency to maintain \(\Delta f/f=1{0}^{-6}\) for the same 893-h observation time. The magenta curve shows the “optimal” upper limit expected for a 2-year, \(100 \%\)-livetime run at Advanced LIGO design sensitivity (“O4-O5"). The dashed curve shows upper limits derived from the Eöt-Wash group24,25. This is a fifth-force experiment, whose constraint does not rely on dark photon (DP) being dark matter (DM). The large spikes of red and blue curves, overlaid on top of each other, are induced by known sources of noise, such as vibrations of mirror suspension fibers.

Future searches in more sensitive data will probe deeper into an unexplored \({\epsilon }^{2}\)\({m}_{{\rm{A}}}\) parameter space. Assuming no discovery and a negligible true GW stochastic background, the magenta curve shows the “optimal” upper limit expected for a 2-year, \(100 \%\)-livetime run at Advanced LIGO design sensitivity (“O4-O5”). This limit looks smoother, as it uses a design sensitivity curve that shows only fundamental noise sources, while the blue curve includes additional, non-fundamental noise artifacts that have not yet been mitigated in LIGO detector commissioning, such as power mains contamination at 60 Hz and harmonics and environmental vibrations. The simulations discussed below uncovered an error of a factor of 4 in the \({\epsilon }^{2}\)\({m}_{{\rm{A}}}\) sensitivity plot in ref. 15. This error has been corrected in the current study. As GW detectors become more sensitive in the future, one expects a stochastic GW background from compact binary coalescence mergers to emerge eventually, with an integrated broadband stochastic signal perhaps detectable as early as the O4-O5 run26. Nonetheless, the stochastic GW strain power from mergers in a single DPDM search bin will remain negligible for years to come.

The inclusion of a global network of detectors, such as Virgo, KAGRA, and LIGO-India, improves the DPDM search sensitivity, in principle. The degree of improvement depends, however, upon the relative alignments among these detectors as well as their sensitivities. The Virgo detector is currently less sensitive than the two LIGO detectors. In addition, its orientation is not well aligned with those of the LIGO detectors. Future third-generation detectors, such as Einstein Telescope and Cosmic Explorer, will have much lower noise, permitting still more sensitive DPDM searches.

Discussion

In this paper, we present a direct DM search using GW detector strain data. We choose the \(U{(1)}_{{\rm{B}}}\) DP as our benchmark scenario; our early results already improve upon prior searches in a narrow DP mass range, and future searches will probe deeper in DPDM coupling strength and wider in mass range. This first analysis uses a non-templated, single-Fourier-bin cross-correlation detection statistic. Refinements to be examined for analysis of future data sets include multiple DFT coherence times, tuned according to search band, templated filtering over multiple Fourier bins, and exploitation of extremely narrow features expected in the DPDM spectrum, resolvable by GW detectors for loud enough SNR.

With more data to be collected by LIGO and other GW detectors in the coming years, as well as with improved search strategies, we expect DPDM searches to probe steadily deeper in DPDM parameter space. The same strategy can be implemented directly in searches for many other ultralight DM scenarios. The use of data from a GW detector demonstrates the versatility of these remarkable instruments for directly probing exotic physics.

Methods

Simulating the DPDM background

The DPDM background is a superposition of many DP wave functions, similar to the axion DM background as studied in ref. 27. In the galaxy frame, each DP has a random polarization direction isotropically distributed. The magnitude of \({\bf{A}}\) is taken to be a fixed number for each DP particle with normalization discussed below. As for the polarization vector, the velocity direction also follows an isotropic distribution. The magnitude of the velocity is obtained from the Maxwell distribution

$$f(v) \sim {v}^{2}{e}^{-{v}^{2}/{v}_{0}^{2}},$$
(6)

where \({v}_{0}\) is taken to be \(0.77\times 1{0}^{-3}c\)28. In the non-relativistic limit, the polarization vector and the velocity vector are independent of each other.

For the \(i\)-th DP particle, the wave function can be written as

$${{\bf{A}}}_{i}(t,{\bf{x}})\equiv {{\bf{A}}}_{i,0}\sin ({\omega }_{i}t-{{\bf{k}}}_{i}\cdot {\bf{x}}+{\phi }_{i}),$$
(7)

where \({\omega }_{i}=\sqrt{{{\bf{k}}}_{i}^{2}+{m}_{{\rm{A}}}^{2}}\equiv 2\pi {f}_{i}\) and \({{\bf{k}}}_{i}={m}_{{\rm{A}}}{{\bf{v}}}_{i}\). The DPDM background can be generated by superposing many, \(N\), of these wave functions

$${{\bf{A}}}_{{\rm{total}}}(t,{\bf{x}})=\sum _{i=1}^{N}{{\bf{A}}}_{i,0}\sin ({\omega }_{i}t-{{\bf{k}}}_{i}\cdot {\bf{x}}+{\phi }_{i}).$$
(8)

Here the phase of the wave function for each particle, \({\phi }_{i}\), is randomly chosen from a uniform distribution from 0 to \(2\pi\).

To simulate DPDM background, we consider \(N=1{0}^{3}\) DPDM particles. We note that having \(N=1{0}^{3}\) suffices to reveal essential features of the DPDM background, such as coherence length and coherence time. Further, this simulation provides a signal PSD, which agrees well with analytic results (\(N\to \infty\) limit). Thus we believe that the result from this simulation is reliable.

Finally, the normalization of \({{\bf{A}}}_{i,0}\) is determined by the local DM energy density. In the non-relativistic limit, the energy density of DM can be calculated as

$$\frac{1}{V}\frac{1}{T}{\int }_{V}{d}^{3}x{\int }_{0}^{T}dt\ {m}_{{\rm{A}}}^{2}{{\bf{A}}}_{{\rm{total}}}^{2}={\rho }_{{\rm{DM}}}\simeq 0.4\;{{\rm{GeV/cm}}}^{3}.$$
(9)

In order to average out the fluctuations in numerical simulation, the temporal and spatial integrals are taken to be much longer than the coherence time and length, i.e.,  \(T\gg {T}_{{\rm{coh}}}\) and \(V\gg {l}_{{\rm{coh}}}^{3}\). Since the DPDM is obtained from a superposition of \(N\) DP particles in an uncorrelated manner, the total amplitude increases as \(\sqrt{N}\). For a fixed DM energy density \({\rho }_{{\rm{DM}}}\), one expects \(| {{\bf{A}}}_{i,0}| \simeq \sqrt{2{\rho }_{{\rm{DM}}}}/({m}_{{\rm{A}}}\sqrt{N})\), consistent with our numerical results.

Interface to LIGO simulations

We use the LIGO Scientific Collaboration Algorithm Library Suite (LALSuite)29 for mimicking of GW detector response of DPDM and for superposing random Gaussian noise. This suite of programs has been developed over two decades for simulating GW signals, detector response, and for carrying out GW analysis, including source parameter estimation.

Below, we give a brief overview of the relevant LALSuite GW response model and explain what is modified to simulate DPDM-induced effects. When the GW wavelength is much longer than the detector’s characteristic size, i.e., \(\lambda \gg L\), one can use the equation of the geodesic deviation in the proper detector frame to calculate the GW-induced effect,

$${\ddot{\xi }}^{i}=\frac{1}{2}{\ddot{h}}_{ij}^{{\rm{TT}}}{\xi }^{j},$$
(10)

where \(\xi\) is the coordinate of a test object in the proper detector frame. At leading order, the relative change of the arm length is

$$R\equiv \frac{\Delta {L}_{{\rm{x}}}-\Delta {L}_{{\rm{y}}}}{L}={h}^{+}{F}_{+}+{h}^{\times }{F}_{\times },$$
(11)

where \({F}_{+}\) and \({F}_{\times }\) are antenna pattern functions, which can be written as

$${F}_{+}= \,\sum _{i,j}{D}_{ij}{({{\bf{e}}}_{+})}_{ij}=\frac{1}{2}\left[{({{\bf{e}}}_{+})}_{xx}-{({{\bf{e}}}_{+})}_{yy}\right],\\ {F}_{\times }= \, \sum _{i,j}{D}_{ij}{({{\bf{e}}}_{\times })}_{ij}=\frac{1}{2}\left[{({{\bf{e}}}_{\times })}_{xx}-{({{\bf{e}}}_{\times })}_{yy}\right],$$
(12)

with polarization tensors

$${({{\bf{e}}}_{+})}_{ij}= \,{({\bf{X}}\otimes {\bf{X}}-{\bf{Y}}\otimes {\bf{Y}})}_{ij},\\ {({{\bf{e}}}_{\times })}_{ij}= \, {({\bf{X}}\otimes {\bf{Y}}+{\bf{Y}}\otimes {\bf{X}})}_{ij},$$
(13)

and detector tensors

$${D}_{ij}=\frac{1}{2}{({{\bf{n}}}^{{\rm{x}}}\otimes {{\bf{n}}}^{{\rm{x}}}-{{\bf{n}}}^{{\rm{y}}}\otimes {{\bf{n}}}^{{\rm{y}}})}_{ij},$$
(14)

where vectors \({\bf{X}}\) and \({\bf{Y}}\) are the axes of the wave frame, and \({{\bf{n}}}^{{\rm{x}}}\) and \({{\bf{n}}}^{{\rm{y}}}\) are unit vectors along the \(x\) and \(y\) arms of LIGO, respectively.

In order to concretely estimate LIGO’s sensitivity to a DPDM signal, we calculate the DPDM-induced relative change of the arm length as a function of time, i.e., \(R(t)\). Then we inject this as the signal into LALSuite. The background is further added as a Gaussian white noise. As a benchmark, the DP oscillation frequency is set to be \(100/\sqrt{2}\) Hz and \({\epsilon }^{2}\) to be \(5\ \times \ 1{0}^{-44}\). For the simulation, we take \({T}_{{\rm{DFT}}}=1800\) s, \({T}_{{\rm{obs}}}=200\) h, and \(\sqrt{{\rm{PSD}}}=1{0}^{-23}/\sqrt{{\rm{Hz}}}\). The signal appears as a negative real number, i.e., SNR \(\simeq -8\). The sensitivity to \({\epsilon }^{2}\) scales as

$$\frac{\epsilon^{\prime 2}}{{\epsilon }^{2}}=\frac{{\rm{SNR}}^{\prime} }{{\rm{SNR}}}\frac{{T}_{{\rm{coh}}}} {T_{{\rm{coh}}}^{\prime} }\sqrt{\frac{{N}_{{\rm{DFT}}}}{{N}_{{\rm{DFT}}}^{{\prime} }}}.$$
(15)

With this scaling, our simulations are consistent with upper limits shown in Fig. 4 based on the search of O1 data.