Cuckoo's Eggs in Neutron Stars: Can LIGO Hear Chirps from the Dark Sector?

We explore in detail the possibility that gravitational wave signals from binary inspirals are affected by a new force that couples only to dark matter particles. We discuss the impact of both the new force acting between the binary partners as well as radiation of the force carrier. We identify numerous constraints on any such scenario, ultimately concluding that observable effects on the dynamics of binary inspirals due to such a force are not possible if the dark matter is accrued during ordinary stellar evolution. Constraints arise from the requirement that the astronomical body be able to collect and bind at small enough radius an adequate number of dark matter particles, from the requirement that the particles thus collected remain bound to neutron stars in the presence of another neutron star, and from the requirement that the theory allows old neutron stars to exist and retain their charge. Thus, we show that any deviation from the predictions of general relativity observed in binary inspirals must be due either to the material properties of the inspiraling objects themselves, such as a tidal deformability, to a true fifth force coupled to baryons, or to a non-standard production mechanism for the dark matter cores of neutron stars. Viable scenarios of the latter type include production of dark matter in exotic neutron decays, or the formation of compact dark matter objects in the early Universe that later seed star formation or are captured by stars.

We explore in detail the possibility that gravitational wave signals from binary inspirals are affected by a new force that couples only to dark matter particles. We discuss the impact of both the new force acting between the binary partners as well as radiation of the force carrier. We identify numerous constraints on any such scenario, ultimately concluding that observable effects on the dynamics of binary inspirals due to such a force are not possible if the dark matter is accrued during ordinary stellar evolution. Constraints arise from the requirement that the astronomical body be able to collect and bind at small enough radius an adequate number of dark matter particles, from the requirement that the particles thus collected remain bound to neutron stars in the presence of another neutron star, and from the requirement that the theory allows old neutron stars to exist and retain their charge. Thus, we show that any deviation from the predictions of general relativity observed in binary inspirals must be due either to the material properties of the inspiraling objects themselves, such as a tidal deformability, to a true fifth force coupled to baryons, or to a nonstandard production mechanism for the dark matter cores of neutron stars. Viable scenarios of the latter type include production of dark matter in exotic neutron decays, or the formation of compact dark matter objects in the early Universe that later seed star formation or are captured by stars.

I. INTRODUCTION
The dawn of gravitational wave (GW) astronomy, commencing with the measurement of GWs produced by coalescing binary black holes (BHs), promises an entirely new era of observational astrophysics [1][2][3]. Recently, coalescing neutron star binaries have also been observed for the first time [4]. Such events are treasure troves of new data for astrophysics and astroparticle physics. They not only feature an associated electromagnetic signal [5,6], but also produce gravitational waveforms that are potentially observable over time periods of order 20 minutes and at length scales which are different from observations of binary pulsars [7].
The ability to track the neutron star binary inspiral over an extended period of time provides increased sensitivity to modifications of the post-Newtonian phase of the inspiral compared to inspirals of O(30M ) black holes. These deviations can be broadly categorized as stemming either from (i) a modification of gravity [8][9][10][11][12][13][14][15][16][17], (ii) a "fifth force" mediated by particles such as axions [18] or ultra-light dark gauge bosons [19], or from (iii) a new long-range force coupled to dark matter (DM) particles accreted inside neutron stars [20][21][22]. Stringent constraints exist on a new fifth force coupled to baryons, see for instance Refs. [23][24][25][26][27] and references therein, which have motivated a particular focus on the possibility of a new long-range force, mediated by ultra-light bosons, that couples only to dark matter or other cosmologia jkopp@uni-mainz.de b ranjalah@uni-mainz.de c opferkuch@uni-mainz.de d shepherd@uni-mainz.de cally stable hidden sector particles residing inside neutron stars. Such scenarios have been studied extensively using probes other than gravitational waves, in particular the properties (and the mere existence) of old neutron stars and white dwarfs  as well as the dynamics of binary pulsar PSR 1913+16 (Hulse-Taylor binary) [53,54]. We note that dark sector particles experiencing a strong long-range force can at most contribute a small fraction to the total DM density of the Universe due to tight constraints on DM self-interactions [55][56][57].
Complementing these results, we will show that DM accumulation in neutron stars throughout the course of ordinary stellar evolution can never be efficient enough for the new force to lead to observable deviations in gravitational wave signals from binary inspirals. Detectable, per cent-level deviations require 10% of the neutron star mass to be made up of dark sector particles. Exploring the constraints on the amount of such particles that a neutron star can contain, we will conclude that only exotic neutron decays or the formation of dark stars can result in such large abundances.
We begin in Section II by reviewing fifth-force effects in neutron star binary inspirals, focusing on the required region of parameter space to observe deviations from general relativity in LIGO/VIRGO data. Of course, our results are easily generalized to other gravitational wave detectors such as ALIA, BBO, DECIGO, the Einstein Telescope, Geo-600, KAGRA, LISA, TAMA, etc. [58,59]. We then derive constraints on these forces, distinguishing between repulsive forces (Section III) and attractive ones (Section IV). We conclude in Section V. Note that we use natural units throughout the text.

II. INSPIRALS AND FIFTH-FORCES
Before turning to the constraints on the interactions of any exotic hidden sector particle, we first illustrate the necessary conditions for an observable new physics signal at LIGO/VIRGO. In doing so, we closely follow Refs. [10,19]. To remain as model-independent as possible, we parameterize the effect of the new dark force using a generic Yukawa potential. We furthermore assume that the new force couples only to dark sector particles. In order for it to act over sufficiently long distances of order 100 km, the mediator must be ultra-light ( 10 −12 eV), and the dark sector particles must neither significantly screen the new force nor efficiently self-annihilate. The most straightforward way of realizing the second constraint is asymmetric DM [60,61].
Let us consider a binary neutron star system, and let us assume that each of the binary partners contains a population of dark sector particles. These populations will affect the inspiral dynamics in two distinct ways: first, the exotic force acting between them will affect the time evolution of the distance between the neutron stars, their orbital frequency, and the time of the merger. Second, radiation of the new light force carriers provides an extra energy loss mechanism.

