Constraining the gravitational coupling of axion dark matter at LIGO

The axion-gravity Chern-Simons coupling is well motivated but is relatively weakly constrained, partly due to difficult measurements of gravity. We study the sensitivity of LIGO measurements of chirping gravitational waves (GWs) on such coupling. When the frequency of the propagating GW matches with that of the coherent oscillation of axion dark matter field, the decay of axions into gravitons can be stimulated, resonantly enhancing the GW. Such a resonance peak can be detected at LIGO as a deviation from the chirping waveform. Since all observed GWs will undergo similar resonant enhancement from the Milky-Way (MW) axion halo, LIGO O1+O2 observations can potentially provide the strongest constraint on the coupling, at least for the axion mass $m_a = 5 \times 10^{-13} - 5 \times 10^{-12}$ eV. Along the course, we also emphasize the relevance of the finite coherence of axion fields and the ansatz separating forward and backward propagations of GWs. As a result, the parity violation of the Chern-Simons coupling is not observable from chirping GWs.

The axion is an important candidate of dark matter. Axions are not restricted to the QCD axion, but a variety of axions are predicted from stringy setups [1]. They are very light pseudo-scalar particles coupling to Chern-Simons terms of some gauge fields F F . Combined with proper cosmological histories, a wide range of axions can be a full dark matter candidate (see e.g. [2]).
However, the axion is very elusive as it couples to standard model particles very weakly, suppressed by its large decay constant f a . Thus, usual direct detection experiments are not sensitive to the axion. A whole new varieties of axion detection experiments and astrophysical probes have been proposed, mainly based on its lightness (due to the pseudo Goldstone nature) and the coherent oscillation (due to the non-relativistic dark matter nature) [3]. They can constrain the axion couplings to photons and electrons, for example through supernova cooling, oscillating electric dipole moments, birefringence of pulsars, quasars, and cosmic microwave background (CMB), and the mixing with the photon inside electron plasma. We refer to [4] for reviews.
But the axion-gravity coupling is relatively weakly constrained. Similarly to the axion-photon aF F coupling, the axion-gravity Chern-Simons coupling can be generically produced [5,6]. Whenever there is a gravitational anomaly, there must exist an associated axion coupling to the gravity. The latest bound on the axion-gravity cou-pling 10 8 km (Eq. (7)) comes from the measurement of frame-dragging effects around the Earth by Gravity Probe B [7].
In the meantime, the chirping gravitational wave (GW) from a binary merger arises as a new tool to probe the Universe. Since it has a well predicted waveform chirping in time and frequency domains in a particular way, even small perturbations to the chirping can be confidently detected. Example studies with dark matter perturbations are [8][9][10][11][12][13], one of which is probing coherently oscillating light dark matter around binary mergers [14].
In this paper, we study how the chirping GW can be perturbed by the coherent axion field as the GW propagates through it. Although the gravitational perturbation is usually very small, a resonant phenomenon can occur when the GW frequency matches with the axion Compton frequency. The resulting signal is a resonance peak in the frequency spectrum.
The resonant phenomenon on the electromagnetic (EM) wave has been studied with various observables. For example, the modification of the EM wave propagating through the coherent axion dark matter field can produce a sharp resonance peak in the frequency spectrum [15][16][17][18][19] or can even produce an echo coming back to us [20]. The resonance can also destabilize axion structures [21,22], possibly leaving some signals in the background or producing an explosive burst [23].
On the other hand, the GW resonance from the coherent dark matter field has not been studied in detail, even though the axion-gravity coupling is well motivated too. Up to our knowledge, the GW resonance was first studied in [24], but it lacks detailed analysis of realistic observables. Our work aims at providing an elaborate analysis for the GW resonance and using it to probe axion-gravity couplings with the LIGO. We will mainly focus on the modification of chirping GWs, but will discuss the instability of axion substructures too. Readers may also refer to [25] for other non-resonant GW observables and [26] for axion-generated GW background.
Our work also improves upon the previous works on the resonance in that we correctly include the finite spatial coherence of the axion field and separate the forward and backward waves. Although similar analyses have been done for EM waves in [16,20], the former ignored the spatial coherence while the latter did not discuss the forward wave. Both treatments are crucial in the LIGO observation, and the absence of parity violation observables is one remarkable consequence.
The paper is structured as follow. We start with a summary of main points and physics of the paper in Sec. II. We derive and solve wave equations in Sec. III, introduce our axion signals on the chirping GW in Sec. IV, and present LIGO bounds and prospects in Sec. V. We provide further details on the resonance with various viewpoints in Sec. VI, and discuss interesting findings on the absence of parity violation in Sec. VII. Then we conclude in Sec. VIII.

