Probing Muonic Forces with Neutron Stars Binaries

We show that gravitational wave emission from neutron star binaries can be used to discover any generic long-ranged muonic force due to the large inevitable abundance of muons inside neutron stars. As a minimal consistent example, we focus on a gauged U(1)$_{L_\mu - L_\tau}$ symmetry. In pulsar binaries, such U(1)$_{L_\mu - L_\tau}$ vectors induce an anomalously fast decay of the orbital period through the emission of dipole radiation. We study a range of different pulsar binaries, finding the most powerful constraints for vector masses below ${\cal O}(10^{-18} {\rm eV})$. For merging binaries the presence of muons in neutron stars can result in dipole radiation as well as a modification of the chirp mass during the inspiral phase. We make projections for a prospective search using the GW170817 event and find that current data can discover light vectors with masses below ${\cal O}(10^{-10} {\rm eV})$. In both cases, the limits attainable with neutron stars reach gauge coupling $g^\prime\lesssim 10^{-20}$, which are many orders of magnitude stronger than previous constraints. We also show projections for next generation experiments, such as Einstein Telescope and Cosmic Explorer.


I. INTRODUCTION
New long range interactions give rise to distinctive signatures in a wide range of observables. Such interactions are, however, strongly constrained by fifth-force tests [1][2][3][4], unless they are either screened or the couplings to the first generation fermions are suppressed. We show that the dynamics of neutron star (NS) binaries provide ideal laboratories to probe long range muonic forces due to the significant abundance of muons inside them (by mass 0.1% M ). The observations of the neutron star merger event, GW170817 [5], and various pulsar binaries [6][7][8][9][10][11][12][13][14][15] give us the opportunity to probe these new exotic forces. These methods of probing muonic forces via neutron star binaries are completely general, applicable to both vector and scalar mediators. However, as a concrete realization we focus on a long-range gauged U(1) symmetry.
Additional U(1) gauge symmetries with masses below the weak scale are simple extensions of Standard Model that can act as a mediator to the dark sector (see e.g., [16][17][18][19]), are common predictions of string theory [20], and can explain experimental anomalies [21][22][23][24][25][26]. The observed matter content limits the linearly independent conserved currents to B − L, Hypercharge (equivalent to kinetic mixing), and (up to neutrino masses) L e − L µ , L µ − L τ , and L e − L τ . 1 The small number of possibilities highlights the need to find all experimental ways to probe these light vectors. While most of the focus when studying the phenomenology of new U(1) gauge symmetries has been above the eV scale, it is interesting to a jdror@lbl.gov b ranjan.laha@cern.ch c toby.opferkuch@cern.ch 1 Relaxing the requirement of anomaly cancellation with just the Standard Model fields greatly enlarges the possibilities but is highly constrained from searches for flavor changing neutral currents [27][28][29].

K E f E g m c C h b l r E A x n 9 A p G b O Q O E r N H T v C F w 6 B i V r V l s P / t F 4 A o 2 o / 5 K 4 f A H N p b I m 1 U S A w e H h Z F g + 5 Z B T E P C a E S h 4 n Y z o
h k l C I X / I n Z Z Z 8 I c r F N c 3 V U u u k X S q a 5 W L p s Z K v X a W F s + g M n a M C M t E l q q E 6 a q A W o u g F L d A n + t I C 7 V V 7 0 9 5 / r R k t 3 T l F f 6 B 9 / A B x c p E l < / l a t e x i t > X = NS, BH, WD, Star < l a t e x i t s h a 1 _ b a s e 6 4 = " M 9 Q d M l K W 5 a t 9 T V Y v p B f C W C 9 s 8 D 8 = " > A A A C G 3 i c b Z B L S w M x E M e z 9 V X r q + r R S 7 A I H k r Z r Y J e l K I e e p J K 7 Q P a U r J p 2 o Z m H y S z Y l n 2 e 3 j x q 3 j