A. Effects of a new Yukawa force
Even in the presence of a new dark sector Yukawa force, the system's gravitational wave emission will follow the predictions of general relativity as long as the distance between the binary partners is significantly larger than the range of that force. However, once their distance drops below that range and the exponential suppression of the Yukawa potential is lifted, observable deviations may occur. The assumption of a Yukawa force is necessary, as the potential effects of a new infinite range force are degenerate with a shift of the neutron star masses. The new force, which can be either attractive or repulsive, results in a modification of the chirp mass is the reduced mass of the inspiraling objects with individual masses M 1 and M 2 . Assuming point-like neutron stars, we can write the magnitude of the force between the two bodies as where ∆ is their spatial separation and m med is the mass of the particle sourcing the Yukawa potential. Here α parameterizes the size of the new force relative to gravity in the regime m med ∆ 1: The sign in this expression determines whether the force is attractive (positive sign) or repulsive (negative sign).
α is the analog of the fine structure constant for the new force, and Q i = q χ N χi are the total charges of the inspiraling objects under the new force, which depend on the charge q χ of a single DM particle and the numbers N χi of captured DM particles. Given that the orbits of the binary system will have circularized by the time the gravitational wave signal becomes observable [62], the orbital frequency of the system is given simply by a modified Kepler's law, 1 The frequency of the gravitational wave signal is f GW = ω/π; therefore if LIGO sensitivity begins at O(10 Hz) [4], Eq. (3) tells us that a binary system consisting of two 1.25 M neutron stars enters the sensitivity band at a spatial separation of O(700 km) (assuming that gravity is the dominant force at these length scales, i.e. | α | 1). A typical neutron star radius being 10 km indicates that a Yukawa force detectable by LIGO must have range m −1 med O(20 − 750 km), i.e. m med 10 −11 -3 × 10 −13 eV.
To determine the inspiral dynamics we require the total energy of the system, The last term is the kinetic energy of the stars, neglecting their spin and internal structure. The power radiated via gravitational waves is [63,64] Equating − dE GW / dt to the time derivative of Eq. (4), we obtain the rate of change of the orbital frequency, where In the massless mediator limit, this reduces to g = −3, and the classical gravity-only result is recovered [63]. In this limit the effect of the infinite range fifth-force results merely in a modification of the apparent strength of gravity, parameterized by the replacement G N → G N (1 + α ) GW (t) as a function of time for a neutron star binary with various values of the fifth-force parameters α and m med . In particular, we show results for m −1 med = {100, 300}km (solid, dashed lines) and three choices of α = {1, 0.25, 0.025} (red, purple, and blue colored lines). We consider here the limit that the system loses energy only by gravitational wave emission, but not by radiation of light mediator particles. This situation is realized approximately if the binary partners have equal dark chargeover-mass ratios. The gray band indicates the typical error on reconstructing the chirp mass, which we take to be 0.4% based on ref. [4]. The boundary condition of the differential equation is the threshold sensitivity of LIGO, fGW(0) = 10 Hz [65]. The neutron star masses are M1 = M2 = 1.25M .
in Eq. (3). Equation (3) can then be used to eliminate ∆ from Eq. (6), yielding We see that an infinite range Coulomb-like force results in a modified chirp mass The system is thus indistinguishable from a purely gravitationally interacting system with a different chirp mass. GW rather than just f GW because this choice yields straight lines in the gravity only scenario (dotted black line in Fig. 1). For m −1 med = 300 km and large values of α , we see that dω/dt is a still constant, but clearly different to the gravity only solution, for the entire frequency range in which LIGO is sensitive. In other words, the Yukawa force looks like a Coulomb force to LIGO, and we would simply reconstruct the event with an alternative, incorrect, value of the chirp mass. There is thus no sensitivity to the new fifth-force for these parameters from GW signals alone; an electromagnetic signal associated with a GW event with reconstructed masses above the Chandrasekhar limit would indicate a modification induced by a dark force.
On the other hand, larger mediator masses and smaller fifth-force couplings lead to deviations that occur partway through the time domain of the LIGO sensitivity band. These scenarios would be observed as prematurely coalescing binaries. Even for α = 0.025 and m −1 med = 100 km, the coalescence time is shifted by ∼ 6 s relative to that predicted by gravity alone (with the chirp mass determined from the first part of the wave form).
The modification of Kepler's law, Eq. (3), changes the initial separation of the neutron stars compared to the gravity-only case at the time when the waveform enters the LIGO frequency band. This in turn can modify the amplitude of the signal. For attractive (repulsive) dark forces one obtains an increase (decrease) in the amplitude. More concretely, the gravitational wave amplitude, A GW (t), is given by [19,66] where d L is the luminosity distance of the GW source. In the pure gravity scenario, A GW × d L is 6.3 × 10 −22 Mpc at the time when the GW emission enters the LIGO frequency band. In the presence of a dark force with | α | = 0.25 and mediator mass m −1 med = 300 km, the amplitude at this time becomes 6.65 × 10 −22 Mpc and 5.96 × 10 −22 Mpc for an attractive and repulsive dark force, respectively.

