Searching for Dark Photon Dark Matter in LIGO O1 Data

TL;DR

This paper reports a direct search for ultralight dark photon dark matter using LIGO O1 data, achieving sensitivity to specific coupling strengths and mass ranges, bridging gravitational wave detection and particle physics.

Contribution

It introduces a novel method to search for dark photon dark matter using gravitational wave data, providing the first direct constraints in this context.

Findings

Achieved sensitivity to squared coupling of ~4 x 10^-45

Set constraints on dark photon mass around 4 x 10^-13 eV

Demonstrated potential for improved future searches

Abstract

Dark matter exists in our Universe, but its nature remains mysterious. The remarkable sensitivity of the Laser Interferometer Gravitational-Wave Observatory (LIGO) may be able to solve this mystery. A good dark matter candidate is the ultralight dark photon. Because of its interaction with ordinary matter, it induces displacements on LIGO mirrors that can lead to an observable signal. In a study that bridges gravitational wave science and particle physics, we perform a direct dark matter search using data from LIGO's first (O1) data run, as opposed to an indirect search for dark matter via its production of gravitational waves. We demonstrate an achieved sensitivity on squared coupling as $\sim 4 \times 10^{-45}$, in a $U(1)_{\rm B}$ dark photon dark matter mass band around $m_{\rm A} \sim 4 \times 10^{-13}$eV. Substantially improved search sensitivity is expected during the coming years of continued data taking by LIGO and other gravitational wave detectors in a growing global network.

Figures

Equations

$$A_{\mu }\simeq A_{\mu ,0}\cos [m_{\mathrm{A}}t-\mathbf{k}\cdot \mathbf{x}+\theta ].$$

$${\mathbf{a}}_i(t,{\mathbf{x}}_i)\simeq ϵe\frac{q_{\mathrm{D},i}}{M_i}{\partial }_t\mathbf{A}(t,{\mathbf{x}}_i).$$

$$S_j=\frac{1}{N_{\text{DFT}}}\sum _{i=1}^{N_{\text{DFT}}}\frac{z_{1,ij}{z}_{2,ij}^{\ast }}{P_{1,ij}{P}_{2,ij}}.$$

$$\displaystyle\sigma^{2}_{j}=\frac{1}{N_{\text{DFT}}}\bigg{\langle}\frac{1}{2P_{1,j}P_{2,j}}\bigg{\rangle}_{N_{\text{DFT}}},$$

$${\mathrm{SNR}}_j\equiv \frac{S_j}{{\sigma }_j}.$$

$$f(v)\sim v^2{e}^{-v^2/v_0^2},$$

$${\mathbf{A}}_i(t,\mathbf{x})\equiv {\mathbf{A}}_{i,0}\sin ({\omega }_{i}t-{\mathbf{k}}_i\cdot \mathbf{x}+{\varphi }_i),$$

$${\mathbf{A}}_{\text{total}}(t,\mathbf{x})=\sum _{i=1}^N{\mathbf{A}}_{i,0}\sin ({\omega }_{i}t-{\mathbf{k}}_i\cdot \mathbf{x}+{\varphi }_i).$$

$$\frac{1}{V}\frac{1}{T}{\int }_V{d}^{3}x{\int }_0^{T}dt{m}_{\mathrm{A}}^2{\mathbf{A}}_{\text{total}}^2={\rho }_{\text{DM}}\simeq 0.4{GeV/cm}^3.$$

$${\ddot{\xi }}^i=\frac{1}{2}{\ddot{h}}_{ij}^{\mathrm{TT}}{\xi }^j,$$

$$R\equiv \frac{\Delta L_{\text{x}}-\Delta L_{\text{y}}}{L}=h^{+}{F}_{+}+h^{\times }{F}_{\times },$$

$$F_{+}$$

$$F_{\times }$$

$$({\mathbf{e}}_{+}{)}_{ij}$$

$$({\mathbf{e}}_{\times }{)}_{ij}$$

$$D_{ij}=\frac{1}{2}({\mathbf{n}}^{\text{x}}\otimes {\mathbf{n}}^{\text{x}}-{\mathbf{n}}^{\text{y}}\otimes {\mathbf{n}}^{\text{y}}{)}_{ij},$$

$$\frac{{ϵ}^{\prime 2}}{{ϵ}^2}=\frac{{\mathrm{SNR}}^{\prime }}{\mathrm{SNR}}\frac{T_{\text{coh}}}{T_{\text{coh}}^{\prime }}\sqrt{\frac{N_{\text{DFT}}}{N_{\text{DFT}}^{\prime }}}.$$

Full text

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 10^{-22}$ eV Hu:2000ke ; Marsh:2013ywa ; Bozek:2014uqa ; Hui:2016ltb , or it may arise from stellar-mass objects, such as primordial black holes Carr:2009jm .

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 mechanism Nelson:2011sf ; Arias:2012az ; Graham:2015rva , parametric resonance production or tachyonic instability of a scalar field Co:2018lka ; Agrawal:2018vin ; Bastero-Gil:2018uel ; Dror:2018pdh , or from the decay of a cosmic string networkLong:2019lwl .