II. OVERVIEW
We consider the MW axion halo, which is a highly coherent superposition of axion waves, with a long spatial coherence ∼ 1/m a ∆v (with a small velocity dispersion ∆v) and a much longer temporal coherence (longer than the duration of the chirping GW in the LIGO band). The long coherence stems from the non-relativistic nature (v ∼ ∆v 1) of the axion dark matter.
The coherent (temporal) oscillation can induce resonant enhancement of the chirping GW, when the GW frequency matches with the axion Compton frequency. Since the waveform of the chirping GW is very well predicted, the resonance peak can be detected. It can be further distinguished from accidental noise because all observed GWs will experience a similar phenomenon from the MW axion halo. We found that the correlation of all 11 GW observations at LIGO O1+O2 can provide one of the strongest constraints on the axion-gravity coupling.
The resonant phenomenon is essentially the stimulated decay of axions, although we treat those waves classically. We present several analyses to make sense of the particlelike interpretation of the solution of wave equations that we actually obtain and use.
Other remarkable technical points: The finite (spatial) coherence does impact the signal. Not only does it reduce the enhancement, but it also broadens the frequency width and induces finite timeduration of the resonance peak.
We distinguish forward and backward-going GWs generated from the propagation through an axion halo. First, only forward waves from a binary merger will be observed. Second, backward waves must be generated by the energy-momentum conservation, if forward waves are to be enhanced. Last, mostly only forward and backward waves are generated, which can be understood from a symmetry consideration well inside a halo.
The distinction of forward and backward waves leads to different observable relations of the parity violation. In our case, the parity violation exists only on backward waves, hence not observable. But in existing studies, parity violation was observable because non-resonant regime was considered and/or stochastic waves were considered where the forward/backward distinction is not possible.

III. PROPAGATION THROUGH COHERENT AXIONS
We solve coupled wave equations between axion fields and GWs by using an ansatz suitable for the propagating GW. Then we discuss the solution near a resonance regime with small enhancement.

A. Coupled wave equations
The gravitational Chern-Simons Lagrangian L = α 4 aRR gives the linearized action in the flat background as (ignoring the cosmic expansion) where α is the gravitational Chern-Simons coupling constant, h ij is the metric perturbation and κ = 1/16πG. Varying this action with respect to the metric perturbation h ij gives the wave equation [27,28] (2) We approximate the axion field a to have only time dependence through its Compton oscillation (spatially homogeneous) where a 0 is the complex amplitude (containing the initial phase information). Axions are non-relativistic (a dark matter candidate) so their small kinetic energy contribution to the Compton frequency is neglected. The amplitude a 0 (hence, the energy density) is assumed to be constant in time, as the energy density of the axion field is much larger than that of the chirping GW (Sec. VI A). For a more realistic axion halo spatial profile, see Sec. IV C.
To solve Eq. (2) for the propagating GW in a finite axion halo, we introduce an ansatz for h ij considering 1 1. Plane waves propagating in theẑ direction.
2. Backward wave. The conservation of momentum enforces the generation of backward propagating waves when the forward wave is enhanced 2 . We will distinguish forward and backward waves, in order to describe forward propagating GWs that we eventually observe. This leads to different observable relations from previous works; see Sec. VII.
3. Circular polarization. The Levi-Civita tensor ijk mixes the + and × polarizations, while right handed (R) and left handed (L) circular helicities are decoupled. 1 In Appendix A, we present another approach of solving the wave equation, giving the same result. 2 From a symmetry consideration, the generation of only forward and backward waves must be true, at least well inside a finite halo. But there can be slight leakage over all directions near the boundary of a halo or a coherent patch, although the boundary still varies smoothly over a large scale. We ignore the leakage.
These conditions give the following ansatz (similarly to the photon ansatz introduced in [16] but in the circular polarization basis): where each helicity mode is expressed as 3 where h F and h B are complex amplitudes for the forward and the backward waves, with the superscript s = L, R denotes helicity ands refers to the opposite to s. The polarization tensorê ij is defined with respect to the direction of +ẑ propagation.
Applying Eqs. (3) -(5) into the wave equation (2) (with |ḧ/ḣ| m a ) gives coupled first-order differential equations for forward and backward waves. They become decoupled in the second order equations as where is the coupling parameter that we use to describe the axion-gravity coupling [24,29], is the fractional deviation of k from the resonance frequency m a /2, and a useful dimensionless combination of parameters is , where m 2 a |a 0 | 2 /2 = ρ a is used. For m a = 10 −13 ∼ 10 −10 eV relevant to the LIGO band, γ 1 for the most range of currently allowed coupling and density. Thus, the enhancement rate will be assumed to be small throughout the paper. The origin of the name is clear from Eq. (6) which describes exponential enhancement when µ is real. These parameters will be used widely in our phenomenology study.

