Prospects for axion searches with Advanced LIGO through binary mergers

The observation of gravitational waves from a binary neutron star merger by LIGO/VIRGO and the associated electromagnetic counterpart provides a high precision test of orbital dynamics, and therefore a new and sensitive probe of extra forces and new radiative degrees of freedom. Axions are one particularly well-motivated class of extensions to the Standard Model leading to new forces and sources of radiation, which we focus on in this paper. Using an effective field theory (EFT) approach, we calculate the first post-Newtonian corrections to the orbital dynamics, radiated power, and gravitational waveform for binary neutron star mergers in the presence of an axion. This result is applicable to many theories which add an extra massive scalar degree of freedom to General Relativity. We then perform a detailed forecast of the potential for Advanced LIGO to constrain the free parameters of the EFT, and map these to the mass $m_a$ and decay constant $f_a$ of the axion. At design sensitivity, we find that Advanced LIGO can potentially exclude axions with $m_a \lesssim 10^{-11} \ {\rm eV}$ and $f_a \sim (10^{14} - 10^{17}) \ {\rm GeV}$. There are a variety of complementary observational probes over this region of parameter space, including the orbital decay of binary pulsars, black hole superradiance, and laboratory searches. We comment on the synergies between these various observables.


Introduction
The importance of the recent direct detection of gravitational waves from black hole and neutron star mergers can hardly be over-emphasized [1][2][3][4]. These observations have confirmed the existence of gravitational waves and black holes, among the most important predictions of General Relativity. The binary neutron star event, GW179817, with the coincident electromagnetic observations, has also yielded insight into the nature of short gamma-ray bursts (GRB) and the production of heavy elements in the Universe. What new discoveries might be on the horizon?
It is clear that existing and future gravitational wave (GW) observatories will enable us to learn a great deal about astrophysics [5]. Some of the expected highlights include insight into production mechanisms from population statistics and constraints on the structure of neutron stars from the GW waveform associated with the end stages of inspiral (where tidal effects become important) and the post-merger phase (where a hypermassive neutron star can form). However, the measurement of gravitational waves from binary mergers also provides an unprecedented opportunity to search for fundamental interactions and particles beyond those of the Standard Model of particle physics and General Relativity; see e.g. [6] for a summary. Several examples include