B. Ultra-light Boson Radiation
Another source of energy loss is radiation of the ultralight mediator particles [10,19]. The leading contribution to the power radiated via vector mediators in the multipole expansion is [53,67] where γ ≡ Q 1 /M 1 − Q 2 /M 2 is the difference in charge to mass ratios of the binary partners. Equation (11) thus illustrates that dipole radiation is only possible if the two stars have different dark charge-over-mass ratio. For scalar mediators, the power of the dipole radiation is [53,67] In what follows, we will focus on vector radiation as a concrete example. The energy loss to both gravitational wave and ultralight boson radiation must equal the orbital energy loss In Fig. 2 we consider the limit where the system loses energy to dipole radiation and to gravitational waves, but gravitational wave emission is not affected by a dark force acting between the binary partners. In other words, we set α = 0, but β = 0 in Eq. (13). This scenario corresponds to the case of only one neutron star being appreciably charged. We obtain where we have defined β parameterizes the magnitude of the radiation effect relative to gravity in the same way as α from Eq. (2) parameterizes the effect of a new force between the binary partners. There is a direct correspondence α and β when the masses of the neutron stars are the same and we take Q 1 = Q 2 for α and Q 2 = 0 for β , i.e.
In Fig. 2 we show solutions of Eq. (14) varying both m med and β . The largest value of m −1 med shown, 12 000 km, is chosen such that dipole radiation is switched on while the waveform is still outside LIGO's frequency band. We see that the effect of radiation is observable even in this case, in contrast to the effect produced by a dark force with the same range, see Fig. 1. The reason is the different ways in which the right hand sides of Eqs. (8) and (14)  GW (t) as a function of time for a neutron star binary with various values of the fifth-force parameters β and m med . We consider here the limit that the system loses energy via dipole radiation and gravitational wave emission, but that the latter is unaffected by a dark force acting between the binary partners. This situation is realized if one of the stars does not contain any dark sector particles. We show results for m −1 med = {4000, 8000, 12 000 }km (dotted, dashed, solid lines) and three choices of β = {10 −1 , 10 −2 , 10 −3 } (red, purple, and blue colored lines). The gray band indicates the typical error on reconstructing the chirp mass, which we take to be 0.4% based on ref. [4]. The boundary condition of the differential equation is the threshold sensitivity of LIGO, fGW(0) = 10 Hz [65]. The neutron star masses are M1 = M2 = 1.25M .
This means that the sensitivity to β in the radiationdominated scenario is about an order of magnitude better than the sensitivity to α in the dark force-dominated scenario from Section II A. This difference arises because the two effects arise at different order in an expansion in ω∆ ∼ 0.1.
The upshot of this section is that for neutron star binaries in which both stars carry a dark charge and have similar dark charge-to-mass ratio, LIGO may be able to establish the existence of a new long-range force if  Of course, radiation of dark force into higher multipoles is possible even for systems that are symmetric in their charge-to-mass ratio. However, such radiation is suppressed by additional powers of ω∆, so we do not expect it to dominate over the effect of the new longrange force between the binary partners. Moreover, the effect of quadrupole radiation, like the effect of the dark force, would be degenerate with the effect of gravity at large m −1 med , making it unobservable in this regime independent of the suppression by ω∆. At smaller m −1 med , however, radiation into higher multipoles would lead to the same kinks as dipole radiation (see Fig. 2), offering an extra handle for establishing the radiation effect. The kinks occur within the LIGO sensitivity window for 4000 km < m −1 med < 12 000 km. In the following sections we will investigate the viability of the parameter values indicated above.
We note that the reach in m −1 med may be significantly enhanced in the future by combining LIGO/VIRGO measurements with results from a space-based experiment like LISA [68,69]. LISA would be able to observe neutron star binary systems long before the merger. If the separation of the binary partners is larger than m −1 med while the system is observable by LISA, but below m −1 med when it enters the LIGO/VIRGO frequency band, the effect of the dark force can be observed in a combined analysis.
To end this section we remark that while the above discussion has been in the context of neutron stars these results are also applicable to black hole or mixed binaries so long as the black holes can carry the appropriate charge.

III. CONSTRAINTS ON REPULSIVE FORCES
We now focus our attention on the case of a repulsive force between the dark particles captured by binary neutron stars. Our results on this scenario are summarized in Fig. 4. The four panels of this figure show, for different assumptions on the particle mass m χ , the limits on the dark fine structure constant α and the total mass of all χ particles in the neutron star, M χ . We normalize α to the effective strength of gravity, G N m 2 χ , and M χ to the mass of the baryons in the neutron star, M b . In the following, we discuss the various constraints shown in Fig. 4 one by one.
Throughout our discussion of constraints we shall quote bounds on the dark force parameter α under the assumption that both neutron stars have identical dark charges and masses. In the opposite limit that one neutron star carries a dark charge while the other does not, equal constraints apply on the parameter β . It is crucial here that we have defined β such that under the assumption that only one neutron star is charged its value is equal to that of α for equally charged neutron stars.