B. Solution for finite propagation
The solutions of the wave equation Eq. (6) can be expressed in terms of initial values h where λ (R/L) = +1/ − 1 and φ 0 denotes the phase part of the axion amplitude as a 0 = |a 0 |e iφ0 . The initial condition relevant to the forwardpropagating chirping GW is h These GW solutions are of the same form as those of electromagnetic(EM) waves in [16,22], even though the wave equations are different. These solutions are valid for both real and complex µ's. We hereafter focus only on the forward wave as it is what we observe from binary mergers. This is overlooked in the previous GW work [24]; see Sec. VII for observational implications.
Now consider finite propagation of GW with µt 1 in Eq. (13a), where t is the propagation time. This limit will be relevant to the finite coherent axion patch. We first express Eq. (13a) in the polar form h (s) and express F (t) and ψ(t) up to their lowest order axion contributions under µt 1 (γ, 1 is always assumed). They are where δ 1 by µt 1 and The leading term in the phase ψ, combined with the phase of the ansatz in Eq. (5), gives the phase velocity equal to the speed of light: (m a /2)t + (m a /2) t = kt so that ω = k. Thus, the second term of Eq. (16) gives the correction to the dispersion relation as will be discussed in Sec. V C. Hereafter, we no longer distinguish the wave number k and the angular frequency ω in the leading order. Similarly, for F (t) in Eq. (15), the 1 refers to the original wave and the second term δ( , t) is the enhancement due to the stimulated axion decay. The resonance shape described by Eq. (15) is different from the naive expectation from Eq. (6). This is due to the finite propagation time, or equivalently the finite coherence of the axion field. In the next Sec. IV A, we discuss physical properties of these solutions with finite propagation time, in comparison to those with infinite propagation.

IV. SIGNAL
We introduce two kinds of axion signals, main one in Sec. IV A and another in Sec. IV B. In the last two subsections, we discuss how to calculate them from the propagation through multiple coherent axion patches of a galactic halo.

A. Signal 1: Resonance with finite coherence
One may use the wave equation in Eq. (6) to describe an exponential growth when µ is real for ma 2 (1 − γ) < k < ma 2 (1 + γ) from Eq. (11). The width m a γ is very narrow (see Eq. (10)) so that k ≈ m a /2. This relation is consistent with the particle interpretation of the phenomenon as a stimulated axion decay into two gravitons. Thus, the growth is also called the 'resonant enhancement'. The maximum enhancement rate from Eq. (11) µ max = m a γ/2 is also determined by γ.
However, the finite coherence of the axion field makes important modifications on the resonance. First, the resonance width is broadened, not simply given by m a γ as above. The resonance in each coherent patch is given by Eq. (15) with substituting t by the coherent patch size L coh ∼ 1/m a ∆v ∼ 1/m a v (∆v ∼ v is the velocity dispersion of axions; see Sec. IV C for more details), where δ patch 1 describes the enhancement at each patch. The frequency width of the enhancement is given by the central peak of the sinc function: −2∆v 2∆v, corresponding to m a /2 − m a ∆v ω GW m a /2 + m a ∆v yielding the peak width ∼ 2m a ∆v. This is different from the estimation above, where the width was determined by γ. The modification can be understood from two perspectives. First, the length of the patch is 1/m a ∆v, so each patch cannot have a frequency resolution better than m a ∆v. This determines the resonance width. Another point of view is that axions have the velocity dispersion ∆v, so that the observed resonance width is Doppler broadened by fractionally ∆v. These two perspectives are essentially the same, since the patch size is determined by the velocity dispersion.
In addition, the broadened frequency width m a ∆v implies the time duration 1/m a ∆v of the resonance, related by the Fourier transform. The time duration equals to the size of a coherent patch L coh , hence the time taken for a GW to pass one patch. Such a long duration, combined with the time-delay of a resonance, may affect detection methods as will be discussed in Sec. V C.
Our main signal is a narrow resonance peak in the chirping GW frequency spectrum, produced by the resonant enhancement. We show an example signal in Fig. 1, which results from the propagation through a 100 kpc axion halo consisting of many smaller coherent patches. Although the exact resonance shape depends on finite coherence and effects from multiple patches in Sec. IV D, the basic properties are as discussed above: the narrow peak at k m a /2 and peak width and height determined largely by γ and ∆v.
Since the chirping waveform from binary mergers is very well predicted and does not usually accompany such a sharp peak, the absence of such a peak in the LIGO observations can constrain the resonance enhancement. The peak can be confused with instrumental noise which also often appears as a sharp peak. But since all GWs arriving at us will experience similar enhancements due to MW axion clouds, one can gain confidence by corre-lating all observed GWs in the frequency and time domains simultaneously. Therefore, if a peak is observed in one GW, a peak with similar properties (frequency, timing, and amplitude) must be observed in all GWs. We quantitatively study this signal with LIGO capability in Sec. V B.