Muonic Vector
< l a t e x i t s h a 1 _ b a s e 6 4 = " C F f P 7 3 G P X Q u 8 6 E P L + 5 b P x g i 4 l m 0 study the constraints for lighter masses, where forces become long-ranged. For B − L, L e − L µ , and L e − L τ , the constraints below this scale become extremely powerful from fifth force tests [1][2][3][4] which constrains such forces to be weaker than gravity once the vector mass drops below O(10 −4 eV) and constraining the gauge coupling g 10 −20 at the lowest masses. Kinetically mixed dark photons do not experience such constraints due to screening of the charge between protons and electrons leading to rich phenomenology at low masses (see e.g., [30] and references therein). Interestingly, L µ − L τ forces also do not experience such constraints since the muon fraction in ordinary matter is negligible.
cleosynthesis [48,49], from SN1987A [50] 2 and neutrino self-interactions [52][53][54] constraining, g 10 −5 . Such searches are inherently much weaker than searches for long range forces as they do not scale with the size of the apparatus. Alternatively if there exists a mass mixing between the Standard Model Z boson and the new vector, L ⊃ ε Z m 2 Z X µ Z µ , it can induce a long range froce that would have been observed in neutrino oscillations unless ε Z g 10 −52 [41]. However since a mass-mixing is forbidden by gauge invariance, it is naturally small and constrained to be very small [55].
Neutron stars are fascinating astrophysical objects exhibiting physics across many different length scales. These objects are born through stellar collapse where the production of a cold degenerate gas of neutrons, via the process p + e − → n + ν e , supports the star from further collapse. Their existence, however, is contingent entirely on the stability of the neutron through sufficient Pauliblocking of the inverted process n → p + e − + ν e [56][57][58][59][60]. As the neutron density increases similar processes involving muons, rather than electrons, become energetically favorable. This leads to the production of a significant number of muons and forbids their subsequent decay [61]. Neutron stars with masses of order a solar mass subsequently have 0.15-0.75% of their mass stored in muons, providing a unique laboratory to test couplings of muons to light new degrees of freedom. This has been leveraged 2 The robustness of the supernova bounds has recently been called into question in [51].
to place constraints on muon-philic dark matter due to its accretion in neutron stars [44,62,63].
A key feature of this muon population is their asymmetric nature, i.e., the production of only muons and not anti-muons. Consequently, the presence of new longrange forces coupled to muons leads to neutron stars acquiring large effective charges. The coupling we are interested in constraining is that of L µ − L τ : where V α and g denotes the new vector and its coupling strength. We will rely on the couplings to muons as neutron stars contain negligible numbers of taus or neutrinos. Extending the constraints we derive here to a scalar force is straightforward, where O(1) changes are expected in the limits on the coupling relative to the vector case. The observation of gravitational waves (GW) from NS-NS mergers as well as timing measurements of binary pulsars provide exquisite sensitivity to not only these scenarios but also many types of dark matter candidates (see [64] and references therein). In both cases GWs are the dominant energy loss mechanism required to describe the dynamics of the system. For the case of mergers, the sensitivity arises from the long observation duration of the post-Newtonian stage of the GW signal [5]. While for the case of binary pulsars, precise measurements of the change in pulsar period yield accurate determinations of both the relativistic corrections to the binary orbits and the orbital decay, (see e.g. [65]). Both of these scenarios are therefore sensitive to additional energy loss mechanisms through the emission of the light vector, while neutron star mergers are also sensitive to an additional force sourced between the muon content of the two neutron stars, see Fig. 1. These effects have been studied in the context of new long-range forces in a hidden sector [66][67][68][69][70][71][72][73][74]. In what follows we consider the impact a gauged U(1) Lµ−Lτ symmetry can have on neutron star binaries.