A. Binding potential for DM particles
It is intuitively clear that a repulsive force between dark particles will limit the maximum number of such particles that can accumulate in or around a neutron star or black hole. Irrespective of the capture or production mechanism, the number of particles is saturated once the net potential of the compact object ceases to be attractive. Consider the potential V (r) felt by a single DM particle χ of charge q χ and mass m χ at a distance r from the center of the star: Here, Q is the total dark charge of the star, and is its total mass (including the total baryonic mass, M b , and the total mass of dark particles, M χ ). M χ can be traded for the number of particles, N χ , of a given mass m χ : Insisting that the potential is attractive (V (r) ≤ 0) at length scales r m −1 med (where the Yukawa force is effectively Coulomb-like) yields the condition This constraint is shown in light blue in Fig. 4. We can use Eq. (20) to derive a constraint on the effective strength α of the new force relative to gravity. Plugging Eq. (20) into Eq. (2), we obtain B. Size of the Dark Core In order for the effect of the dark force on the neutron star to be predictable using the calculations of Section II,  it is necessary that the core of dark sector particles be significantly smaller than the screening radius of the force. Otherwise, the interactions of the neutron stars will be primarily with a cloud of dark particles in which they are both embedded, rather than remaining a simple centralforce problem. This would require a much more sophisticated treatment. We leave the dynamics of such a system for future work, and here aim to enforce the requirement that the dark particle core be smaller in extent than the screening radius m −1 med . We remark however that the far-outlying particles at the edge of the neutron stars' dark matter halo will be very easily stripped off in the presence of the second neutron star because they feel the dark force field created by that star's DM halo long before the inspiral dynamics can be affected. Thus, those outlying particles would not affect the gravitational wave signal, as they would no longer be bound to just one of the two stars.
Using the virial theorem, we can calculate the stable core radius for given α , m χ , and N χ . The general form of the virial theorem (valid for non-relativistic and relativistic systems) is , where i runs over all particles, and r i , p i , F i are their positions, their momenta, and the forces they experience, respectively. The notation . . . indicates time averaging, which here translates to the requirement that the system is in a stable configuration. For an isotropic potential, the virial equation in its relativistic form simplifies to [71,72] i p 2 The potential energy for a single particle at the edge of the dark core (radius R χ ) is given by The first term is the combination of the gravitational self-interactions of the exotic particles and the new dark interactions in the zero screening limit R χ m med 1. The second term is the gravitational interaction between the baryons and the test particle defined separately in the two cases where the χ-core is either smaller or larger than the radius R of the baryonic matter. We assume that the dominant contribution to the average momentum on the left hand side of Eq. (22) arises from Pauli repulsion. (We will justify this assumption shortly.) The p i are then of the order of the Fermi momentum of a degenerate 3dimensional system of free spin-1/2 fermions (with no additional internal degrees of freedom) where V is the volume of the χ-core. For repulsive forces, we find that the system is always in a non-relativistic configuration in the parameter regions of interest to us (unshaded portions of Fig. 4), allowing the simplification     2)) at each point in the plot. We indicate in red the estimated minimum value of α or β required for the dark force to leave an observable imprint in LIGO data. Regions of parameter space where the neutron star no longer has an overall attractive potential for DM, i.e., Nχ exceeds the bound of Eq. (20), are shaded in blue. Note that this bound is independent of the DM production or capture mechanism. Bounds shown in black are model-dependent: in the region above the black-dashed line, the DM population inside the neutron cannot be generated by capture from the halo, based on the geometric arguments that lead to Eq. (28). Above the black solid line, the hypothesis that the DM population is produced during the supernova that created the neutron star is not tenable because the energy available in a supernova (as measured in SN 1987A) is not sufficient. Inside the yellow shaded region, the dark core has a radius larger than the range of new Yukawa force, Rχ ≥ m −1 med /2. The solid, dot-dashed, and dashed yellow lines correspond to m −1 med = {100 km, 1000 km, 1 × 10 5 km}, respectively. Within the yellow regions, the DM particles would be stripped away from their host neutron star before they can have a significant impact on the inspiral. The neutron star parameters are chosen in accordance with Table I, and we normalize α such that qχ = 1.
Inserting Eq. (25) and Eq. (23) into the virial theorem, Eq. (22), yields The different contributions to this equation are represented graphically in Fig. 3, where dashed lines indicate repulsive forces and solid lines indicate attractive forces. The plots also show the contribution from ordinary thermal pressure, which, however, is relevant only for relatively heavy DM (m χ TeV), weak dark forces (α ∼ G N m 2 χ ), and small dark cores (M χ 10 −20 M b ). This parameter region is not of interest to us as it corresponds to tiny α ; this justifies our previous assumption that degeneracy pressure dominates. The intersection between the sum of all repulsive contributions and the sum of all attractive contributions corresponds to the solution of Eq. (26) and thus gives the radius of the dark core as a function of m χ , α , and N χ .
For light DM mass, Fermi-repulsion results in radii larger than the size of the neutron star. The regions of parameter space where R χ > m −1 med /2 are shaded in yellow in Fig. 4. A χ-core with such an extent results in inspiral dynamics very different to those assumed in our estimates so far. In particular, as argued above, particles at radii larger than m −1 med will be efficiently stripped from the neutron star. The yellow curves, therefore, can be taken as an indication of the maximum amount of charge that fits geometrically within the domain of influence of a single neutron star, and thus can contribute to a dark-force-induced change in inspiral dynamics and gravitational wave signals.

C. Constraints on Dark Core Production via Particle Accretion
There exists an additional upper bound on N χ under the assumption that the particles giving the neutron star its dark charge have been accreted from the host galaxy's DM halo. An upper limit on the number N capt χ of DM particles (χ) accreted can be obtained by assuming that any DM particle that passes through the neutron star over its lifetime t NS is captured. (A more realistic estimate based on the dynamics of baryon-χ scattering is discussed in Appendix A.) The upper limit on N capt χ is based on the effective geometric cross section of the neutron star, σ eff geom , given by where v esc is the escape velocity at the surface of the neutron star andv is the average χ velocity relative to the neutron star at infinity. The velocity-dependent term by which σ eff geom differs from the geometric cross section σ geom = πR 2 accounts for the focusing of the DM wind by the neutron star's gravitational field [73]. In the second equality in Eq. (27), we use v esc >v, an approximation that is well satisfied for typical galactic DM velocities v ∼ 200 km s −1 and escape velocities at the surface of a neutron star v esc ∼ 100 000 km s −1 . The upper limit on N capt χ is then where ρ χ is the mass density of χ in the vicinity of the neutron star. This constraint is shown as a horizontal black dashed line in Fig. 4. As the constraints on N χ from Eqs. (20) and (28) scale inversely to one another with respect to m χ , the value m opt χ that maximizes N χ and thus α occurs when both constraints are equal. We find and, by plugging m opt χ into Eq. (21), The second equality follows because the mass of the accreted DM particles is much smaller than the baryonic mass M b . For the numerical estimate in the second line of Eq. (30), we have used the neutron star parameters given in Table I, we have (very conservatively) assumed that ρ χ saturates the observed DM density in the solar neighborhood ρ χ = 0.3 GeV/cm 3 , and we have taken v = 220 km s −1 . We note that this limit on α is independent of both m χ and α . The limit is relaxed in regions of very low DM velocityv.