B. Signal 2: Explosion
Another constraint comes from the existence of certain dark matter substructures. If µ is too large, the stimulation becomes quicker and quicker so that a coherent axion patch becomes unstable and decays almost entirely into GWs. Such happens when [21] which essentially means that the enhancement rate µ is larger than the passing time within a coherent patch L coh = 1/(m a ∆v) in Eq. (22). Thus, this may happen for small enough ∆v (long enough coherence) and high enough density ρ a .
If there existed such substructures that could explode (satisfying the above condition), they must have almost disappeared by today because there are background photons and GWs everywhere with essentially any frequencies. Produced photons and GWs might have been dissipated enough or became a part of the stochastic backgrounds so that they might not be observable today. Instead, it is the observation of certain dark matter substructures surviving today which can impose an upper limit on the Chern-Simons coupling.
The observed dark matter substructures with possibly the largest enhancement rate are likely dwarf galaxies. They have small velocity dispersion ∆v = O(1 − 10) km 10 −5 (thus, the long coherence length) and large dark matter density ρ 10 3 × 0.3 GeV/cm 3 at its central core [30] albeit some uncertainties. We conservatively use these values to estimate the upper bound on the Chern-Simons coupling, from the existence of dwarf galaxies; in any case, the bound on the coupling is not so sensitive to the density as it scales with ρ 1/4 a in Eq. (10). As shown by the red solid in Fig. 2, this constraint is weaker than that from the resonance peak. Also, the approximation with small µt 1 for the MW axion halo is thus good.
As an aside, which axion substructures could lead to explosion? The axion minicluster [31] has long coherence and high density. Virialized within its Jeans length, the axion minicluster has the Jeans length [4] r J = 2π which is of order of the de Broglie wavelength. Here, ρ mc is the density of a minicluster. Thus, the explosion µ max L coh µ max r J 1 (Eq. (18)) happens when As expected, this value of is much smaller than the bounds coming from dwarf galaxies and resonances (cf. Fig. 2). Although this estimate can be subject to small gravitational redshifts due to the minicluster itself and axion self interactions [22], we conclude that axion miniclusters are irrelevant to our work. If miniclusters had existed, they would have almost disappeared by today by explosion, or the axion coupling is too weak to be probed by any methods.

C. Modeling an axion halo with multiple coherent patches
A realistic axion halo is not infinitely coherent. The coherence property varies among axion dark substructures. As discussed in Sec. IV B, it is good enough to consider an axion halo without miniclusters; such a scenario is motivated by the misalignment production mechanism [32][33][34].
Such axion halo is virialized with the Milky-Way (MW) whose total mass is ∼ 10 12 M in a radius of 100 kpc. The virial velocity v ∼ 10 −3 with the Maxwellian dispersion ∆v ∼ v ∼ 10 −3 leads to the superposition of axion fields with the long spatial coherence length 1/m a ∆v ∼ 1/m a v ∼ 10 3 /m a (with k = m a v) [3]. This length is the size of a coherent patch in that a halo has a spatially oscillating profile with the oscillation length scale of L coh . This is essentially the random-walk superposition of N axion waves (with random phase) which leads to the total amplitude |a 0 | ∝ √ N consistent with the energy density given by ρ a = m 2 a |a 0 | 2 /2 (see Eq. (25) for the value of ρ a ). In addition, ω = m a (1 + v 2 /2) so that the temporal coherence is broken only after a long time 1/m a ∆v 2 ∼ 10 6 /m a , much longer than the GW propagation time in each coherent patch. Thus, we ignore the temporal incoherence while taking into account the spatial incoherence.
Thus, an axion halo is composed of many smaller patches of sizes about the coherence length. As the GW propagates through an axion halo, it passes through multiple coherent patches. As the resonant effect grows only within a coherent patch, the total enhancement will be the sum of the individual patch's effect. We discuss how to sum them up in Sec. IV D.
In Sec. III, we have solved wave equations by assuming the spatially homogeneous and infinite axion field. We apply this solution to each coherent patch, which is actually of finite size and spatially varying. In effect for simplicity, we approximate each coherent patch as a Heaviside profile with the length L coh and the amplitude satisfying ρ a = m 2 a |a 0 | 2 /2. The solution is thus good enough well inside the patch, but our calculation does not include the entrance and exit of GWs through the boundary of a patch. Nevertheless, this approximation can still capture the main physics of the phenomenon. We defer more accurate calculations to the future.

D. Summing effects from multiple patches
A dark matter halo in a galaxy consists of many smaller coherent patches. The resonant enhancement occurs only within a coherent patch. Therefore, we need to sum the effects from each patch.
There is a subtlety here. As discussed in Sec. IV A, the resonance has a time duration of L coh due to the finite frequency width of a resonance. Thus, not all resonance stimulates axion decays simultaneously. It is complicated to account for the fraction of GWs participating in the stimulation at each moment. But it is the original chirping GW which is largest and dominantly stimulating the axion decay; while at later time of the propagation, when the enhanced signal grows larger than the original chirping one, this issue becomes more relevant.
Rather than figuring out an accurate method, we estimate the range of the maximum and minimum possible summation. The enhancement in one patch is 1 + δ(f ) from Eq. (15) (and Eq. (17)), where δ 1 is peaked at the central resonance frequency f 0 . What is the enhancement after passing N patches? The maximum summation assumes that all the axion signals from one patch contribute to the stimulation in the next patch, yielding the maximum total enhancement On the other hand, the minimum summation assumes no axion signals but only original GW stimulates in the next patch. This yields the minimum total enhancement We use these two estimations to obtain an uncertainty band of our estimation (see Fig. 2, for example). A more realistic summation is likely to be between them.
In both cases, the summation depends on the N δ(f ). As δ depends linearly on ρ (δ ∝ γ 2 ∝ ρ), we can use the line-averaged density along the line of sight (LOS), for each coherent patch. For a LOS toward outside the galactic halo, the following line-averaged density is obtained for both NFW and Burkert profiles ρ(r) of the MW dark matter halo, where 100 kpc is the assumed halo radius and 8 kpc is our distance from the MW center. We have taken best-fit parameters for both profiles from [35]. We have checked that 100 kpc 8 kpc ρ(r)dr 0.98 ∞ 8 kpc ρ(r)dr for both profiles, confirming that the 100 kpc radius is sufficient. We use this average density value for ρ a in our numerical study.
Last, the factor N δ makes the importance of finite coherence in yet another manifest way. From Eq. (9) and Eq. (17), we have δ patch (f ) ∝ ρ 4 m 2 a /∆v 2 ×sinc 2 ( /2∆v). For the travel through the MW axion halo of size R, there are N = Rm a ∆v number of coherent patches. So the whole enhancement depends on the combination The linear dependence on the R and ρ is reasonable, and the overall dependence on 1/∆v implies that the enhancement is greater for the longer coherence from the smaller velocity dispersion. Thus, the effect of finite coherence indeed suppresses the size of the enhancement while broadening the frequency width.