• Self-interactions beyond the Einstein-Hilbert action
It is important to understand how the gravitational sector might be modified by new selfgraviton interactions. The possible form of self-graviton interactions is strongly limited by diffeomorphism invariance, as well as causality and analyticity arguments [7]. The effect of additional graviton interactions in the Cosmic Microwave Background (e.g. [8]) and on compact binary mergers (e.g. [9]) has received some attention in the literature. However, there are a number of challenges associated with the well-posedness of timeevolution in the non-linear theory (e.g. [6,10,11]) which thus far precludes a full picture of binary mergers. New diffeomorphism invariance breaking graviton interactions can also be introduced to modify binary dynamics as well as gravitational wave propagation. Various forms of massive gravity theories can be tested through its modification of graviton dispersion relations [12][13][14][15][16][17][18]. However, none of the massive gravity theories can predict a calculable change to the generation mechanism of gravitational waves due to various instabilities [18,19].
• Exotic compact objects Several proposals exist for compact objects made out of new particles, for example boson stars/axion stars (see Ref. [20] for a review). These new compact objects will have masses and sizes (compaction) different from black holes and neutron stars. The measurement of gravitational and electromagnetic radiation resulting from the merger of such objects provides one means of constraining the associated new particles and interactions (see e.g. [6]).
• Light states coupled to gravity Scalar fields are a ubiquitous feature of physics beyond the Standard Model of particle physics and many extensions of General Relativity. Light scalars that couple to gravity can be probed by black hole superradiance [21][22][23]. In this case, the large gravitational field in the proximity of black holes and their rapid rotation can source the clustering of large numbers of light bosons, which in turn extract angular momentum from the black hole. Indirect observations of the spin distribution of black hole binaries by Advanced LIGO will shed light on the existence of these light particles [23], and searches for continuous wave signals at Advanced LIGO and future gravitational wave detectors might observe these light particles directly [22].
• New force mediator If coupled to matter, light scalars can mediate new long range interactions between compact objects, commonly termed "fifth forces" (see [24] and references within). These interactions have been constrained by laboratory experiments [25][26][27] as well as astronomical observations of the solar system (e.g. [28]) and beyond (e.g. [29]). Laboratory experiments constrain universally coupled fifth forces to be much weaker than gravitational strength if the force has a range that is longer than a few micron [25][26][27]. New scalar forces that arise only in a strong gravity or high density environment, however, are unconstrained and can be looked for with Advanced LIGO.
In this paper, we focus on this last category, building upon previous work [30,31] suggesting that binary neutron star (NS-NS) and neutron star-black hole (NS-BH) mergers can provide powerful new probes of light scalar force mediators. In particular, we assess the sensitivity of advanced gravitational wave detectors, such as Advanced LIGO and VIRGO, to the effects of axions on the GW waveform in binary mergers. Before proceeding, we review the properties of axions.
The QCD axion was originally introduced as a solution to the strong CP problem [32][33][34][35]. Experimental searches for a neutron electric dipole moment (EDM) suggest that the strong CP angle is much smaller than 10 −10 [36], while CP angles in the Cabibbo-Kobayashi-Maskawa (CKM) matrix have been measured to be O (1). The puzzling smallness of the strong CP angle can be resolved by introducing the axion particle a with the coupling where g s is the strong coupling constant, G µν is the gluon field strength withG µν = 1 2 µνρσ G ρσ its dual, and f a the axion decay constant. At low energies, the axion field a gets a potential from the coupling to gluons where m π and f π are respectively the pion mass and decay constant, and m u,d stands for the mass of the up, down quarks. The mass of the QCD axion is related to the axion decay constant by m a = 5.7 × 10 −12 eV 10 18 GeV fa , with m a 10 −12 eV if we require f a m Pl [37]. Recently, it was suggested that if the axion sector has a discrete shift symmetry, the potential of the axion can be much shallower, and the axion mass can be exponentially small [38], opening up the parameter space over which one should search for a QCD axion. In addition to its coupling to gluons, the QCD axion can have model-dependent couplings to photons and derivative couplings to standard-model fermions (see e.g. [39]).
There are a number of other motivations for considering pseudo-scalar particles with many of the same properties as the QCD axion, typically referred to as axion-like particles (ALPs). In the following, we will generally refer to ALPs as "axions", which can have any mass m a and decay constant f a as well as any subset of the interactions possessed by the QCD axion. For example, string theory compactifications generally predict a number of light axions [40]. Axions might be the dark matter [37] (or comprise a significant fraction of it) or provide a candidate for dynamical dark energy [41].
Axions have been constrained by various experiments through their couplings to photons. The axion dark matter experiment (ADMX) published the first constraint on the QCD axion parameter space, assuming that the axion is the dark matter of the Universe in the µeV mass range [42]. Many experimental searches for axions through their couplings to nucleons, electrons, and photons have recently been proposed to cover a much wider range of masses and coupling strengths [43][44][45][46][47][48][49][50][51]. Beside laboratory measurements, indirect measurements of energy loss and energy transport in various astrophysical objects, for instance SN1987 [52], have set the most stringent constraint on QCD axions in the large mass/strong coupling regime. One can also derive constraints on axions from black hole superradiance [22], while for axions with a nuclear coupling one may impose constraints from the measurement of the CP properties of nearby stellar objects [30].
For axions with a nuclear coupling of the form Eq. (1.1), it has been shown [30] that axions can be sourced by compact objects with a high nuclear density, such as neutron stars, thus endowing compact objects with a scalar charge. Such a scalar charge has important implications for NS-NS or NS-BH binary mergers, leading to the emission of axion radiation and an axion-mediated fifth force. Preliminary estimates of these effects on the orbital dynamics and GW waveform were presented in Refs. [30,31], which demonstrated that in principle there can be a significant and detectable effect to target.
For theoretical predictions to match the exquisite data quality of GW179817 and future detections, it is necessary to understand the importance of relativistic corrections, which are typically characterized by the post-Newtonian (PN) expansion in v 2 , the characteristic velocity associated with the orbit. A very useful tool for developing waveforms to arbitrary order in the PN expansion is the EFT framework developed by Goldberger and Rothstein [53,54]. The EFT framework has been used to calculate post-Newtonian corrections to the gravitational potential and quadrupole moments of binary systems, and, as a result, GW waveforms (see [54] for a review), as well as new observable effects beyond General Relativity [9]. One of the merits of an effective field theory approach is that it can be easily extended to include new degrees of freedom and new interactions. In this paper, we extend the effective field theory of gravity for binary systems to include couplings to an axion, and calculate at next-to-leading order (e.g., to 1-PN order) the axion forces between neutron stars as well as axion radiation, both of which will be crucial for computing GW waveform corrections from new interactions. Our result also applies to theories which include an additional massive scalar degree of freedom coupled to neutron stars. To our knowledge, this result for a massive scalar does not appear elsewhere in the literature.
A principle additional result of this paper is a forecast demonstrating the potential for Advanced LIGO to look for massive scalars, and in particular axions, with an event similar to GW179817. We find that Advanced LIGO is a very sensitive probe of the scalar charges of neutron stars and the range of the scalar force (or equivalently, the mass of the axion). Translating this into constraints on {f a , m a } for axions, we find that a GW179817-like event could look for axions in a large region of the parameter space of significant theoretical interest. This region of parameter space is also the focus of efforts by binary pulsar measurements, black hole superradiance, and laboratory experiments, opening the window for interesting joint analyses. In the optimistic scenario of a detection, these other efforts would provide a means for an independent confirmation of the existence of a new fundamental particle. We hope that this analysis motivates a systematic observational effort to constrain axions with existing and future detections by Advanced LIGO and VIRGO, as well as with next-generation gravitational wave detectors.
The paper is organized as follows. In Section 2, we summarize the main effects discussed in [30] and discuss qualitatively the observable consequences of axions on GW waveforms emitted during binary mergers. In Section 3, we adapt the EFT framework to analytically calculate the corrections to the GW waveform from axion mediated forces and axion radiation, and in Section 4, we forecast the observable reach with Advanced LIGO at design sensitivity. In Section 5, we conclude and discuss the implications of an Advanced LIGO discovery or exclusion.