II. MUONS IN NEUTRON STARS
The presence of muons in neutron stars arises due to chemical equilibrium and charge neutrality maintained via the processes n → p + e − /µ − + ν e/µ and e − → µ − + ν e + ν µ , referred to as beta equilibrium. The existence of muons in neutron stars follows from estimating the Fermi energy of a neutron gas, E F = (3π 2 n n ) 2/3 /2m n (where m n and n n are the nucleon mass and number density of the neutrons respectively). Neutron star masses around a solar mass and radius of 10 km give typical Fermi energies ∼ 100 MeV, suggesting a significant muon abundance.
The exact muon content of a neutron star depends upon the QCD equation of state (EoS) relating the energy density to the pressure in the interior of the neutron star. We base our estimates of the muon content on the most recent Brussels-Montreal EoS [75], which is an update of older works based on the two-and three-nucleon force calculations of [58,75]. These EoS are all compatible with recent limits on the tidal deformability constraints from GW170817 [5,76]. Note however, that two of these EoS, BSk22 and BSk26, are disfavoured due to neutron star cooling measurements. Measurements suggest that relatively few neutron stars exhibit large cooling rates associated with the direct Urca process [77,78]. The former EoS admits direct Urca processes for neutron star masses in excess of 1.151M , and therefore anomalously large cooling rates in the majority of the neutron star population [79], while the latter EoS cannot support these processes at all which is in conflict with recent results in [80], suggesting the presence of this process in the neutron star MXB 1659-29. Nevertheless, we take the envelope (shown in Fig. 3 of the supplementary material) of all four EoS from [75] as a conservative estimate of the neutron star muon content as a function of the neutron star mass relevant for the different binary systems we consider. Finally, we note that similar EoS are expected to hold for the description of rapidly rotating neutron star pulsars, while such rotating neutron stars in binary mergers are expected to aid in the extraction of the underlying EoS [81].

III. NEUTRON STAR BINARY MERGERS
A new muonic force can have a dramatic effect on NS-NS and NS-black hole (BH) binaries. In the absence of exotic forces and in environments where dynamical friction can be neglected, the dynamics of these inspirals are determined by gravitational attraction and emission of gravitational waves. Any new exotic force changes the dynamics in two different ways: (i) the Yukawa force between the muon cores can accelerate or decelerate the merger, and (ii) the emission of the mediator particle increases the energy loss of the system and accelerates the merger. For concreteness, we consider a repulsive force mediated by a vector boson and follow the techniques advocated in [69]. In this section, we outline the technique via which we derive our new constraints. We postpone the detailed formulae to the supplementary material.
If the muonic charges carried by the two astrophysical objects are denoted by q 1 and q 2 , then the Yukawa force between them can be written as [68] where α ≡ g 2 q 1 q 2 /(4πG N m 1 m 2 ) =q 1q2 > 0 and r denotes the distance between the two astrophysical objects. The masses of the neutron stars are denoted by m 1 and m 2 , and the mass of the mediator vector boson is denoted by m V . The presence of such a new force modifies Kepler's law and the total energy of the system, E tot . The energy loss rate of the system, dE tot / dt, is then determined by the energy loss via gravitational waves, dE GW / dt, and the energy loss via the emission of the new vector particle, dE V / dt: is the charge-to-mass ratio. Due to the presence of this exotic force, both the plus and the cross polarizations of the GWs are affected. We analytically calculate the GW amplitude and its phase to first-order in α and γ. We add post-Newtonian corrections following [82].
To derive the upper limits on α and γ, we follow the standard Fisher information matrix analysis [83][84][85][86]. The frequency range used in the Fisher information matrix is denoted by [f low , f high ]. We choose f low following [69] and f high = f ISCO , where f ISCO denotes the frequency of the innermost stable circular orbit [68]. While deriving the limits on α and γ from the observation of GW170817, we use the aLIGO sensitivity curve. We also show projected constraints of the limits that can be achieved by Cosmic Explorer and the Einstein Telescope. For all these cases, we use the analytical form of the sensitivity curve as provided in [69].