D. Constraints on Dark Core Production in Supernovae
The upper limit on N χ from Eq. (28), and the resulting limit on α , Eq. (30), were based on the assumption that the DM population bound to the neutron star arises from capture of DM particles from the halo. Relaxing this assumption and invoking more speculative production mechanisms can lead to significant boosts in the allowed values of α . We first note that the largest α values compatible with the constraints of Section III A (blue regions in Fig. 4) occur in the region of parameter space where the dark force is not much stronger than gravity and where M χ ∼ M b .
One possibility for producing a massive dark core without accreting DM from the halo arises if DM particles are produced via bremsstrahlung in a new-born hot neutron star [21,22]. The upper limit on the total mass of exotic particles produced this way arises from the measured neutrino burst of SN 1987A [74][75][76][77][78]. Adopting the 'Raffelt criterion' [79] that χ production should not reduce the duration of the neutrino burst from SN 1987A by more than 50%, we require that the instantaneous χ luminosity should be L χ ≤ 3 × 10 52 ergs s −1 . Various different effects need to be taken into account to evaluate L χ [80][81][82][83][84]. A near future detection of supernova neutrinos in all flavors will robustly strengthen these bounds [85][86][87].
We conservatively estimate that the total available energy to produce χ during the supernova is 10 53 ergs, which, under the further conservative assumption that these particles are produced at rest, translates to M χ 0.045M . This constraint is very weak due to the current scarcity of supernova neutrino measurements, and corresponds to approximately 1/3 of the total energy of the supernova being used to produce DM. This constraint is shown as a solid black horizontal line in Fig. 4. We can turn it into a constraint on α by requiring that it is saturated simultaneously with the the upper bound on N χ from Eq. (20) (blue shaded region in Fig. 4). The resulting equation can be solved for α , and the result can be plugged into the definition of α , Eq. (2). If we assume that the two neutron stars have similar mass and dark charge, we obtain the m χ -independent limit Note that the known thermal mechanisms to generate exotic particles in supernovae are inefficient for particle masses 10 MeV. In Fig. 4, we nevertheless plot the supernova constraint for m χ 10 MeV as an absolute upper limit on the amount of energy available to produce exotic particles in the supernova. We are, however, not aware of any mechanism that could possibly be efficient in converting an O (30%) fraction of the supernova energy into heavy particles. Any supernova-driven mechanism for DM production would lead to all neutron stars being similarly charged, and thus the LIGO sensitivity to these mechanisms is best modeled in Fig. 1, as the charge-tomass ratios of the binary stars will be similar.

E. Constraints on Dark Core Production via Neutron Decay
An alternative production mechanism for DM particles inside a neutron star has been proposed in Ref. [88] and further investigated in Ref. [89]. In this scenario, χ is abundantly produced inside the neutron star via a hypothetical exotic neutron decay mode. In this case a repulsive dark force would have an additional benefit: without it, the neutron star equation of state (EOS) would be softened because of decreased Fermi repulsion in the nuclear matter. This would preclude the existence of high-mass neutron stars with M 2M , in conflict with observations [90,91]. A repulsive dark force, however, stiffens the equation of state again. This places a lower bound on the size of the dark fine structure constant α . Ref. [89] found a bound on the strength of the needed repulsive force as a function of the mediator mass. In our notation, this constraint reads Note, however, that the calculations of Ref. [89] were performed in a scenario with a partially screened dark force, whereas scenarios in which the dark force could significantly impact inspiral dynamics and gravitational wave emission would require it to be unscreened inside the neutron star. This will weaken the constraint on α from EOS arguments compared to Eq. (32). Moreover, the neutron decay scenario requires careful tuning of m χ : stability of 9 Be places an upper bound on the mass splitting between χ and the neutron of 1.59 MeV [92,93], and the lifetime of 11 Be further strengthens this bound to require splittings of less than 0.50 MeV [94].
To estimate the maximum quantity of hidden sector particles that can be produced via neutron decay, we consider the chemical potentials of neutrons and DM. We describe the neutrons by a Fermi gas with chemical potential where n n is the number density of neutrons (which decreases over time in this scenario), and E B is a typical nuclear binding energy, taken here to be E B = 9 MeV [89]. The second term under the square root corresponds to the square of the neutron Fermi momentum. We assume the χ particles experience an unscreened potential, i.e., m −1 med R, resulting in the chemical potential Neutron decays occur until µ χ = µ B . Therefore, the maximum N χ can be determined by equating the two chemical potentials. Choosing for m χ the smallest value allowed by the lifetime of 11 Be, m χ = 939.06 MeV [94], and for α the smallest value which can still be comparable in its effects to those of gravity, α = G N m 2 χ , we obtain M χ /M b ≤ 0.36. Here, we have used the neutron star parameters given in Table I, and we have assumed the radius of the dark core to be similar to the radius of the neutron star.
The above constraint on M χ /M b corresponds to α = 0.071, which would be marginally detectable in LIGO/VIRGO data, see Fig. 1. Such a modification could be more clearly observable by comparing LIGO/VIRGO measurements with results from a spacebased experiment like LISA, see Section II. Note that the upper bound on α in the neutron decay scenario has been calculated for the choices of m χ , α which maximize the potential effect on inspirals; either increasing m χ or increasing α leads to a tighter upper bound on the mass ratio M χ /M b . An increased α thus leads to a decrease in the dark charge, so that overall the strength of the dark force between the binary stars is decreased. Note also that, in scenarios where dark particles are produced mainly via exotic neutron decay, there is no way to create an uncharged neutron star, so the dominant effect of the dark force is always a modification of the potential in which the neutron stars orbit as shown in Fig. 1, never dipole radiation of the force mediator.
We finally note that the coupling that induces decays of the form n → χ + X (where X can, but does not have to be, the dark force mediator A ) will always lead to a coupling of neutrons to A at 1-loop level. This shows that the assumption of a dark force coupling only to hidden sector particles must be relaxed to produce an observable effect in the inspiral waveform. Then, however, new physics constraints from modified gravitational wave signals need to be discussed in the context of constraints arising from more generic fifth-force searches [23][24][25]. The interplay between these various constraints is model-dependent and would require dedicated studies.