Neutron stars with axions
In this section, we summarize the main properties of axions discussed in [30], which lead to axion mediated forces as well as axion radiation. The coupling of the axion that we search for is the axion nuclear coupling in Eq. (1.1). At low energies, and when the axion is the dark matter, this coupling induces an oscillating electric dipole moment of the nucleus, and has been used to look for dark matter axions in the CASPEr experiment [47]. Note, however, that our setup does not require the axion to be the dark matter. It was recently suggested [30] that for axions with non-vanishing nuclear coupling, there are corrections to the axion potential when the nucleon number density is non-zero. Neutron stars, and to a lesser extent, white dwarfs and stars, can have large enough nucleon number density to significantly change the shape of the axion potential. Over a wide range of axion parameter space, these corrections can lead to phase transitions in large and dense objects, like neutron stars, implying new constraints on the axion parameter space, and providing new opportunities to look for such axions in Advanced LIGO and future gravitational wave experiments. The axion will become tachyonic at a = 0 inside the neutron star in the parameter space where the axion mass in vacuum (m a ) satisfies the condition where the parameter σ N ≈ 59 MeV can be determined from Lattice simulations and n N stands for the number density of neutrons inside a neutron star. For axions satisfying the condition with R NS denoting the radius of the neutron star (NS), the axion will develop a field profile inside/outside of the neutron star approximately given by The axion potential outside the neutron star has minima at a = 0, 2πf a , . . ., and therefore this profile connects the inside of the star, where a = ±πf a , to the local minimum of the potential (in vacuum) at a = 0. The axion profiles of neutron stars in a binary will interact, leading to changes in the strength of the long range interaction, and therefore modifying the power radiated in gravitational waves. Before we delve into details about the calculation of the GW waveform both analytically and numerically, we will first summarize the main observable effects and where they come from.
Axion mediated forces The axion mediates a force between neutron stars when the axion Compton wavelength λ a = 1/m a is comparable to, or larger than, the separation between neutron stars. The force is at leading order, where Q = ±4π(πf a R NS ) are the scalar charges of the neutron stars. Such a force can be either attractive or repulsive, depending on whether the axion field value is the same or opposite sign on the surface of the two neutron stars, respectively. Such a force can be of comparable strength to the gravitational force when the axion decay constant f a is comparable to the Planck scale m Pl . The existence of such a short range interaction can significantly modify the orbital motion of the neutron stars, and therefore the gravitational waveform. At short distances (r R NS ), the axion mediated force deviates from the inverse square law due to the induced axion charges and dipole moments of the neutron stars, which can also change the gravitational wave waveform.
Axion radiation The other major observable effect comes from axion Larmor radiation during the inspiral. The axion radiation turns on when the orbital frequency of the inspiral becomes larger than the mass of the axion. The total power radiated in a neutron star binary inspiral has the form at leading order, where µ = M 1 M 2 M 1 +M 2 is the reduced mass of the system, Ω is the orbital frequency and r denotes the distance between the two neutron stars. The axion radiation is sourced primarily by a time-dependent scalar charge dipole while the gravitational wave radiation is sourced primarily by a time-dependent mass quadrupole. The axion radiation has a weaker frequency dependence when ω m a compared to the gravitational wave radiation, and therefore it is more important at longer distances compared to gravitational wave radiation. Observationally, this implies a specific change to the GW waveform. In absence of an axion force -for instance such a force does not exist at leading order in a NS-BH merger -one gets an additional contribution to df /dt that scales as at leading order, compared to df /dt ∝ f 11/3 for gravitational radiation.
In the following, we discuss in more detail how to calculate corrections to the gravitational wave waveform due to axion mediated forces and axion radiation, and then present how one can use the observation of binary mergers by Advanced LIGO/Virgo in order to constrain the axion parameter space. We will consider NS-NS mergers as well as NS-BH mergers and make use of the phenomenological parameters defined below. The charge of the individual compact objects, which determines the size of the axion mediated force, is for a black hole. (2.7) The dipole moment of the system, which determines the axion radiation, is where r 12 = r 1 − r 2 is a vector that points from charge Q 2 to charge Q 1 , and P = | P | is the magnitude of the dipole moment. In the case of a NS-BH merger, the axion mediated force is zero and the axion radiation is non-zero, while for a NS-NS merger, both the axion mediated force and axion radiation can be present.
To gain a qualitative understanding of the effects of axions on the GW waveform, we show a cartoon plot of the strain versus time in Fig. 1. The effect of axions on the waveform is negligible at times before the objects in the binary are separated by roughly a Compton wavelength of the axion. As the orbit decays within the Compton wavelength, scalar radiation can become an important source of orbital energy loss, especially for large Compton wavelengths. This has the effect of increasing the frequency of the GW, and hastening the merger. Scalar radiation will be present both for NS-NS and NS-BH binaries. For NS-NS binaries, the effect of the scalar force also becomes important once the orbit has decayed to within the Compton wavelength, and can has a strong effect on the orbital dynamics up to the merger. For neutron stars with the same sign scalar charge, the scalar force is attractive, increasing the frequency of the GW and hastening the merger. For neutron stars of the opposite sign scalar charge, the force is repulsive, decreasing the frequency of the GW, and delaying the merger. In the next section, we discuss these effects in more detail.

The effect of axions on binary systems
In this section, we will study the effects of the axion field on the inspiral GW waveform. The inspiral dynamics are usually studied using a PN expansion, in which solutions of the  Figure 1: Schematic plot of the strain versus time for a GW waveform emitted during a binary merger in the presence of an axion. The arrows indicate whether the effects of the axion hasten or delay the merger, and therefore shorten or lengthen the chirp (and increase or decrease its pitch, respectively).
Einstein equations are expanded in the characteristic velocity of the system v. The inspiral waveform can be obtained to a very high accuracy provided the inclusion of sufficiently high order terms. The PN equations of motion can be derived using different methods, all of which lead to the same results at the same PN order. In this paper, we utilize the EFT approach proposed in [53]. We will first review the properties of the EFT and then generalize it to include axions.
A neutron star binary simultaneously involves many scales, i.e. the size of the neutron star R NS , the separation between two neutron stars r, and the wavelength of the emitted gravitational waves λ GW . During inspiral, these three scales have size R NS r λ GW and are related to the velocity through R NS /r ∼ r 2 /λ 2 GW ∼ v 2 1. The smallness of v during the inspiral phase allows us to calculate PN corrections with EFT methods order by order.
To obtain an EFT in the infrared (IR), one can write down an action with all possible terms that respect the symmetries of the system. For example, to calculate the instantaneous potential between binary neutron stars, we represent the neutron stars by two point-like particles, while the mass, spin, and finite size effects of the neutron stars are encoded in the series of couplings between gravitons and the particle world-lines. The value of these couplings can be obtained by utilizing a series of "matching conditions": comparing the physical quantities, for example, the Newtonian potential, calculated with an EFT approach to the quantities one can directly compute easily in the ultraviolet (UV) limit (e.g., General Relativity).
An infrared EFT can also be obtained by "integrating out" the heavy degrees of freedom in the ultraviolet EFT. Specifically, for a neutron star binary, off-shell gravitons mediating long range interactions between two neutron stars (potential gravitons) typically carry momentum k ∼ 1/r ω, while on-shell gravitons that are emitted by the binary (radiation gravitons) typical carry momentum k ∼ ω ∼ v/r and are therefore "lighter" than potential gravitons.
The effective action of the low energy radiation graviton can be obtained by integrating out the "heavy" potential graviton: whereh µν denotes the radiation gravitons, H µν stands for the potential gravitons.
The EFT approach has the advantage of manifesting power counting in the expansion parameter of the theory, which in the case of neutron star binaries is precisely the relative velocity of the binary neutron stars v, and therefore makes it easier to track the PN order. As demonstrated in [53], in the EFT framework, the instantaneous potential as well as the gravitational radiation can be systematically calculated to any order in the PN expansion by including the relevant couplings and "matching conditions", and working out the corresponding Feynman diagrams.
In the following, we will consider binaries consisting of two scalar charged neutron stars or one scalar charged neutron star and a black hole. Similar to the case of pure gravity, we will first write down a series of operators which encode the interactions between the axion and the neutron stars. In particular, we will include operators that characterize the charges and induced dipole moments of the neutron stars. We will then calculate the axion mediated forces between the neutron stars, and utilize several matching conditions to determine the couplings in the EFT. We then treat the effects caused by the axion field perturbatively, and calculate the leading order effects of the axion field on the 1PN Newtonian potential, as well as on the 1PN gravitational radiation. In the EFT with a scalar, as we will demonstrate, we can treat both the scalar charge and the orbital velocity as separate expansion parameters and keep the leading corrections in each.
We will consider two scalar charged neutron stars with mass M 1 and M 2 and charges Q 1 and Q 2 , and with their positions being x 1 and x 2 respectively. As usual, we define and work in the center of mass frame defined at the corresponding PN order.