IV. BINARY PULSARS
Binary pulsars are a powerful probe of ultralight vectors. In a binary system, the motion of the pulsar and its companion (which in rare cases can also be a pulsar) is imprinted in the pulsar time-of-arrival data as an oscillation with a period, P b , which is typically of O(days). This is much larger than the pulsar period, which is O(msec − 10sec). While simple Keplerian orbits can explain the qualitative motion of a binary pulsar system, the pristine precision of pulsar measurements allow the detection of deviations from classical mechanics due to a few different relativistic effects.
The deviations of a binary pulsar system from simple orbital motion can be described in terms of post-Keplerian parameters [87,88] (see also [65] for a review), the periastron precession,ω, the combination of gravitational redshift and Doppler shift, γ d , as well as the secular drop in the binary period,Ṗ b , which is O(10 −12 ). The three parameters depend on different combinations of the pulsar and companion masses. Since measurements ofω and γ d typically carry a much smaller uncertainty thaṅ P b , it is natural to use these to fix the masses and usė P b to set constraints on new physics as advocated in [89] (alternatively, one can use the so-called γ d −ω −Ṗ b test where one allows for deviations in all 3 parameters and varies the masses however assuming that the new physics deviations are small this does not make a meaningful difference in the results).
For a binary pulsar system, the large muon abundance leads to the emission of the light L µ − L τ vector. 3 The subsequent rate of change in the energy of a pulsar relative to gravity is given by [89] Ė where the elliptic correction functions for an eccentricity, , are: where J n is the Bessel function of n-th order and J n is its derivative, while the sum begins at n 0 = m V P b /2π. The observable for each binary pulsar is the ratio of the intrinsic change in the binary orbital period, 4 P int b , to the prediction from GR, P GR b . Rewritten in terms of the energy ratios from above yieldṡ 3 In addition there are relativistic corrections to the force between two neutron stars which are prominent at slightly higher vector masses than radiation, however these are subdominant to constraints from ensuring the new force does not eclipse gravity during neutron star mergers. 4 Note that P int b is the observed value of the period drop which requires subtraction of the effects due to galactic rotation.
We are now in a position to set constraints using current pulsar systems. There are many binary pulsars which have by now observed gravitational radiation including pulsar-(non-pulsating) neutron star binaries [8,[13][14][15]90], a pulsar-pulsar binary [6], pulsar-white dwarf binaries [7,[10][11][12], as well as a pulsar-Oe-type star binary [9]. In setting limits different systems have different advantages. Since the energy ratio in Eq. (4) is proportional to P 2/3 b , large orbit binaries are ideal at probing low vector masses. On the other hand, binaries can only emit radiation efficiently if m V 2π/P b leading to smaller binaries becoming effective at probing larger m V . Furthermore, we note that in order to have significant emission, the binary must carry a dipole moment, requiring the pulsar and its companion to have differing charge-to-mass ratios (as is apparent from the linear dependance of Eq. (4) on the charge-to-mass ratio). This quantity is maximized in pulsar-white dwarf or pulsarvisible star binaries where the charge of the white dwarf is negligible.
To set our constraints on gauged L µ − L τ using pulsar binaries we take the 2σ limit onṖ int b /Ṗ GR b provided by experiments with the parameters summarized in Fig. 5 of the supplementary material.

V. RESULTS AND DISCUSSION
The constraints on a gauged U(1) Lµ−Lτ derived in Sections III and IV are shown in Fig. 2 for the m V -g plane. The lines marked NS-NS merger, show the estimated sensitivity that could be achieved with a dedicated LIGO/VIRGO analysis using the observation of the neutron star merger, GW170817. The sensitivity curves in blue rely on dipole emission of the vector while the green curve arises from the new force between the two neutron stars. The right-hand boundaries of both these constraints are functions of f ISCO . These boundaries lie at different L µ − L τ vector masses due to the different dependencies of the vector emission (force) on the angular velocity (binary separation) and therefore frequency of the binary. The dark blue shaded region indicates parameter space where the GW170817 system would be unbounded. The dipole emission constraint asymptotes to a constant values at small vector masses as the only dependence on the vector mass arises in the step-function [see Eq. (B12)].
The constraints shown in red are obtained using current binary pulsar data. We show the envelope of the constraints from all pulsars as the shaded purple region in Fig. 2, while we defer the results for individual pulsars to Fig. 5 in the supplementary material. For eccentric binaries, the typical distance between the binaries changes significantly across the orbit and allows vector emission across different masses leading to the observed step like pattern. Given the pulsar binary separations, the dipole radiation emission is only active for vector masses below O(10 −18 eV). The sensitivity to the gauge coupling is weaker than the equivalent constraints from NS-NS mergers, however they do not require additional analysis by LIGO and are a present constraint. As such the binary pulsar constraints serve as a robust alternative to the GW measurements from coalescing neutron star binaries.

VI. CONCLUSION AND OUTLOOK
We demonstrate the discovery reach to a new muonic force using neutron stars (NS) binary systems. The significant muon abundance inside neutron stars leads to two new effects: (i) dipole emission of the force mediator, and (ii) an additional force between binaries comprising of two neutron stars. These effects lead to changes in the dynamics of both inspiraling neutron stars and pulsar binaries. Measuring the deviation in the gravitational waveform and the pulsar period, respectively, are powerful tests of the presence of these new forces. Based on current data and assuming a U(1) Lµ−Lτ force as an example, Fig. 2 shows the discovery reach to gauge couplings as small as O(10 −20 ), orders of magnitude better than current probes.
Given the estimated sensitivity of current gravitational wave experiments to these scenarios, we advocate for dedicated studies using improved calculations of the gravitational wave waveform confronted with data from the measured NS-NS merger event GW170817. In addition, future detection of BH-NS events would allow for the isolation of the vector emission.

ACKNOWLEDGMENTS
We thank Kfir Blum, Raghuveer Garani, Admir Greljo, Edward Hardy, Joachim Kopp, Marko Simonovic and Yotam Soreq for useful discussions. In particular we would especially like to thank Evan McDonough for his hugely helpful correspondence. This work was initiated and performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY1607611. TO has received funding from the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 637506, 'νDirections') awarded to Joachim Kopp. JD is supported in part by the DOE under contract DE-AC02-05CH11231.
Note Added: During the preparation of this work, we became aware of [91] which also considered using the muon content in neutron stars to set constraints on ultralight vectors using pulsar timing data. the relationship between the pressure and density of the neutron star. Based on these fitting functions for the EoS (see Eq. (C4) in [75]), the Tolman-Oppenheimer-Volkoff (TOV) equations can be solved yielding a relation between the mass and density of the neutron star as a function of the radius. To determine the muon abundance, charge neutrality is assumed as well as equilibrium between the muons and electrons, resulting in chemical potentials that are the same µ e = µ µ . Here Y i is the abundance i = p, e, µ of protons, electrons and muons defined as Y i ≡ n i /n where n is the total density and n i is the density of the species in question. From the chemical potential, under the assumption of the muon and electron gases being degenerate we have where m e and m µ are the electron and muon masses, respectively. This yields the muon number density To evaluate this expression we take the number density of electrons as a function of neutron density from Eq. (C17) in [75], allowing for the total mass of muons M µ inside the neutron star to be determined for the different EoS. The results of which are shown in Fig. 3. We observe that neutron stars with masses greater than a solar mass have M µ ≥ 1.5 × 10 −3 M NS , while a two solar mass NS would have muon content in the range 0.7 × 10 −3 M NS ≥ M µ ≥ 0.24 × 10 −3 M NS . For reference we also show the two heaviest observed neutron stars; the NS-WD binary PSR J0348+0432 and the millisecond pulsar J0740+6620 [92].