IV. CONSTRAINTS ON ATTRACTIVE FORCES
We now switch gears and consider attractive dark sector forces. In this case, important constraints arise from black hole formation inside neutron stars. Additionally, the gravitational binding of DM particles to the neutron star is now small compared to their attractive selfinteractions, leading to the possibility that the dark core can migrate out of the neutron star and stop affecting the observed inspiral dynamics. In the following we provide analytical estimates of these constraints and then present the full numerical results in Figs. 5 and 6.

A. Black hole formation
As a neutron star accumulates a large number of DM particles close to its center, it runs the risk of the dense dark core collapsing into a black hole which eventually consumes the whole star [37,44]. The mere existence of old neutron stars tells us that black hole formation cannot be too efficient, and this in turn restricts the DM parameter space [95,96]. The most optimistic scenario (i.e. the one in which black hole formation is least likely) assumes that the DM particle χ is a fermion. In this case, Pauli repulsion helps to stabilize the dark core at radii larger than the Schwarzschild radius.
For the following estimate we follow Refs. [34,37,44], but assuming that the range of the dark force is always larger than the radius R χ of the χ core. Using the virial theorem (see Section III B), we determine if a stable radius exists for given values of α , m χ , and N χ . In contrast to the case of repulsive forces, it is possible for the DM particles to be either non-relativistic, fully relativistic, or in between. Inserting the Fermi momentum p F = (9πN χ /4) 1/3 /R χ from Eq. (24), and the potential energy given by Eq.
This expression is of course fully analogous to Eq. (26), and the different terms in it can be visualized in Fig. 3, except that now the dark force is attractive, i.e. the cyan lines should now be solid rather than dashed. If, for a given N χ , a solution for R χ exists, the dark core is stable. If no solution exists, the core will collapse into a black hole. We have checked that in the parameter regions of interest to us, the solution for R χ , if it exists, is always larger than the Schwarzschild radius.
For the parameter ranges that we are most interested in, namely those featuring the largest possible N χ and α , minimum stable R χ values typically occur at relativistic p F . In this case, Eq. (35) can be simplified to the constraint This limit is saturated if the baryonic gravitational potential is completely negligible. If we assume that the maximum value of N χ is determined by DM capture from the halo (see Eq. (28)), we can turn Eq. (36) into a constraint on α . Using the definition of α from Eq. (2), we obtain Here, we have used the neutron star parameters from Table I for both neutron stars.
The parameter values where no stable dark cores can exist, i.e. where Eq. (35) has no solution, are shown in purple in Figs. 5 and 6. In Fig. 5, no assumptions are made on how the dark particles came to be within the neutron star. In Fig. 6, the same constraints are shown as a function of m χ and α , assuming that the dark core has been accumulated via capture. The two panels correspond to different assumptions on the capture cross section σ nχ . We observe from Fig. 6 that at m χ 1 GeV, larger m χ implies more stringent constraints on α . This is because of the lower DM number density at large m χ , which implies smaller total dark charge Q χ and thus a smaller attractive force. Of course, Fermi repulsion is also weaker when less DM is accreted, however the corresponding contribution to the virial equation (35) scales less strongly with m χ than the dark force potential. At m χ 1 GeV, the purple exclusion curve and the α contours in the left panel of Fig. 6 flatten because of Pauli blocking, which limits the number of final states available to the neutrons participating in scattering processes and thus limits the number of target particles that DM particles can scatter off and be captured. Constraints on an attractive dark force from binary neutron star systems as a function of the dark fine structure constant α (expressed here relative to the gravitational force) and the amount of DM bound to the neutron stars (assuming both stars carry the same amount of DM). The four panels correspond to different choices of the DM mass mχ, as indicated in the plots. The dotted contours show the value of α (in the case of neutron stars with equal dark charge-over-mass ratio) or β (for the case that only one neutron star carries a dark charge, see Eq. (2)) at each point in the plot. We indicate in red the estimated minimum value of α or β required for the dark force to leave an observable imprint in LIGO data. In the purple region, the dark force induces collapse of the neutron stars' dark cores to black holes. Inside the gray shaded region, the dark cores are expelled from the neutron stars as soon as the screening of the force is lifted. We show this region for m −1 med = 100 km (solid), 1000 km (dot-dashed), and 1 × 10 5 km. Bounds shown in black are model-dependent: in the region above the black-dashed line, the DM population inside the neutron cannot be generated by capture from the halo, based on the geometric arguments that lead to Eq. (28). Above the solid black line, the hypothesis that the DM population is produced radiatively during the supernova that created the neutron star is not tenable because the energy available in a supernova (as measured in SN 1987A) is not sufficient. Within the yellow regions, the DM particles would be stripped away from their host neutron star before they can have a significant impact on the inspiral. The benchmark values of m −1 med are the same as for the gray DM expulsion bound. The neutron star parameters are chosen in accordance with Table I, and without loss of generality we have chosen qχ = 1.
The bounds from black hole formation shown in Figs. 5 and 6 have been derived assuming the dark core is a degenerate Fermi gas, i.e. that Fermi repulsion is maximal. It is possible that the dark core collapses to form a black hole before the required number of particles are captured to satisfy the degeneracy condition. This would lead to a more stringent bound in comparison to Figs. 5 and 6. We have also neglected the question of whether the black hole formed would indeed consume the entire neutron star in finite time or if it evaporates. If the black hole evaporates before consuming the star, no residual dark core would be left, so in this case no modifications to the inspiral dynamics are expected. However, using the results of ref. [97], we conclude that evaporation is faster than the collapse of the neutron star only if M χ 1 × 10 −20 M .
Solving Eq. (35) also yields the radius of the dark core when such a stable radius exists. For light DM mass, Fermi repulsion results in core radii R χ larger than the neutron star radius R. For very light DM, R χ becomes even larger than m −1 med , so that the dynamics of the dark core and the binary inspiral would be very different to that of Section II. The parameter region where this happens is shaded in yellow in Figs. 5 and 6. As previously discussed in Section III B, these yellow curves correspond to the maximum amount of charge that fits within the domain of influence of the neutron star. Particles at radii larger than m −1 med will be efficiently stripped off the neutron star.