Binding Energy
Let us start with the effective action of binary neutron stars in pure gravity [53] where Γ µ = Γ µ αβ g αβ . The first term is the Einstein-Hilbert action, while the second term fixes the harmonic gauge. The dynamics of the two-body system is described by the third term using the world line approximation. In principle, one could have more generic couplings between gravitons and world lines, which appear at high PN order. Such terms are omitted for the moment. Now we consider a massive scalar field φ with We assume a reflection symmetry of V (φ), as in the axion case, which eliminates couplings such as φ 3 and φh 2 00 . Similar to self-interactions of gravitons, self-interaction vertices such as φ 4 and higher powers only contribute at higher order in the PN expansion 1 . For these reasons, it is enough to consider V (φ) = 1 2 m 2 s φ 2 . For the charged neutron star solutions discussed in Section 2, we should consider all possible couplings between the scalar and the world lines that respect the symmetry of the full theory, and therefore the last term in Eq. (3.3) becomes where q n and p n are the scalar couplings to the neutron star to be determined by utilizing matching conditions. Both p i and q i have mass dimension one. Here we only show the terms that contribute up to 1PN. Note that we also do not include u µ ∂ µ φ (where u µ is the 4velocity), which is proportional to the equation of motion (up to a total derivative) at leading order, and therefore is a redundant operator.
To calculate the binding energy as well as radiation in GR, we first expand the metric around Minkowski space then interactions between the point-like particles and the gravitons as well as the scalars are obtained by Taylor expanding action (3.5) in v. For example, We also have couplings between the scalar field and gravitons from Eq. (3.4), where η ij = −δ ij and the dot product between momenta k and q is defined as k · q = δ ij k i q j . Further more, we decompose h µν = H µν + h µν as well as φ = Φ +φ in the way that H µν (Φ) represents the off-shell potential graviton (scalar), whileh µν (φ) is the long-wavelength radiation graviton (scalar). The graviton propagator, which stems from the expansion of the Einstein-Hilbert action with gauge fixing conditions, is given by: where P µν,αβ = 1 2 (η µα η νβ + η να η µβ − η muν η αβ ). Given P 00,00 = 1/2 and P 00, Using the power counting rules in [53], we can find the Feynman diagrams, as shown in Fig. 2, Fig. 3, Fig. 4, and Fig. 5, that contribute to the binding energy up to 1PN. Among these diagrams, Fig. 2b, Fig. 3b, Fig. 3c, Fig. 4c and Fig. 5c represent the binding energy from the GR sector [53]. Together with the kinetic term, they give the Lagrangian for the binary neutron stars in pure gravity: with being the Einstein-Infeld-Hoffmann Lagrangian [55]. Corrections from the scalar field are represented by Fig. 2a, Fig. 3a, Fig. 4a, Fig. 4b, which contributes a Yukawa potential. At 1PN, the corrections are given by The first term in Eq. (3.17), 5 ms r , comes from the renormalization of the mass of the neutron star from axion mediated interactions at one loop, and can therefore be absorbed by redefining the mass of the neutron star. In the following, we will neglect this term in the axion potential since it is not observable. Collecting all the terms gives us the effective Lagrangian from the scalar sector up to 1PN order: For simplicity, we define the following dimensionless parameters where we have −1 ≤ α ≤ 1. Note that q > 0 if the scalar force between two neutron stars is attractive, and vice versa if repulsive. We also define: The 1PN correction to the Newtonian potential is given by and  For circular orbits with orbital frequencyΩ, the modified Kepler's law reads Before we move on, let us briefly comment on the magnitude of effects if the scalar field φ is the axion a in [30]. The leading order correction to the potential from an axion field can be matched to the axion potential between neutron stars calculated with the image charge method in [30]. This contribution, as we will show, is at order q R NS GM v 2 and the leading correction that distinguish an axion mediated force from a gravitational force. In Tab. 1, we summarize the dimension and the magnitude of the parameters as they are in the axion neutron star model studied in [30].
The leading corrections from the axion are given by Fig. 2a and Fig. 4a, and therefore These constitute a minimal model for the effects of the axion on the binary.