V. LIGO BOUNDS AND PROSPECTS
We use 11 GWs observed in LIGO O1 and O2 to obtain constraints on the axion coupling. As discussed in Sec. IV A, every observed GWs will exhibit a common resonance peak due to the MW halo. In this section, assuming that the correlation of GWs can be made to find the common peak, we focus on individual GW properties in estimating the LIGO sensitivities.

A. Detection criteria
We measure the likelihood L of the existence of a resonance peak using the peak strength as where i is summed over all observed GWs. The peak strength ∆SNR peak is defined as the SNR in the resonance region (−2π∆v ≤ ≤ 2π∆v from Eq. (17)) subtracted by the original SNR of the chirping GW; this roughly measures the significance of the deviation from smooth chirping. The ∆SNR peak can be calculated from Eq. (14) or by multiplying Eq. (23) or Eq. (24) to the original waveform.
We require the log-likelihood to be larger than 100: This is the only requirement in our simplified analysis. Although this simple requirement can be mimicked by a strong peak in single GW, in real analysis the correlation of all the GWs (about the resonance shape in both frequency and time domain) must be made for further consistency. Assuming that such correlation can be made, we use the requirement χ ≥ 100 to estimate the LIGO bounds and prospects. The 100 is arbitrary but conservative requirement. The original SNR in the resonance region is ∼ O(1) (for the 11 LIGO observations). Even if a somewhat larger frequency bin is used, ∆SNR peak 10 might be good enough to be confidently detected; the fractional measurement uncertainty of the overall amplitude estimated by the Fisher information matrix is ∼1/SNR [36] so that our requirement is well above this sensitivity. We want to be conservative as real analysis including matched filtering and correlation may bring additional uncertainties. But the conservative estimation can be good enough because the signal strength ∆SNR peak depends on strongly (N δ ∝ γ 2 ∝ 4 from Eq. (26)), thus a mild improvement on the requirement does not bring large improvement on the bound. Therefore, while encouraging a more dedicated analysis, we are content with estimating conservative bounds and prospects based on our simplified analysis; see Sec. V B for the results and Sec. V C for other realistic aspects.

B. Results
In Fig. 2, we show the LIGO bound on the axion Chern-Simons coupling as a function of the axion mass m a (the corresponding peak GW frequency f 0 is shown on the upper horizontal axis). The gray shaded region is excluded, from the absence of a resonance peak in the 11 LIGO observations so far (Signal 1 in Sec. IV A); each GW is considered up to its innermost stable circular orbit. This region is obtained by the most pessimistic summation of multi-patch effects as in Eq. (24). The hatched region indicates ambiguities in the summation method; this is the region that could be excluded if a somewhat more optimistic summation can be used. This region extends to the lower range of obtained by the most optimistic summation in Eq. (23). A more realistic bound may lie somewhere in this band (Sec. IV D). The dot-dashed extension of the bounds are the expected bounds with one more NS-NS observation so that a correlation can be made with existing NS-NS data in the highest frequency range O(1000) Hz; heavier binaries merge at lower frequencies. The existing bound from Gravity Probe B satellite measurement of the frame dragging effect [7] is shown as the horizontal dashed. The bound from the existence of dwarf galaxies (imposing that such systems are not exploded by resonant enhancement) is shown as the red solid (Signal 2 in Sec. IV B). This is weaker than the previous two.
The LIGO bound from the absence of a peak is stronger than the existing established bound, at least for the axion mass range 5 ×10 −13 eV m a 5×10 −12 eV. The bound can be stronger if a more aggressive summation can be used, and the heavier mass range up to m a 5 × 10 −11 eV can be constrained if more NS-NS mergers are observed, as discussed.
How will the bound improve with more data and smaller noise? For example, a 10 times smaller noise (achievable with, e.g., Einstein Telescope) will enhance SNR by 10 and observe 10 3 times more GWs, yielding ≈ (10 × √ 1000) 1/4 4.2 times stronger bound on . Similarly, n times smaller requirement on χ means n 1/4 times stronger bound on . The measurement of lower frequency range from future GW detectors can also provide new constraints on the lower range of m a .
One can also note that the bound becomes stronger for the heavier axion. This is basically because γ ∝ m a for a given axion energy density ρ a (see Eq. (10)), giving (N δ) halo ∝ m 3 a as in Eq. (26). This strong dependence on m a overcomes the frequency dependences of the noise curve and chirping GW spectrum; but slight mass dependence of the bound comes from these.
In all, LIGO is potentially able to improve the bound on the axion-gravity Chern-Simons coupling. We encourage a careful reanalysis of the currently available data.