B. Expulsion of the Dark Core
We can estimate whether the dark cores of two inspiraling neutron stars remain trapped at the center of their host stars, or whether they are ejected from them by virtue of the strong dark sector force acting on them.
If the dark cores are ejected as soon as the dark force between the two neutron stars becomes unscreened the observed signal will be unchanged compared to the predictions of pure general relativity.
Consider the case where the dark core of one of the neutron stars is about to be expelled from the star and is located at its boundary. We compare the gravitational force experienced by the DM particles (which is strongest at the boundary of the neutron star) to the opposing dark force induced by the dark core of the other neutron star in the binary system: Here the indices i and j label the two neutron stars (i = 1, j = 2 or i = 2, j = 1). As before, R i are the neutron star radii, ∆ is the distance between them, M b,i and M χ,i are the total masses of their baryonic and dark constituents, respectively, M i = M b,i + M χ,i are their total masses, and Q i are their dark charges. The dark core is expelled from the neutron star if |F grav,i | < |F dark,i |, i.e. if In the regime m med ∆ 1, and for DM particles carrying unit charge under the dark force (q χ = 1), this condition simplifies to In the last equality, we have assumed that both neutron stars carry equal amounts of DM, and we have used the definition of α from Eq. (2). This implies that the gravitational potential can be overcome already at very large distances, leading to expulsion of the dark cores long before the actual merger. In Figs. 5 and 6 we show the bounds arising from requiring that the expulsion distance be less than the force range as gray contours. These contours are calculated without neglecting Yukawa screening effects. We consider three different values for the range of the dark force, m −1 med : at m −1 med = 100 km, the onset of the dark force would happen within the LIGO/VIRGO frequency band; at m −1 med = 1000 km, the dark force would be suppressed while the binary system emits gravitational waves in the LISA frequency band, but unsuppressed when it emits in the LIGO/VIRGO band; at m −1 med = 1 × 10 5 km, the dark force would be important for both LISA and LIGO/VIRGO.