Radiation Power
Let us now study the corrections from the scalar field to the gravitational radiation at 1PN order. Our goal is to get the corrected radiation power, which is a necessary ingredient for calculating the inspiral waveform. The EFT for the radiation gravitons can be obtained by integrating out the potential H µν and "potential axion" Φ defined by φ = Φ +φ: where S full = S GR + S φ + S pp .

Gravitational wave radiation
Formally, the source term of the effective action of radiation gravitons can be written as where T µν is the pseudo-energy-stress tensor that can be read off from the path integral (3.30).
To manifest the PN order, it is not enough to just keep T µν at the right PN order; one should also expandh µν to the right PN order, which is achieved by performing multipole expansions around the center of mass [53]. Multipole expansion of actions is discussed in detail in [56]. Schematically, one can divide S source eff into two parts: the conserved part and the radiation one. The former hash 00 ,h 0i and their spatial derivatives coupled with conserved quantities, such as the ADM mass and momentum, and therefore does not radiate. The latter one has the form S rad eff = dt where R µνσρ is the linearized Riemann tensor defined by the metricḡ µν = η µν +h µν /m Pl and ijk is the Levi-Civita symbol. The I ij g , I ijk g , and J ij g are the mass quadrupole, mass octupole and current quadrupole, respectively, which, after extensive use of the Ward identity, doing integration by parts, and using the wave equation, are related to the pseudo-energy-stress tensor though where the dots denote time derivatives, and [· · · ] STF denotes the symmetric trace free components. Note that in the above equations, we have omitted terms that contribute at order higher than 1PN. (We refer the reader to [56] for more complete and compact expressions.) Finally, the power of gravitational radiation can be calculated using the optical theorem [53], Now, let us get back to the path integral (3.30) and find the expression for T µν . We only need to calculate T µν to finite PN order. According to S pp , we have T 00 ∼ vT 0i ∼ v 2 T ij at leading order. On the other hand, we have r∂ i R µνρσ ∼ vR µνρσ , since radiation gravitons carry a typical momentum k ∼ v/r. With these power counting rules, we conclude that I ij g ∼ vI ijk g ∼ vJ ij g at leading order. Thus, at leading order the gravitational radiation is simply with T 00 = n=1,2 M n . We find that the axion field has no effect on the gravitational radiation at leading order. Substituting I ij g into Eq. (3.36), one gets the well-known quadrupole formula where the brackets denote a time average. Calculation of the gravitational radiation power to next-to-leading-order needs the leading mass octupole, the leading current quadrupole, and the mass quadrupole up to O v 2 . According to Eq. (3.33), we only need to calculate T 00 , T kk and T 0i in I ij g up to O v 2 . The leading corrections from the scalar field are shown in Fig. 6, from which we find that all corrections have a magnitude of qv 2 . Given the smallness of q, they can be simply neglected at 1PN and therefore the gravitational wave radiation power is the same as in the case of pure gravity: where r is related to Ω through the modified Kepler's law (3.25) at 1PN. Note that we do not expand Eq. (3.38) in v 2 at this point. Φh ij (c) Figure 6: Corrections from the axion field on T 00 and T kk at 1PN.

Scalar radiation
In addition to gravitational radiation, there could be scalar radiation in the presence of the scalar field. Actually, a scalar field with a Compton wavelength much larger than the binary star separation leads to scalar radiation that dominates the energy loss, and therefore it should be severely constrained by observations. Similarly to gravitational radiation, the source term for scalar radiation can be written as where J is calculated in a PN expansion. In principle, to get the scalar radiation power at 1PN, we need to calculate J to 2PN order. This is because the power of dipole radiation is usually one PN order lower than that of quadrupole radiation. However, given the smallness of the scalar charge, we only need to calculate J to 1PN for dipole scalar radiation and at leading order for quadrupole scalar radiation. Diagrams that contribute to J up to 1PN are shown in Fig. 7, where In some cases, one may also want to include terms proportional to pv 4 , which come from the diagram in Fig. 8 and contribute  Collecting all the diagrams shown in Fig. 7, we have J = − nq n /m Pl with The radiation power can be calculated using where the multipole moments I L , which arise from multipole expanding the source action (3.40), are given by [56] Here L denotes a collection of index i 1 i 2 ...i l , and x L = x i 1 x i 2 ...x i l . For l = 0, we have x 2 a ∝ r, therefore dI/dt ∝ṙ which vanishes at 1PN. Thus, there is no monopole scalar radiation for circular orbits. For l = 1 we obtain, Similar to the case of axion potential, in the following, we keep only the leading contribution to the dipole axion radiation.

Matching to Axions
Consider the axion model of [30], which can lead to scalar charged neutron stars if To make use of the results above, we have to fix the parameters q and p, defined in Eq. (3.21), by matching with the full theory. This can be done in the limit of m a → 0, since the finite size effects of the neutron stars and therefore q and p should not depend on m a , as long as condition (3.51) is satisfied. In this case, the pure scalar potential between two static sources according to the EFT is According to [30], the charged neutron stars have constant axion field value at the edge of the stars. In this case, the scalar potential between two charged neutron stars of radius R NS can be calculated using the image charge method; it reads Note that despite of the simple relation, q i and Q i are different since the former is the Wilson coefficient we introduced in the EFT as well as the free parameter in the waveform, while the latter is the scalar charge of the neutron star in the specific axion model. The parameter p is therefore bounded from below by 1/8 in the limit R (i) NS = 2GM i and resides in the range (0.25, 0.4) for neutron stars that are consistent with various constraints (see reference [57] and reference within). Such a requirement ensures that corrections to the axion potential and radiation, enhanced by 16p and 8p compared to the corresponding GR corrections for potential and radiation, respectively, are the leading corrections that help distinguish axion mediated interactions from gravity.