Appendix B: Modifications to the gravitational waveform
Given the force between the two astrophysical objects as denoted by Eq. (2), one can derive the orbital frequency, ω, of the system: where the negative sign before α denotes that the force is repulsive. While writing the above equation, we neglect the spins of the astrophysical objects, and assume them to be point objects. The orbital frequency is related to the frequency of the gravitational waves, f GW , as f GW = ω/π. The total energy of the system can be written as where µ denotes the reduced mass of the system. The energy loss rate due to the emission of gravitational waves can be written as The energy loss due to the radiation of a light vector particle can be written as where γ is defined below Eq. (3). Due to the presence of this extra new force and a new way to lose energy from the system, the gravitational wave signature from the merger of two compact astrophysical objects change. The two polarizations of the gravitational waves can be written as [69] h The inclination angle is denoted by ι, and the time and phase of the system at coalescence is denoted by t c and φ c . The orbital phase of the system is denoted by φ and we will mention how to calculate it later. The amplitude of the signal is denoted by A(t) where where the luminosity distance to the source is denoted by D L . Given the detector responses, F + and F × , the strain detected by the detector is given by where z = 2 (φ c + φ(t − t 0 )) and the time in the detector frame when the coalescence is detected is given by t 0 .
Defining D eff and φ 0 such that we can write the strain as In order to determine the upper limits on the parameters α and γ from current observations, we need to determine the Fourier transform of h(t). Using the stationary phase approximation and restricting our calculations to first order in α and γ, we get [69] where x ≡ G is defined as where the error function is denoted by erf(x). Finally, the PN corrections from the gravity-only contribution take the form where η ≡ m 1 m 2 /m 2 is the symmetric mass ratio and the coefficients of the sum are given in [82]. We do not include corrections with n > 4 as these additional terms do not significantly alter the presented results. In order to determine the upper limit on the new physics parameters, α and γ, we follow the Fisher information matrix analysis [83][84][85][86]. We denote the dimensionless spin parameters of the two astrophysical objects by χ 1 and χ 2 . We define the symmetric and anti-symmetric dimensionless spin parameter as χ s = (χ 1 + χ 2 )/2 and χ a = (χ 1 − χ 2 )/2 respectively. For the Fisher information matrix, we take the underlying parameters to be θ = {log A, t c , φ c , log M c , log η, χ s , χ a , α, γ} . (B17) Using these 9 parameters, we can construct the 9 × 9 Fisher information matrix Γ whose components are given by where a and b runs from 1 to 9. The 9 parameters, defined in Eq. (B17), are denoted by θ a,b . The inner product is defined as where S n (f ) denotes the spectral noise density of various gravitational wave detectors that we use in our analysis. We use analytical forms of S n (f ) as given in [69] for our calculations. The signal-to-noise ratio, ρ, is given by We define the covariance matrix, Σ ≡ Γ −1 . The rootmean squared error, that can be determined from an observation, for a given parameter θ a is given by the square root of the (a, a) component of Σ. The same expression also gives the 1σ upper limits on the new physics parameters, α and γ: where θ a new denotes α and γ. For the extraction of the 1σ upper-limits we have chosen the parameter values to mimic the observed GW170817 event assuming slowly spinning neutron stars.  This corresponds to the choices m 1 = 1.46M , m 2 = 1.27M with χ 1 = 0.01 and χ 2 = 0.02. While for the effective luminosity distance we use D eff = 40 Mpc. In the left-hand panel of Fig. 4 we show the resulting sensitivity curves for both aLIGO [solid lines] as well as the next generation ground based experiment ET [dashed lines]. In much of the parameter space ET will improve sensitivity in g by at least an order of magnitude. We note that the sensitivity curve and therefore resulting sensitivity for Cosmic Explorer is parametrically similar to ET. In the right-hand panel of Fig. 4 we also show the projected sensitivity in the case of the observation of a NS-BH binary merger. Here we assume that m 1 = M BH = 5M and m 2 = M NS = 1.46M , again assuming small spins of both compact objects (χ 1 = 0.01 and χ 2 = 0.02). This type of merger is particularly sensitive as the charge-tomass ratio is maximized given that that BHs carry zero charge under U(1) Lµ−Lτ . Subsequently the constraints are also less sensitive to the uncertainty in the muon content of the neutron star.

Appendix C: Pulsar binary data
The data used to set the binary pulsar constraints is shown in Fig. 5 (Left). For a given binary the necessary parameters are the binary period, the change in the period relative to that predicted by gravity, the eccentricity, and the masses. We use the neutron star equation of state to compute the muon abundances of the neutron stars as described in the text. In Fig. 5 (Right) we show the constraints from individual pulsar systems with optimistic and pessimistic assumptions on the muon abundance as described in Appendix A.