C. Constraints on Dark Core Production Mechanisms
As in the case of a repulsive dark force (Section III), further constraints exist on specific dark core production mechanisms. In Fig. 5, we indicate by the horizontal black dashed lines the maximum amount of dark matter that can be accreted from the halo assuming a geometric capture cross section (see Section III C and Appendix A). Figure 6 also assumes DM particle capture, with either a fixed DM-nucleon scattering cross section (left panel) or a geometric capture cross section (right panel). The greatest possible value of α that can be realized assuming the dark core is produced via accretion from the halo (Eq. (28) Producing dark sector particles in neutron decay is not a viable option for attractive dark forces. In the absence of such a dark force the neutron star equation of state is already too soft to allow for 2M neutron stars [88,89], and an attractive force will only further soften it.
The possibility that the dark core is produced radiatively during the supernova explosion that creates the neutron star is constrained in the same way as for repulsive dark forces, see Section III D. This constraint is shown in Fig. 5 as a solid horizontal black line. The greatest possible value of α for radiatively produced DM is 0.23. It occurs for m χ = 0.25 GeV under the assumption of perfectly efficient production in the supernova. Note again that these masses are well above the supernova temperature, thus it is implausible to expect a large fraction of the supernova energy to be converted to dark particles in order to achieve this upper-limit value.
In summary, we conclude from Fig. 5 that dark sector forces cannot significantly affect gravitational wave observations of neutron star inspirals as long as we assume that the neutron stars acquire their dark charge during or after their creation in a supernova. In principle, this no-go theorem could be avoided in tiny regions of parameter space (upper left corner of the plots in Fig. 5) by a mechanism that endows neutron stars with dark cores that account for at least ∼ 1% of their mass. Such a massive DM core could be realized if the DM is already collapsed into a very dense structure, and subsequently the entire DM core either seeds the formation of a star or is captured en masse by a star or neutron star in what must be an incredibly rare event. Of course, such a mechanism is unavailable to repulsively-interacting DM, as the dark sector interactions would prevent such collapse.
Both of these scenarios are constrained by the requirements that the collapse of pure DM systems not damage the observed overall halo structure of galaxies. In particular, collapsing and asymmetrically charged DM objects should only account for a small fraction of the DM in the Universe, making the envisioned capture events even more rare. A detailed analysis of the abundance bounds in this case is left to a future publication. Assuming nonetheless that the neutron star consists of equal parts DM and baryons, i.e. M χ ∼ M b , the maximum value of α ∼ 0.3 occurs for m χ ∼ 1 GeV and α ∼ G N m 2 χ , with the amount of DM particles just below that which would cause collapse to a black hole.
Since compact dark objects that have collapsed in the early Universe are scarce, it is likely in the above scenarios that a neutron stars carrying a large dark charge forms a binary system with an uncharged partner. Therefore, in these scenarios, it is likely that the dominant signature of the dark force would be due to dipole radiation (see Section II B).

V. CONCLUSIONS
We have investigated the conditions under which a new long-range force acting on dark matter can affect the dynamics of neutron star inspirals and can lead to observable modifications to gravitational wave signals observed in detectors such as LIGO/VIRGO and LISA. This scenario is based on the assumption that a large population of DM particles exists within neutron stars.
For repulsive forces, a crucial constraint arises from the requirement that the dark core of the neutron star remains stable (i.e. that gravitational attraction remains stronger than the repulsive dark force). Moreover, such scenarios are severely limited by the fact that a suffi-ciently large dark core cannot be produced via DM accretion from the halo; the greatest effect on inspirals which is possible through accretion has strength normalized to gravity of α = 2 × 10 −15 . Considering more exotic production mechanisms for the DM particles, we conclude that the strongest possible signal in the case of a repulsive dark force corresponds to an α value of 0.071, and arises through neutron decays similar to those which might explain the neutron lifetime measurement disagreements.
For attractive dark forces additional constraints arise from the possible collapse of a massive dark core to a black hole that would consume the whole neutron star. In addition, a strong attractive force could lead to the expulsion of the dark cores from their host neutron stars in a binary system. In this case, the dark cores and their self-interaction cannot contribute significantly to gravitational wave signals any more. A similar effect is the stripping of DM particles from the neutron stars in a binary system. In fact, if the DM population bound to a neutron star consists of relatively light particles, it may extend to radii much larger than that of the baryonic matter. Such DM particles are easily freed from the star's binding force and thus stop contributing to the inspiral dynamics.
In the case of attractive dark forces, the maximum effect on the inspiral dynamics and the gravitational wave signals also depends on the mechanisms by which the dark cores of the neutron stars form. For accretion from the halo, the maximum possible effect is characterized by α = 10 −14 . If the DM population is produced in the supernova that also gives birth to the neutron star, the strongest possible signal is given by α = 0.04. Even stronger modifications to gravitational wave signals are possible up to α = 0.3. This value, however, can only be realized if the DM particles have either seeded the formation of the neutron star's progenitor star or have been accreted as an already collapsed core in a very rare process.
We conclude that, if any deviations from the predictions of general relativity are found in gravitational wave signals from neutron star inspirals, new exotic mechanisms of DM production in neutron stars are required. This would likely herald the existence of large, compact DM structures that either seed star formation or are captured by stars or neutron stars later in their lives.
We finally note that, beyond the observables discussed in this paper, the dynamics of the dark sector could also affect gravitational waves from binary neutron star inspirals in very different ways. Namely, rather than altering the force between the two inspiraling bodies, DM particles could instead alter the properties of the neutron stars. This includes, for instance, alterations in their tidal deformability due to the presence of dark sector particles [21] or alterations arising from particles produced by black hole superradiance effects [98]. Neither of these effects would contradict the results presented here for probing dark sector forces.
The capture rates for dark particles in the scenarios under consideration in this paper have been calculated in refs. [34,44,99,100]. Here, we summarize these calculations, fixing the notation and highlighting the various assumptions made. We emphasize that, while we have utilized standard DM halo parameters throughout this article as a benchmark, our results do not require that the particles captured in neutron stars make up all the DM in the Universe.
Firstly, we approximate that the average DM velocity, escape velocity, and baryon density are uniform throughout the body of the neutron star and that the DMnucleon scattering cross section is velocity-independent. Then, the capture rate of DM on a neutron star is The parameters and functions that enter this equation and have not already been defined are: • v(r) is the infall speed of the dark particles, which we approximate by the escape velocity at the surface of the neutron star. For compact objects like neutron stars, this is a good approximation because the kinetic energy a dark particle acquires during infall is much larger than the initial kinetic energy at r = ∞.
• ξ is the fraction of neutrons that contribute in the scatting after including the effects of Pauli blocking. ξ = 1 for DM masses above a GeV, but for smaller values where the momentum of the DM particle becomes smaller than the neutron Fermi momentum p F = (3π 2 ρ n /m n ) 1/3 (given here in terms of the neutron mass density ρ n and the neutron mass m n ) the fraction of neutrons which have enough momentum to have an accessible final state to scatter in to is Here m r is the reduced mass of the neutron-DM system, m r = m χ m n /(m χ + m n ).
For the saturation cross-section we use the result of Ref. [29] σ sat where N b is the number of baryons in the neutron star.
• B is a function accounting for the minimum energy loss necessary to capture a DM particle: This function is smaller than one only when the DM mass is 10 6 GeV.
Using C χ from Eq. (A1), the number of DM particles captured will be Note that we here neglect the possibility of DM selfannihilation, co-annihilation [44] or semi-annihilation [101], which would reduce N χ .