Waveform and prospects for detection with Advanced LIGO
In this section, we will first calculate NS-NS and NS-BH merger waveforms with axion induced corrections based on the axion mediated force and axion radiation found in the previous section. We will then compare this GW waveform to the one within General Relativity, and assess the detectability of these corrections. The result of such a comparison is presented as a projected constraint on the axion parameter space. The method used in this section can be adapted to the study of any other theory where a light massive scalar is coupled to the neutron star or similar compact objects with an almost gravitational strength. In this section, we will only keep the leading corrections to the potential and radiation to an order that is relevant for breaking the degeneracy between axion induced corrections to the gravitational waveform and post-Newtonian corrections.

Waveform
The inspiral waveform measured in a gravitational wave detector is of the form of [58,59] with D L being the luminosity distance to the source and where Q encodes the detector response as a function of the angular position and orientation of the binary. For convenience, we neglect the cosmological red-shifting of the observed frequency of gravitational radiation (motivated by the limited horizon for neutron star binary mergers with current interferometers). In addition, we assume an ideally oriented binary and set Q = 1. If d ln h 0 /dt dφ/dt and d 2 φ/dt 2 (dφ/dt) 2 , the Fourier transform of the time-domain waveform can be computed using the stationary phase approximation, where In the above two equations, t should be thought as a function of f and defined as the time at which dφ/dt = 2πf . Usually one can solve for r(Ω) from the modified Kepler's law (3.25), and then get the analytical frequency domain waveform at 1PN. However, in the presence of a massive scalar field, Ω 2 is not analytical in terms of the PN parameters, and therefore we cannot solve r(Ω) in general. For this reason, we first calculate H and φ in terms of r, and translate them to f using a numerical interpolation function r(Ω) when we generate the waveform. The system we solve is given by where E = 1 2 M ηr 2 Ω 2 + V GR + V a , P = P g + P a and we have used Ω = πf . Together with Eqs. The potential and radiation terms we calculated in sections 3.1 and 3.2 contain all corrections at leading order in axion charge, and up to next to leading order in the PN expansion. These corrections will all be needed if we were to extract information about axions from LIGO data. The waveform we calculated numerically, and subsequently use to estimate the reach for the axions using a Markov Chain Monte Carlo (MCMC) sampler, however, only take corrections Eq. (3.23), Eq. (3.28) for the potential, and Eq. (3.49), Eq. (3.50) for the radiation into consideration. The potential terms include the Newtonian and 1PN corrections to the gravitational potential, the leading order axion potential and its leading correction due to "image charges". The radiation terms include the leading and next leading order gravitational wave quadrupole, as well as the axion dipole and quadrupole radiation and its leading correction from axion radiation from induced dipole. These terms are sufficient to break the degeneracy between axion induced corrections to the gravitational waveform and post-Newtonian corrections 2 .

Forecast
Given a high signal-to-noise (SNR) detection of a neutron star merger, it is possible to use the measured inspiral waveform not only to infer the parameter of the binary, but also to derive constraints on parameters in the scalar sector, such as q 1,2 , p 1,2 , and λ. A measured signal s(t,θ) consists of a noise realization n(t) and a merger waveformh(t,θ) depending on the "true" parametersθ, namely s(t,θ) = n(t) +h(t,θ). For a set of template waveforms g(t, θ), which depend on a set of candidate parameters θ, the likelihood function is where N is a normalization factor [60]. Given two signals h(t) and g(t), the inner product ( h | g ) on the vector space of signals is defined as where S n (f ) is the detector noise spectral density andh,g are the Fourier transforms of h, g. The inner product is defined so that the probability of having a noise realization n 0 (t) is p(n = n 0 ) ∝ exp[−(n 0 |n 0 )/2]. To find the average ∆χ 2 , one then marginalizes the logarithm of the likelihood over many noise realizations (e.g., [61]) ∆χ 2 (θ) ≡ 2 log L(s|θ)/L(s|θ) (4.11) where L(s|θ) is the likelihood evaluated at g =h with H and Ψ the amplitude and phase of the waveform in the stationary phase approximation. Assuming a Gaussian likelihood, one can interpret ∆χ as the number of "sigmas" at which the parameter set can be constrained given the noise model.

Forecasted constraints on q i and λ
To give an idea of the constraints on the parameters in the axion sector, we consider two fiducial scenarios. In the first scenario, we assume a neutron star binary with masses M 1 = 1.2M and M 2 = 1.24M , evolving at a luminosity distance of D L = 40 Mpc, in pure GR. We also assume the radii of the two neutron stars to be R where A ≡