It was recently proposed in Graham:2015ifn ; Pierce:2018xmy 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 reported LIGOScientific:2018mvr . 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 $10^{-13}\sim 10^{-11}$ eV. The DPDM is an oscillating background field, for which the rest-frame oscillation frequency satisfies: $f_{0}\approx\left[{m_{\rm A}\over 10^{-12}\>{\rm ev}}\right]$(241 Hz). We assume the DP is the gauge boson of the gauged $U(1)_{\text{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)_{\text{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 10^{-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(10^{-6})$ deviations.

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

[TABLE]

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_{\text{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_{\text{obs}}<t_{\text{coh}}\simeq\frac{4\pi}{m_{\rm A}v_{\text{vir}}^{2}}$ and $l_{\text{obs}}<l_{\text{coh}}\simeq\frac{2\pi}{m_{\rm A}v_{\text{vir}}}$. For example, if the DP field oscillates at 100 Hz, the coherence time is only $10^{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 planewaves 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:

[TABLE]

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)_{\text{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 Dror:2017ehi impose very strong constraints on the gauged $U(1)_{\text{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)_{\text{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)_{\text{EM}}$ coupling constant. We emphasize that we choose DP to be a $U(1)_{\text{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)_{\text{B}-\text{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 Estimation. We approximate the DPDM field as a planewave within a coherence region. For a DP field oscillating at frequency $O(100)$ Hz, the coherence length is $O(10^{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_{\text{DFT}}$, with a total observing time $T_{\text{obs}}=N_{\text{DFT}}T_{\text{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 (PSD) for two interferometers are related to the raw powers as ${\rm PSD}_{1(2),ij}=2P_{1(2),ij}/T_{\text{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 disturbances Covas:2018oik , which are excluded from the analysis. Detailed descriptions of broadband LIGO noise contributions may be found in TheLIGOScientific:2016zmo , 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

[TABLE]

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

[TABLE]

where $\langle\ \rangle_{N_{\text{DFT}}}$ denotes an average over the $N_{\text{DFT}}$ DFTs, which may have slowly varying non-stationarity. The signal-to-noise ratio (SNR) can be defined by

[TABLE]

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 SNRj 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 Pierce:2018xmy is $\Delta f/f=10^{-6}$, based on a Maxwell velocity distribution. The binsize in frequency space is set by $T_{\text{DFT}}$, i.e., $\Delta f=1/T_{\text{DFT}}$, which is optimal at only $f_{\text{opt}}\simeq 10^{6}/T_{\text{DFT}}$. For a frequency higher than $f_{\text{opt}}$, the relative frequency binning is finer, implying loss of signal power in single-bin measurements. At frequencies lower than $f_{\text{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. Lentz:2017aay ; Necib:2018iwb . 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 2f_{\text{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.

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. Fig. 2 shows the resulting efficiencies (power fractions) for $T_{\text{DFT}}$ set to be 1800s. 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.

Data selection and analysis. The strain data used in this analysis was downloaded from the Gravitational Wave Open Science Center (GWOSC) web site Vallisneri:2014vxa and transformed to create 1786 1800-second 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 gravitational wave 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 trials 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 gravitational wave 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 Covas:2018oik , 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.

Fig. 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 greater than 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 greater than 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.

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 (CL) 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) formalism Feldman:1997qc 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=10^{-6}$ for the same 893-hour 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 binsize. The mean achieved upper limit is generally worse than the optimal sensitivity because, with fixed binsize 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$\ddot{\textrm{o}}$t-Wash group based on Equivalence Principle tests using a torsion balance Su:1994gu ; Schlamminger:2007ht . 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 10^{-13}$ eV.

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 Pierce:2018xmy . 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 run LIGOScientific:2019vic . 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 gravitational wave detector strain data. We choose the $U(1)_{\text{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 gravitational wave 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. This novel use of data from a gravitational wave 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 Foster:2017hbq . 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

[TABLE]

where $v_{0}$ is taken to be $0.77\times 10^{-3}c$ Smith:2006ym . 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

[TABLE]

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

[TABLE]

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=10^{3}$ DPDM particles. We note that having $N=10^{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

[TABLE]

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_{\text{coh}}$ and $V\gg l_{\text{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_{\text{DM}}$, one expects $|{\bf A}_{i,0}|\simeq\sqrt{2\rho_{\text{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) LAL 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,

[TABLE]

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

[TABLE]

where $F_{+}$ and $F_{\times}$ are antenna pattern functions, which can be written as

[TABLE]

with polarization tensors

[TABLE]

and detector tensors

[TABLE]

where vectors ${\bf X}$ and ${\bf Y}$ are the axes of the wave frame, and ${\bf n}^{\text{x}}$ and ${\bf n}^{\text{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 10^{-44}$. For the simulation, we take $T_{\text{DFT}}=1800$ s, $T_{\text{obs}}=200$ hours and $\sqrt{\rm PSD}=10^{-23}/\sqrt{\rm{Hz}}$. The signal appears as a negative real number, i.e., SNR $\simeq-8$. The sensitivity to $\epsilon^{2}$ scales as

[TABLE]

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

Data availability

The upper limits in Fig. 4 can be found at https://doi.org/10.5281/zenodo.3525343 and the LIGO O1 data that support the findings of this study are available from the GWOSC website https://www.gw-openscience.org Vallisneri:2014vxa .

Code availability

The source code for the analysis is available from the corresponding author upon request.

References