C. Time-delay of a resonance from dispersion
The resonant enhancement also modifies the group velocity of a resonance peak, delaying the arrival of the peak relative to other frequency parts of chirping GW. As the original chirping GW has almost one-to-one correspondence between the frequency and arrival time, this time-delay produces an observable change of the timedomain waveform of the GW.
The dispersion relation can be obtained from the correction term in Eq. (16) and ansatz (5) for the case of small enhancement (µt 1) in the vicinity of = 0 as This means that ω = k at t = 0 (not enhanced yet) but starts to deviate from k as GW propagates through a coherent patch. Note that the dispersion is parity independent; see Sec. VII for usual parity-dependent dispersion. For k > m a /2, ω decreases from k toward m a /2, and opposite for k < m a /2. As µt grows larger than 1, referring back to more general equation Eq. (13a), we find that the phase converges to some constant which implies (by ansatz Eq. (5)) ω = m a /2 regardless of k. This behavior is approximately understood because the GW produced from axion decays has ω = m a /2 by the energy conservation, while its spatial mode is determined by initial chirping GW with the wavenumber k = m a /2. As the enhancement grows, ω = m a /2 = k dominates a whole GW. Back to Eq. (29) with µt 1, the group velocity of the axion signal v g = dω/dk at = 0 is This again means that ω = k and v g = 1 at t = 0 (not enhanced yet) starts to deviate with t. The group velocity is less than 1 so that the axion signal arrives later than the chirping GW. This could complicate the search because too large time delay will conceal the correlation between the appearance of axion signal and the arrival of chirping GW. Thus we estimate the time delay. The average group velocity during the propagation through one coherent patch (and this is the average group velocity in the galactic axion field) is given bȳ Eq. (31) gives the time delay of the resonance peak with respect to the chirping GW. In Fig. 3, we plot this timedelay contours. In the parameter space that can be probed at LIGO, the time delay is 1 -100 seconds, which is also about the duration of a resonance ∼ 1/m a ∆v = 1 -100 seconds. This time scale is, however, longer than the typical duration (seconds or less) of chirping GWs in the LIGO band. Thus, we need to include longer timeseries of data in order to capture the peak which may not be much time-overlapped with the chirping GW. The contours (solid) of the arrival-time delay of a resonance peak with respect to the original chirping GW, induced from the propagation through a 100 kpc axion halo. The overlaid are the exclusion plots in Fig. 2.

D. Similar bounds on the axion-photon coupling
We briefly comment on the axion-photon coupling. Since the solution of the coupled EM wave equations is in the same form as Eq. (13a) [16,22], we can readily apply the same analysis done here to the photon case. The signal would be the extragalactic EM waves with a common peak. By simply requiring the maximum total enhancement of a single good EM signal to be greater than 10 (as light measurements are more precise), we estimate the bound on the axion-photon coupling to be g aγγ 10 −8 -10 −2 GeV −1 for the axion mass range 10 −11 -1 eV. This is similar or slightly weaker than the laboratory bounds, while much weaker than the Helioscope bound by about 2 -8 orders of magnitudes [37]. We defer more detailed analysis and comparison to a future project.

A. Energy conservation and axion backreaction
The energy conservation implies that the amplitude of the axion field should decrease as the GW amplitude is enhanced. However, dark matter energy density is much greater than any reasonable GW energy density. We can estimate the chirping GW energy density as the following. GW150914 emitted ∼ 3M of energy at 400 Mpc, and by assuming all the energy was released in the last 0.1 second of chirping, we have ρ GW = 60 eV/cm 3 . This is incomparably smaller than the dark matter energy densityρ DM = 0.04 GeV/cm 3 in Eq. (25). Thus, we can assume that axion fields do not decrease in our work.
But if somehow energies of both waves become similar, the coupled wave equations describe the energy transfer between them through the time-evolution of both amplitudes. For example for the EM wave case, Eq.(8) and Fig.2 of [16] show such time-variation. Back to a general point of view, the absence of explicit time dependence of the Lagrangian guarantees the energy conservation for a dynamically-evolving axion field. After all, the backreaction of axion fields will stop the exponential growth of GWs (explosion) at some point.

B. Stimulated axion decay rate
We obtain another insight on the stimulated decay by calculating the axion decay rate from the energy gain of the GW, which equals to the energy loss of the axion. Since the spatially averaged energy density of GW is given by ρ GW = ω 2 GW (|h (R) | 2 + |h (L) | 2 )/64πG, the energy density gain of the forward and backward waves are where the last equality is due to the momentum conservation (this can also be explicitly derived from Eq. (13b)). From the energy density loss of the axion 2×∆ ρ GW F , we obtain the decay rate of the axion as (for µt 1; otherwise, the rate increases exponentially) As it should be, this is proportional to the energy density of the GW and independent on the axion energy density. Remarkably, the form of sin 2 ((ω −ω 0 )t)/(ω −ω 0 ) 2 is nearly identical to the probability of stimulated emission in quantum mechanics [38]. This consideration supports the physical picture of the resonant enhancement as the stimulated decay.