24
G 5/6 M 5/6 π 2/3 D L and M ≡ µ 3/5 M 2/5 denotes the chirp mass. The parameter t c stands for the time at which the separation goes to zero in the Newtonian limit of the coalescing particles and φ c is the corresponding phase. Given our assumption of no axion field, the "true" values of the parameters in the scalar sector are q 1,2 = 0, p 1,2 = 10/16 and λ = ∞.
We perform a MCMC sampling of the likelihood function using the emcee package [62] on the full 10-dimensional parameter space. We use the forecasted noise curve for Advanced LIGO at design sensitivity ("Design") based on the Zero Det, High Power scenario [63]. This provides an idea of the noise-limited constraints that could be obtained by Advanced LIGO for a nearby NS-NS inspiral event. In the left panel of Fig. 9, we show the marginalized 3σ forecasted constraints in the (q − λ) plane where q = q 1 q 2 /M 1 M 2 (see, Eq. (3.21)). As it can be seen in the plot, there is a degeneracy between q and λ. Constraints on q become tight as λ increases. The GR limit can be achieved as q → 0 or λ → 0. Hence, it is not a single point in the parameter space. In principle, the 3σ constraint contours should approach a non-zero constant as λ goes to infinity. However, sampling this infinite ridge in the likelihood function in the large λ limit requires a prohibitively large number of samples. Therefore, Fig. 9 only shows the forecasted constraints in the small λ limit.
One may expect that the contours approach a constant in q as λ goes to infinity (massless axion limit), allowing us to fix this asymptotic constraint by sampling the likelihood function for λ → ∞. We obtain the following 3σ constraint on q 3 : A notable feature about the neutron star solutions discussed in this paper is the induced ("image") charge effects on the axion profile. Given that 8p > 1 (8p = 1 corresponds to the compaction of black holes, see Eq. (3.54)), q can be tightly constrained due the induced charge effects described by the last term in Eq. (3.28), especially in the large λ limit where the exponential suppression associated with the Yukawa potential is less important.
As a second scenario, we consider a binary system that consists of a 1.2M neutron star and a 1.24M black hole at D L = 40 Mpc in pure GR. For NS-BH binaries, the only effect the axion has on the inspiral dynamics is through scalar radiation, which can be characterized by q 1 and λ. Thus, the waveform template is in this case is parameterized by (4.13) We sample the likelihood function using the same method and noise curve as in the first scenario to derive noise-limited constraints that could be obtained by Advanced LIGO for 3 The 3σ constraint on negative q directly obtained from the MCMC sampling is −7.2 × 10 −9 < q, which is tighter than that on positive q. This asymmetry, also visible in the left panel of Fig. 9, is due to the way we cut the inspiral waveform. In the MCMC sampling, we cut the inspiral waveform at fixed separation. A negative q in the binding energy will delay the phase of waveform, leading to more cycles compared to a positive q. Therefore, the constraint on negative q is tighter than on positive q. Because the exact constraint is dependent on the cut we choose, we report |q| < 6.1 × 10 −8 for simplicity. A more complete model of the merger itself is necessary to evaluate the asymmetry between constraints on positive and negative charges.
a nearby NS-BH merger event. The marginalized 3σ constraints in the (q 1 /M 1 − λ) plane are shown in the right panel of Fig. 9. Because scalar radiation can be emitted only if the scalar wave frequency is larger than the mass of scalar field, the constraints on q 1 become weaker when λ becomes less than the typical wavelength corresponding to 10Hz, i.e., the lower bound of the LIGO observational band. Analogously to the first case, we perform the MCMC sampling in the limit of λ → ∞ to resolve the constrains on large λ. We find the following 3σ constraints on q from a NS-BH inspiral event: Let us briefly compare the constraints from the NS-NS and NS-BH mergers in the λ → ∞ limit considered above. The axion influences the NS-NS merger through an attractive or repulsive scalar force and scalar radiation (note however that for an attractive force, the charge dipole for the example we have chosen is small, and therefore the scalar radiation is negligible), while it influences the NS-BH merger only through the presence of scalar radiation. For the roughly equal mass binaries that we have considered, we can directly compare the constraint on q for the NS-NS event to |q 1 /M 1 | 2 < 3.2 × 10 −7 from the NS-BH event. It can be seen that a stronger constraint can be obtained from the NS-NS event, implying that the scalar force is driving the constraints more than the contribution from scalar radiation.

Forecasted constraints on the axion parameter space
Using Eq. (2.7) and Eq. (3.54), we map q i to f a , and thus project the constraints in Fig. 9 to the axion parameter space. The result for the two fiducial binary systems we have studied above are shown in Fig. 10. As above, the constraints on q i are sampled in two regimes, the small λ regime and the large λ regime. The constraints on the regime in between are found by simple interpolation and are plotted in dashed line. The NS-NS binary is more constraining than the NS-BH binary due to the stronger constraint on the scalar charge given in the NS-NS case.
Interestingly, our forecasted constraints for Advanced LIGO (the blue shaded region in the plot) are complementary to existing constraints on the axion parameter space, e.g., constraints from direct measurements of the Sun, from measurements of the orbital decay of binary pulsar systems, or of black hole super-radiance (see, Ref. [30] for a complete description of these constraints). From our analysis, we find that Advanced LIGO has the potential to pin down the axion mass and decay constant within the range m a 10 −11 eV, f a (10 −14 − 10 −17 ) GeV, (4.14) or, in the absence of a detection, to exclude axions with m a and f a in this region of parameter space.
Axions with masses and decay constants in the above range are most likely in significant tension with the detected gravitational wave signal from the binary neutron star event GW179817 [4]. After having developed the methods and tools to constrain the axion parameter space from a given waveform in this work, our aim is to proceed with an analysis of data  Figure 10: Forecasted marginalized constraints on the axion parameter space. Colored regions below the curves can be constrained by analyzing data of the LIGO detections of NS-NS (blue) and NS-BH (red) mergers. The forecasted constraints are found by performing MCMC sampling in the small m a limit (the horizontal solid lines) and in the large m a limit (the curly solid lines). We connect these forecasted constraints by linear interpolation (the dashed lines).
For comparison, we also show the existing constraints (in gray) from direct measurements of the Sun, from measurements of the orbital decay of binary pulsar systems and from black hole super-radiance [30]. The region above the dotted purple line are parameter spaces where an axion profile can be sourced by a neutron star. The black line shows the parameters of a QCD axion, while the dotted gray horizontal line marks the value of the reduced planck scale M pl .
from existing and future gravitational wave detections next. In particular, we would like to apply our analysis pipeline to the GW179817 [4] data in a follow-up work.

Conclusions and discussion
In this paper, we have examined the exciting possibility for Advanced LIGO to detect new light scalar particles through their influence on the gravitational waveform produced in NS-NS and NS-BH binary mergers. Employing an EFT approach, we have calculated the first relativistic corrections to the binary orbital dynamics and gravitational waveform in the presence of a light scalar coupled to neutron stars. We use this waveform to forecast the constraints from Advanced LIGO for an event similar to GW179817 on the parameters in the EFT, which in the scalar sector include the charges of the neutron stars and a relativistic correction corresponding to image-charge effects. This result, summarized in Fig. 9, is applicable to theories where a light scalar couples to neutron stars with near gravitational strength. We then specialize to a particularly well-motivated light scalar, the axion. If there are in fact axion(s) with mass(es) and decay constant(s) in a region of parameter space where Advanced LIGO has a good sensitivity, then as it can be seen from Fig. 9, the parameters of the EFT can be measured with high precision. Such a scenario is most likely under significant pressure from GW179817 [4], however one can speculate about the implications should a detection be imminent.
The discovery of a new particle using an entirely new observable would, of course, be an incredible development in its own right. One of the first follow-up questions would be: what other physical phenomena could this new particle be related to? Axions with nuclear couplings and masses in this range can potentially solve the strong CP problem of the Standard Model [38]. In addition, the detected axion could in principle be a dark matter candidate. In the region of parameter space accessible to binary NS mergers, the axion must be produced non-thermally, implying evidence for a non-trivial cosmological history.
There are a number of avenues for finding corroborating evidence to a detection of axions with LIGO. If the axion is the dark matter, the same nuclear coupling that leads to a force between neutron stars also leads to a time-dependent nuclear electric dipole moment that can be targeted by precision magnetometry [47,64] as well as various resonant experiments that look for photon couplings of the axion. Precise knowledge of where to look in parameter space can greatly improve the prospects for detectability using such techniques. A precise knowledge of the masses (and couplings) of the axion will significantly narrow down the range of axion masses to scan, while the sensitivity to the axion coupling improves as (m a /δm a ) 1/4 since more time can be allocated into the frequency range where the mass lands (sensitivity scales as t 1/4 ). The region of axion parameter space covered by binary NS mergers is also accessible to probes of black hole superradiance, e.g, gaps in the distribution of black hole spin or gravitational waves from rotating axion clouds [21,40].
The axion would also provide an interesting additional probe of the structure of the merging neutron stars. From Eq. (2.4), the scalar charge of the individual neutron stars is dependent upon the compaction (recall that the compaction is defined as GM/R) and the decay constant f a . Also note that from Eq. (3.54), the EFT parameter p is sensitive to the compaction of the neutron stars. For an event with sufficiently high SNR, the axion force therefore provides a new way to constrain the compaction of neutron stars. With future detectors, it may also be possible to use the post-merger waveform associated with the hypermassive neutron star resulting from the merger event to provide further knowledge on both the properties of the axion and nuclear equation of state [31]. We leave further investigation of the post-merger signal for future work.
In the absence of a detection, it is possible to set stringent constrains on the region of parameter shown in Fig. 10, for axions that possess a nuclear coupling. However, it is not necessary for all such particles to possess a nuclear coupling. Nevertheless, the lack of a detection would imply that laboratory experiments relying on such couplings, such as CASPEr-Electric [47], should also fail to make a detection over the same region of parameter space. Knowing where not to look could be useful in guiding such searches. The effects associated with superradiance could in principle be found even in the lack of a detection from binary neutron stars. In this case, one would strongly constrain the nuclear coupling of the axion, and potentially the QCD axion (the only target for laboratory experiments looking for an axion through its nuclear couplings).
Let us also briefly comment on the possibilities with future gravitational wave detectors. For the purpose of constraining light scalars from neutron star mergers, Advanced LIGO is limited by its overall sensitivity and its frequency coverage. An increase in sensitivity over the Advanced LIGO band, as would be provided by third generation gravitational wave detectors such as Einstein Telescope [65], would yield a number of advances. The projected constraints on scalar charge would become tighter as the SNR per event would be higher. Greater detection rates would allow for a joint analysis (e.g. "stacking" the events) of many events that could provide stronger projected constraints than individual events. In addition, greater sensitivity at high frequencies could provide access to the end-stages of the inspiral and the ring-down of the hyper-massive neutron star or black hole that can form as a result of the merger. This would provide new information about the scalar sector through additional relativistic corrections, and through effects on the structure and evolution of post-merger objects (for example, as explored in Ref. [31]). A space mission such as LISA [66] will provide sensitivity at lower frequencies. For individual events, this would provide access to scalars with a lower mass as the binary evolution could be tracked at larger separation. In addition, the projected reach on the charge dipole of the binary would improve since orbital energy-loss due to scalar radiation is more important at lower frequencies. Finally, it will be possible to observe the merger of white dwarfs (either individually or as a stochastic background), which would allow one to examine the nature of the coupling between axions and compact objects. In particular, it would be interesting to examine the density-dependent coupling invoked in the axion model we have studied here.
In summary, the observation of binary neutron star mergers provides a novel opportunity to search for new light scalar particles, including axions. The waveforms presented in this paper, and the forecasted constraints, provide the technical basis and proof-of-concept necessary to proceed with an analysis of data from existing and future events. In particular, we hope to perform an analysis using data from the existing event GW179817 in future work. The results of such an analysis will greatly inform other observational and laboratory efforts to search for light scalars, and provide constraints over an extensive and well-motivated region of parameter space for axions.

A Degeneracy with higher PN corrections
The waveform obtained above considers only 1PN corrections to GR. In principle, one can improve the waveform by simply replacing the 1PN expressions of the gravity sector with higher PN expressions. In the limit we consider in this paper, where spin and tidal effects are neglected, in this section we estimate the impact of including higher PN corrections in the gravity sector on our constraints. The constraints on the axion can be characterized by the phase difference caused by the axion field. Specifically, we consider a binary system composed by 1.2M and 1.24M masses and assume the two stars carry the same scalar charge. We calculate the total phases, Ψ GR and Ψ A , by integrating the phase over (10 -1000) Hz in the cases with and without the axion. The constraints on the axion field can be characterized by the differences of total phases, ∆Ψ = |Ψ A − Ψ GR |. From Eq. 4.11, it can be seen that ∆χ 2 can be significant only once the phase difference is order one. We calculate ∆Ψ using different PN order expressions in gravity sector, and the plot the contours of ∆Ψ = 1 in Fig. 11. We find that the region of parameter space over which the phase difference is order one, as shown in Fig. 11, does not change significantly when including higher PN terms in the gravity sector, especially for the parameter range we are interested in. It is therefore justified to just use the 1PN correction in order to forecast constraints on the axion parameters.