C. Effective 'graviton' mass
Even without the axion-gravity Chern-Simons coupling, the GW experiences a dispersion due to the intervening mass density, similarly to the photon's plasma mass in the electron medium. Following the EM wave case in [21], we check that such effect is negligible for the GW.
The GW refractive index is given by n = 1 + 2πGρ/ω 2 [39]. This gives the dispersion relation ω 2 = k 2 − 4πGρ. The ratio of the matter-induced dispersion to the effect of gravitational Chern-Simons coupling in Eq. (A4) is (34) at the resonance. This is an incredibly small number; for instance, ρ = 1 GeV/cm 3 , m = 10 −12 eV, and = 10 8 km give 6.3 × 10 −30 . Thus, the dispersion of GWs due to intergalactic matter can be ignored.

D. Axions in the source galaxy and intergalactic region
The sharp peak was unambiguously associated with the axion signal because every GWs will exhibit a common peak from the propagation through the Milky-Way axion halo. What about axions in other galaxies (in particular, the one that hosts the source of the GW) and in intergalactic region?
First, the cosmological redshift of a source galaxy varies among different sources. So does the observed peak frequency. Such peaks may not be confidently identified.
The effect from intergalactic dark matter may not be strong enough due to low density and redshift. While a half of total dark matter resides in the intergalactic region, the density there which can be estimated as Ω DM ρ c ∼ 10 −6 GeV/cm 3 is 10 −4 times smaller thanρ a in the MW halo (Eq. (25)). Thus, the enhancement factor N δ ∝ Rρ a /v (Eq. (26)) either needs more than a few Gpc propagation (10 4 times longer than the MW halo size) or much smaller velocity dispersion, to produce similar size of total enhancement. The small dispersion is unlikely; for example, the Local Group velocity with respect to the CMB and the escape velocity of the galaxy are all O(100) km/s 10 −3 . In addition, continuously varying redshift through the intergalactic region will further hinder the generation of a sharp and large axion signal peak. Thus, we ignore the intergalactic contributions.

VII. COROLLARY: ABSENCE OF PARITY-VIOLATION OBSERVABLES ON THE CHIRPING GW
Many previous works have studied the parity-violation signals in photons due to the axion Chern-Simons coupling. However, our solution shows that the parityviolation is not observable in the resonance regime of the chirping GW. They are not in contradiction as we explain in this section.
The parity violation in previous works are manifest in two ways: one is through the dispersion relation [15,25,[40][41][42][43][44] and the other through the enhancement [15,24]. For our case, the parity-dependent dispersion is absent because we consider the resonance regime (k ≈ m a /2), and the parity-dependent enhancement is absent because we consider the forward propagation of waves (not stochastic waves).
First, the parity-dependent dispersion is obtained from the wave equation in the form of Eq. (A4) as ω(k) 2 = k 2 + 2λ (s) m a kγ sin(m a t). (35) This shows the usual parity-dependent (λ (s) -dependent) dispersion relation, making the phase velocity deviate from c and oscillate oppositely for opposite polarizations. But the deviation (the second term) oscillates in time with the frequency 1/m a . In previous works with k m a /2, this oscillation was much slower than the high frequency of photons. But in our case in the resonance regime, they are comparable (k ≈ m a /2) so that the deviation almost averages out in one period of the GW/photon. Instead, near the resonance, there arises the dispersion relation which is parity independent, as discussed in Sec. V C.
Second, the usual parity-dependent enhancement arises from the time evolution of spatial Fourier modes, h(k, t) = h(x, t)e −ikx [15,24]. But these spatial modes are the sum of forward and backward propagating waves. By separating the propagation direction, as in our solution Eq. (13a) and 13b, we find that the parity-violation exists only in backward waves due to the initial condition h B = 0. This effect is not observable in our case because what we observe is only h F . The previous solutions are suitable for stochastic backgrounds, like CMB or stochastic GWs, where waves with all directions are mixed up. Such waves can exhibit the parity violation as polarization-dependent enhancements, as studied in previous works.

VIII. CONCLUSIONS
We have shown that the LIGO observation of chirping GWs can constrain the axion-gravity Chern-Simons coupling, through a resonance peak of the GW induced by the coherently oscillating axion dark matter field. As all the observed GWs will have a peak with common properties (frequency, duration, and height), the correlation among them can confidently detect or reject the peak. We have found that 11 GW observations at LIGO O1+O2 can already provide the strongest bound on the coupling, at least for m a = 5 × 10 −13 ∼ 5 × 10 −12 eV (see Fig. 2). With more LIGO observations, the range can be extended and the bound can be stronger. A careful reanalysis of existing data is encouraged.
The resonance phenomenon is essentially the stimulated decay of the axion. Not only does the resonance condition f 0 m a /2 support this particle-like interpretation, but also the decay probability estimated from the energy gain and loss of the fields agrees with the quantum mechanical description of stimulated emissions and absorptions. This is remarkable as we have never quantized these waves.
The finite coherence of the axion field determines the resonance (axion signal) shape in large part. First, it suppresses the height of the resonance peak in the chirping GW spectrum, while broadening the peak width. The broadening in the frequency domain also makes the signal persist as long as the size of a coherent patch. The resonance-produced axion signal is also time-delayed compared to the original chirping GW, and this timedelay is also affected by the finite coherence. Resonance searches must account for these effects.
A proper ansatz treating forward and backward-going waves separately is crucial for our work. It is because only the forward-going chirping GW can be observed. This is different from the stochastic background of GWs and CMB, where waves with all directions are mixed up. As a consequence, our solution does not exhibit the parity violation from the axion Chern-Simons coupling in the forward wave at the linear order, but this is not in contradiction with previous studies.
Last but not the least, the resonant effect can sometimes become so efficient that an axion substructure may not exist today. This happens when the axion structure has small velocity dispersion (hence, long coherence) and high density. The axion minicluster is one example that might have exploded by today, but the coherent axion field virialized with a whole galaxy does not explode given the current bound on the coupling. Although we assumed that the signal of explosion had diffused away, it would be interesting to study if any observable signals remain.
In all, we have studied one way to probe and constrain the axion-gravity Chern-Simons coupling, which is generic and well motivated. A careful reanalysis of LIGO data may provide one of the strongest constraints on this coupling. Various other types of axion-gravity couplings may also be probed in a similar way. In this appendix, we solve Eq. (2) by another method. By expressing h ij in Eq. (2) as a Fourier transform, the wave equation for each Fourier amplitude becomes [24] h (s) + 4λ (s) γ cos(m a t) 1 + 4λ (s) k ma γ sin(m a t) whereh (s) represents the spatial Fourier mode with wave number k (so different from the amplitude h By following the transformation in [25] with cosmic expansion neglected, we define ψ as and we havë Ψ + k 2 + 2λ (s) m a kγ sin(m a t) Ψ = 0 (A4) which is the ordinary Mathieu equation. The exponential factor in Eq. (A3) is a non-resonant term since the exponent oscillate with small amplitude. Any analytic method solving the Mathieu equation tracks only the resonant term, so this factor does not affect the result. Also physically, this factor will be canceled out in average, due to its dependence on the relative phase between the axion field and the gravitational wave (for such cases, the argument of cos's and sin's should have a constant phase term, like mt + ψ 0 ). To solve Eq. (A4) by the two variable expansion method [45,46] we look at the solution behavior near the resonance at k = m a /2. To do this, we use (4γ) for an expansion parameter, and the expansion will done up to the first order. The two variables in the expansion are the ordinary time ξ = t and the slow time η = 4γt. ξ is the time scale of wave oscillation, while η is the time scale of amplitude change. We regard Ψ as a function of the two independent variables, as Ψ = Ψ(ξ, η). And the time derivative operator becomes Similarily, k and Ψ are expanded as k = m 2 1 + 1 (4γ) + 2 (4γ) 2 + · · · (A6) and Ψ(ξ, η) = Ψ 0 (ξ, η) + (4γ)Ψ 1 (ξ, η) + (4γ) 2 Ψ 2 (ξ, η) + · · · . (A7) By putting Eqs. (A5)-(A7) into Eq. (A4), we can obtain series of equations assorted by the order in (4γ). The 0th order equation is and the 1st order equation is (note that ξ = t) Eq. (A8) gives Ψ 0 in the form of and putting this into Eq. (A9) with using trigonometric identities gives where the higher frequency terms in RHS whose resonance appear only in higher orders are not shown. This gives the slow flow equations for A and B as We put these to Eq. (A10) and recover the notations we used in the main paper. Since was used to denote the fractional difference between k and m a /2, we have 1 = /(4γ) (see Eq. (A6)). Then, we have νη = ±µt and the solution we obtained is (recall ξ = t) Ψ(t) = C 1 cos( m a 2 t) − β (s) sin( m a 2 t) e µt where the arbitrary constants are redefined. Recalling that we started from Eq. (A1) about the spatial Fourier mode ∼ e ikx , the e i ma 2 t part denotes the backward wave and the e −i ma 2 t part is for the forward wave. We now put the initial condition of vanishing backward wave at t = 0. We first write the D coefficients as D 1 = D(1 − iβ (s) ) and D 2 = −D(1 + iβ (s) ), and normalize by the initial amplitude of the forward wave h 0 = 2D(1 − (β (s) ) 2 ). This gives Ψ(t) = h 0 1 + (β (s) ) 2 1 − (β (s) ) 2 sinh(µt)e i ma 2 t +h 0 cosh(µt) − i 2β (s) 1 − (β (s) ) 2 sinh(µt) e −i ma 2 t . (A19) To further simplify the coefficients, we recall Eq. (A17). Since λ (s) = ±1, we have This is the same solution in Eqs. (13a) and (13b), considering the phase factors in the ansatz, Eq. (5). Since we assumed initial phase of the axion field to be zero, it does not appear here. Also, here the initial amplitude of the forward wave h 0 is assumed to be real. Thus, we have obtained identical solution via solving the Mathieu equation.