Probing ultralight isospin-violating mediators at GW170817

Gravitational wave (GW) signals arising from binary neutron star mergers offer new, sensitive probes to ultralight mediators. Here we analyze the GW signals in the GW170817 event detected by the LIGO/Virgo collaboration to impose constraints on the ultralight isospin-violating mediator that has different couplings to protons and neutrons. Neutron stars, which primarily consist of neutrons, are the ideal places to probe the isospin-violating mediator. Such a mediator can significantly alter the dynamics of the binary neutron star mergers, through both the long-range Yukawa force and the new dipole radiation. We compute the gravitational waveform by taking into account the new physics effects due to the isospin-violating mediator and use the Bayesian inference to analyze the gravitational wave data in the GW170817 event. We find that although the current fifth force experiments (including MICROSCOPE and EW) often provide more stringent constraints than the GW170817 data, in the parameter space where the isospin-violating force is completely screened by the Earth (namely, the Earth is charge neutral under this force), the GW170817 data offer the leading constraints: the upper bound on the neutron coupling is $f_n \lesssim 10^{-19}$ in the mediator mass range of $\simeq(3\times10^{-16},\,5\times10^{-14})$ eV.


Introduction
Ultralight bosons appear in a number of well-motivated new physics models beyond the standard model [1][2][3][4].These particles are especially fascinating due to their potential to mediate a new long-range force which is often referred to as a fifth force.Extensive experimental efforts have been devoted to searching for the effects of the fifth force across a wide range of distances and couplings [5].One notable aspect of the fifth force searches is the exploration of the composition-dependent interactions [6], which could lead to violations of the weak equivalence principle (WEP) [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].The isospin-violating mediators, which have different couplings to neutrons and protons, are an intriguing type of composition-dependent interaction in fifth-force experiments and also of great interest in dark matter direct detection [23][24][25][26].Moreover, ultralight isospin-violating mediators can also be searched for through gravitational wave signals in binary neutron star (BNS) mergers.
Thus, we use gravitational wave (GW) data from the BNS merger event GW170817 [27] to probe the isospin-violating force.The GW170817 event, the first BNS merger event detected by the LIGO/Virgo collaboration [27], has a signal-to-noise ratio (SNR) of 32.4,making it an ideal place to study various topics, including the properties of the NSs [28], the equation-of-state (EoS) of nuclear matter [29], and predictions from general relativity (GR) and beyond [30].The GW170817 event was detected by the two LIGO detectors (Hanford and Livingston) and the Virgo detector (Cascina) [27].In contrast, the other BNS merger event, GW190425, was only detected in the LIGO Livingston detector, with a lower SNR of 12.9 [31].Moreover, optical counterparts of the GW170817 event [32,33] and the gamma-ray burst event, GRB 170817A, [34,35] have been detected.The electromagnetic counterparts are instrumental in locating the host galaxy of the BNS system, facilitating the parameter estimation of GW signal [28].Thus, we focus on the GW170817 event in our current analysis.
The ultralight mediator has two major effects on the dynamics of the BNS: First, it leads to a new long-range force so that the orbital frequency of the BNS can be modified significantly.Second, it provides a new radiation channel through which the BNS system can lose energy, if the mediator mass is small compared to the orbital frequency.This then leads to modifications to the GW signals arising from the BNS system; see e.g., Refs.[36][37][38][39][40][41][42][43][44][45][46] for previous studies.Orbital decays due to the radiation of ultralight mediators also have important effects in pulsar inspiraling [40,[46][47][48][49][50][51].To our knowledge, the effects of ultralight isospin-violating mediators on the GW170817 data have not been studied before.To probe the isospin-violating force, we perform Bayesian analysis on the GW170817 data [52], by using the PyCBC inference [53,54].We find that the isospin-violating mediator with mass ≲ 10 −11 eV can be probed by the GW170817 event, which consists of data with frequency up to ∼2000 Hz [27].
The parameter space of the isospin-violating mediators with mass ≲ 10 −11 eV is constrained by a number of experiments, including WEP, motions of asteroids [55,56] and planets [18,57], and black hole superradiance (BHSR) [58].The dominant WEP experiments include the Eöt-Wash (EW) torsion balance experiment [9,10], the lunar laser ranging (LLR) experiments [19][20][21][22], and the MICROSCOPE experiment [11][12][13][14][15].Note that the signals in the WEP experiments depend on the charges of both the test mass and the attractor (gravity source).As emphasized in Ref. [10], measurements with varied test masses and/or attractors are necessary, as charges may vanish for certain test masses or attractors.For example, EW has used Be-Ti and Be-Al pairs as test masses [10], and MICROSCOPE has used Ti-Pt and Pt-Pt pairs [13].While substituting different test masses can be relatively straightforward, changing attractors often presents more of a challenge due to the limited number of nearby celestial bodies.The most commonly used attractors include the Earth [9][10][11][12][13][14][15] and the Sun [10,[19][20][21][22]. Another critical factor in the fifth force experiments is the distance.For example, experiments using the Sun as the attractor probe mediators with mass ≲ 10 −18 eV, as the distance to the Sun is ∼ 1.5 × 10 8 km.In contrast, experiments that use the Earth as the attractor, such as MICROSCOPE [11][12][13][14][15], can probe mediators in a larger mass range, m ≲ 3 × 10 −14 eV, which is determined by the radius of the Earth.Consequently, the scenario where the charge of the Earth is zero and the mediator mass is in the range of ∼ (10 −18 , 3 × 10 −14 ) eV cannot be probed by WEP experiments using the Sun as the attractor.It is thus of importance to employ alternative methods to explore this scenario.
Neutron stars, primarily composed of neutrons, differ significantly in composition from the Earth, making BNS mergers an ideal place to probe the parameter space where the charge of the Earth is zero.Because the Earth consists of mainly heavy elements, its proton-to-neutron ratio is close to unity.In contrast, we find a proton-to-neutron ratio of ∼ (7 − 12)% for NS with mass in the range of ∼ (1 − 2) M ⊙ , if the BSk24 EoS is used [59].Moreover, the separation of the two neutron stars in the GW170817 event is ∼ (20, 400) km, which is smaller than the Earth radius by one to two orders of magnitude, so that it has the potential to probe mediators with an even larger mass range of m ≲ 10 −11 eV.In the mass range of ∼ (10 −16 − 10 −13 ) eV, the dominant WEP constraints come from the WEP experiments that use Earth as the attractor, such as the MICROSCOPE experiment [11][12][13][14][15], and the EW experiment [9,10].We find that the GW170817 data can offer leading constraints for this mass range, if the Earth charge vanishes.
To correctly interpret the GW170817 constraints, as well as other experiments that use celestial bodies as the source of the external field, it is important to properly take into account the effects due to the finite sizes of these celestial bodies when the Compton wavelength of the mediator becomes small compared to the size of the celestial bodies [60].To this end one needs to integrate over the charge density distribution within the celestial bodies.In our analysis, we use the Tolman-Oppenheimer-Volkoff (TOV) equation [61,62] together with an EoS [59] to determine the proton and neutron density distributions for NSs with different masses, and then use these distributions to compute the charge of the NSs for the isospin-violating force.We also provide a simple analytic expression of the Earth charge, correcting a previous erroneous formula in Ref. [14].
The rest of the paper is organized as follows.In section 2 we discuss the effects of the ultralight isospin-violating mediator on the BNS dynamics, and provide calculations of the NS charge by taking into account the charge density distributions in the NS.In section 3 we compute the GW waveform in the presence of the fifth force, and then use it to perform Bayesian analysis of the GW170817 data in section 4. We discuss various fifth force constraints in section 5, including constraints from asteroids and planets data, from WEP experiments, and from BHSR.We present the GW170817 constraints on the ultralight isospin-violating mediators along with other experimental constraints in section 6, where we also discuss two interesting cases: the baryon number case, and the case where the Earth charge is zero.We summarize our findings in section 7. The detailed calculations on the neutron star properties are given in Appendix A.

Ultralight isospin-violating mediator and BNS
In this study, we utilize the data from the GW170817 gravitational wave event [27] to constrain the isospin-violating force.We consider the following interaction Lagrangian where V µ is the vector field that couples to proton p and neutron n with couplings f p and f n , respectively.If V µ is an ultralight particle, it can affect the BNS dynamics in two different ways: First, its long-range Yukawa force can alter the orbital frequency of the BNS significantly.Second, in addition to the gravitational radiation, there is a new dipole radiation channel through which the BNS can lose energy.

Yukawa force and orbiting frequency
If the vector field V µ has a sufficiently small mass, it generates a new long-range force.If the mass m V ≲ 10 −11 eV, the long-range force has a characteristic scale λ = 1/m V ≳ 20 km, thus having the potential to significantly alter the orbiting speed of the binary neutron stars, which is primarily governed by gravity. 1 The new long-range force can increase (decrease) the orbiting speed of the BNS, if it is attractive (repulsive) between the two neutron stars.At the leading order, the total force acting on the two neutron stars is the sum of the Newtonian gravity and the Yukawa-type force mediated by the V µ field: [36,37,40] where G is the gravitational constant, m 1 and m 2 are the masses of the neutron stars, r is the distance between the two neutron stars, m V is the mass of V µ , and α is the parameter that characterizes the relative strength of the new Yukawa force induced by V µ .The parameter α is given by where Q 1 and Q 2 are the charges of the two neutron stars.A positive (negative) α indicates a repulsive (attractive) Yukawa force.
In the presence of such a Yukawa-type force, the orbital frequency of the BNS system is given by the modified Kepler's law [36,37,40]: where M = m 1 + m 2 .The α term resulting from the Yukawa force exhibits a distinct radial dependence as compared to the gravity term, which is crucial in distinguishing between the Yukawa and gravity forces in the waveform analysis of the gravitational wave data.

Dipole radiation
If the orbital frequency of the BNS system is larger than the mass of the light mediator (in natural units), the BNS system can lose energy by radiating the V µ field, thus presenting a new energy-loss channel in addition to the gravitational radiation channel.Analogous to the electrodynamics, the radiation due to the massive V µ field can be expressed as a set of multipole expansion terms [47], with the leading order contribution from the electric dipole radiation 2 The radiation power (energy emitted per unit time) in the electric dipole radiation is given by [36,37,47,48], where ω is the orbital frequency of the BNS system, m V is the mediator mass, µ ≡ m 1 m 2 /M is the reduced mass of the BNS system, and γ is given by (2.6) Note that we have used a zero orbital eccentricity to obtain Eq. (2.5), which is justified since the effect of the orbital eccentricity in GW170817 is negligible [63].
To compute the orbiting frequency of the BNS system, we use the energy conservation to write where E G is the total energy of the BNS system that is predicted by GR, P G is the gravitational radiation power, and E V is the Yukawa correction to the potential energy of the binary system: In our analysis we use the binary energy at the Newtonian level [40] and the quadrupole gravitational radiation power [37,64] Gµ 2 r 4 ω 6 . (2.10)

NS Charge
Both the parameter α, which enters the new Yukawa force, and the parameter γ, which enters the new radiation power, depend on the charge Q of the NS.When the wavelength of the mediator is sufficiently large so that the NS can be treated as a point particle, the charge Q of the NS can be obtained via a simple algebraic summation.However, when the size of the NS is comparable to or larger than the wavelength of the mediator, one has to properly take into account the effects due to the finite size of the NS and the charge distribution within the NS.For a spherically symmetric object with a uniform density, the charge is given by [60] where Q pt is the charge when the object can be treated as a point particle, Φ(x) ≡ 3(x cosh x − sinh x)/x3 , R is the radius of the object, and λ is the Compton wavelength of the mediator. 3It is thus straightforward to obtain the charge of a spherically symmetric object with a charge density distribution: where r is the radial coordinate, and q(r) is the charge density distribution.Note that in the case of constant q, one has Q pt = q(4πR 3 /3).
For the isospin-violating mediator, the charge density distribution q(r) in Eq. (2.12) is given by q where n p (n) is the proton (nucleon) number density.In the long-wavelength limit, the charge of the NS becomes where Z (N ) is the proton (neutron) number of the NS.To obtain the charge density distribution, we first determine the total nucleon number density n(r), by solving the TOV equation [61,62] together with an EoS.We then use the proton-to-nucleon fraction Y p (n), to obtain the proton and neutron distributions.See Appendix A for the detailed analysis.We note that throughout our analysis, we will mainly use BSk24 [59] (the default EoS unless explicitly specified) for the EoS.For comparison, we also consider an alternative EoS, BSk22 [59], in Figs.(2,7,8).

GW170817
In our analysis we focus on the GW170817 event, where the primary NS has a mass of m 1 = 1.46 +0.12 −0.10 M ⊙ and a radius of 11.9 +1.4 −1.4 km, and the secondary NS has a mass of m 2 = 1.27 +0.09 −0.09 M ⊙ and a radius of 11.9 +1.4 −1.4 km [29,65].For a mediator with mass of m V ≳ 10 −11 eV (λ ≲ 20 km), 4 one has to take into account the effects of the charge distribution of the NS.The total proton and neutron numbers for the two NSs are: Fig. 1 shows the proton and neutron charges of the two NSs in GW170817 as a function of the mediator mass.In the left panel figure of Fig. 1, the charges are normalized with respect to the long-wavelength limit Q pt ; the ratio Q/Q pt starts to creases from unity where the mediator mass becomes ≳ 10 −11 eV.In the right panel figure of Fig. 1, the charges are normalized with respect to eV, which shows that the integral over the nucleon density distribution provides a substantial contribution.
Note that the isospin-violating mediator with mass ≲ 10 −11 eV is also constrained by other experiments that use Earth or Sun as the target.To correctly interpret the limits from these experiments, one has to compute the charge of the Earth/Sun by taking into account the charge distributions within these objects.We discuss these constraints in section 5.

Inspiral waveform
The gravitational waveform in the frequency domain (the Fourier representation of the GW strain) can be written as [66,67] where f is the GW frequency, i = (+, ×) denote the "plus" and "cross" polarizations of the GW, A i (f ) is the amplitude, and Ψ(f ) is the phase.In our analysis, we compute the amplitudes and the phase via where A i,G and Ψ G are the amplitudes and phase in the TaylorF2 model in the LALsuite [68], A i,V and Ψ V are the contributions from the new vector mediator.

Amplitude
In the stationary phase approximation (SPA), the waveform amplitude is given by [37,66] A where θ JN is the inclination angle such that cos θ JN = Ĵ • N with J being the angular momentum of the BNS system and N being the line of sight, η ≡ m 1 m 2 /M 2 , D L is the luminosity distance of the BNS system, t s is defined as the stationary point such that ω(t s ) ≡ φ(t s ) = πf with ϕ being the orbital phase of the BNS system, and r(t s ) and ω(t s ) are the r and ω values at time t s , respectively.Because ω = πf at t s , we determine r(t s ) and ω(t s ) by evaluating r(ω) and ω(ω) at ω = πf .We determine r(ω) from Eq. (2.4).In our analysis we treat α as a small parameter and expand r(ω) in powers of α [37]: which is then inserted into Eq.(2.4).We thus find where β ≡ GM m V and v ≡ (GM ω) 1/3 .We compute 1/ ω(ω) via where P V is given Eq.(2.5), E V is given Eq.(2.8), E G is given Eq.(2.9), and P G is given Eq.(2.10).In the last step, we have used the energy conservation relation in Eq. (2.7), and the prime denotes the derivative with respect to ω.We next expand r 2 ω−1/2 to find the leading terms in α and γ, which are then used in Eqs.(3.4-3.5) to obtain A +,V (f ) and A ×,V (f ).We add these two NP terms to the TaylorF2 amplitudes.

Phase
In the SPA, the phase of the GW waveform is given by [37,66] We first compute t(ω) and ϕ(ω) via [67] where ω c , t c and ϕ c are the orbital frequency, time, and the orbital phase at coalescence, respectively.Following the LALsuite's TaylorF2 routines [68], we identify the coalescence as the moment that the orbital frequency approaches infinity, ω c → +∞.We compute We determine t s and ϕ(t s ) by evaluating t(ω) and ϕ(ω) at ω = πf .To determine the NP contribution to the phase, we expand 1/ ω to find the leading terms in α and γ, which are then inserted into Eqs.(3.10-3.11)and further in Eq. (3.9) to obtain the NP phase Ψ V .We add Ψ V to the phase term in the TaylorF2 template.
We note that the TaylorF2 template is obtained at the 3.5 post-Newtonian (PN) order [67].The PN order in Ref. [67] refers to the power of x = v 2 beyond the leading order contributions, where v = (GM ω) 1/3 = (GM πf ) 1/3 is the characteristic velocity of the binary system.For example, at the GW frequency f = 100 Hz, where the LIGO detectors reach their maximum sensitivity [69], one has v ≃ 0.16 for the BNS system in the GW170817 event, and the 3.5PN order expression of E G deviates from the Newtonian expression only by ∼ 2%.To obtain new physics contributions, however, we have avoided using the 3.5PN expansion.This is because the new physics contributions at the 3.5 PN order are so complex that the phase term in Eq. (3.9) cannot be analytically obtained; a numerical treatment on the integrals in Eqs.(3.10-3.11) is not favored, as they could produce significant uncertainties or even hidden errors in our MCMC analysis.Therefore, following Refs.[40,42], we have used Newtonian level expressions to obtain the new physics contributions in this section.For example, the E G formula in Eq. (2.9) and the P G formula in Eq. (2.10) are derived at the Newtonian level (0PN) and then used in Eq. (3.8) and Eq.(3.12) to obtain the new physics contributions to the amplitude and to the phase.
To estimate the difference between our approach and the one where a 3.5PN analysis can be reliably carried out, we compute the new physics contributions that are proportional to α in Eq. (3.12), by using E G and P G at both the 0PN and 3.5PN orders.We find that our approach underestimates the new physics contributions by ∼ 10% at f ∼ 100 Hz for m V ≲ 10 −13 eV.Because of the dominance of the BHSR constraints for 5 × 10 −14 eV ≲ m V ≲ 2 × 10 −11 eV, our GW analysis only probes new parameter space outside this mass range; see section 5 for details.Fig. 3 shows that for m V ≲ 5 × 10 −14 eV, the 3σ upper bound on α is in the range of (0.3 − 0.4).Therefore, the next-leading-order (NLO) corrections of new physics, which are of order α, seem rather significant.Thus, in the parameter space of interest (namely m V ≲ 5 × 10 −14 eV), the uncertainty due to the use of the 0PN expression is smaller than that of neglecting the NLO corrections of new physics, by a factor of ∼ (3 − 4).We leave the analysis of new physics contributions with high-PN orders and/or high-α orders to a future study.

Parameter estimation
To search for the isospin-violating force, we use the Markov-Chain Monte Carlo (MCMC) method to sample the parameter space.We select 40 different mediator masses in the range of [10 −14 , 10 −10 ] eV, equally spaced in the log-scale, as shown in Fig. 3. 5 For each mediator mass, we perform the Bayesian analysis of the GW170817 data [27,52], by using PyCBC [53,54].We evaluate the posterior probability p( ⃗ θ|d, H), where d is the GW170817 data, H is the signal model, and ⃗ θ denotes the parameter space of the signal model H.We obtain the signal model H by modifying the TaylorF2 template [28] to take into account the following new physics effects due to the isospin-violating mediator: (1) the effect on the BNS inspiral dynamics due to the additional Yukawa type force, as given in Eq. (2.2); (2) the new energy loss channel due to the radiation of the light mediator, as given in Eq. (2.5).
In our analysis the parameter space ⃗ θ is spanned by 11 different parameters, which can be grouped into four categories: • Neutron star parameters: the masses of the two neutron stars, m 1 and m 2 , -the aligned spins of the two neutron stars, χ 1z and χ 2z , where the z direction is along the orbital angular momentum of the BNS system, -and the tidal deformabilities of the two neutron stars, Λ 1 and Λ 2 ; • BNS parameters: the inclination angle θ JN , where cos θ JN = Ĵ • N with J being the angular momentum of the BNS system and N being the line of sight; the polarization of BNS ψ, which is the angle between the natural polarization basis of the GW and the reference polarization basis [70]; • The coalescence time t c of the merger; • New physics parameters: α, the relative strength of the Yukawa-force, as given in Eq. (2.3); γ, the parameter for the new radiation channel, as given in Eq. (2.6); In our analysis, we fix the following parameters for the GW170817 event: luminosity distance D L = 40.7 Mpc and the sky location (RA, Dec) = (197.450374• , −23.381495 • ) [28,32,33].Additionally, to better converge the stochastic samplers, we marginalize the coalescence phase ϕ c in estimating the likelihood function [28,71].Table 1 shows the priors of the 11 parameters that are being sampled.
To compute the posterior probability p( ⃗ θ|d, H), we analyze the GW170817 data in the frequency range of f low ≤ f ≤ f ISCO , where f low = 20 Hz [42], and f ISCO is the frequency at the innermost stable circular orbit (ISCO) of the BNS system [37]: Parameter prior Parameter prior The typical prior distributions for the sampled parameters.We use flat (uniform) distributions for the parameters with boundaries that are shown in the brackets.∆t c means we sample the coalescence time around the (LIGO) gps time 1187008882.4 with a range ∆t c .For certain sampler runs, we may adjust the priors of the fifth force parameters {α, γ} since they vary significantly when changing the value of the mediator mass.
where in the last step, we have used m 1 = 1.46 M ⊙ and m 2 = 1.27 M ⊙ [28].The TaylorF2 template no longer accurately describes the BNS dynamics when the frequency exceeds f ISCO [28].

WEP
The WEP experiments test the difference in accelerations between two test masses in an external gravity field [5].The experimental results are usually expressed as limits on the Eötvös parameter, which is defined as the normalized difference of acceleration between two bodies A and B in the same gravity field [72] η The WEP experiments also provide stringent constraints on the fifth-force.For a Yukawa force with a Compton wavelength λ = 1/m V , where m V is the mediator mass, the Eötvös parameter in the external field of the source S is (assuming the fifth-force is extremely small compared to gravity) where r is the distance between the experiment and the source S, and Q i (m i ) denotes the charge (mass) of the object i with i being A, B, and S, respectively.In our notation, the charge of a point mass in the isospin-voilating case is given by Eq. (2.14), namely Q pt = f p Z + f n N , where f p (f n ) is the mediator's coupling to proton (neutron), and Z (N ) is the proton (neutron) number of the point mass.We note that Refs.[10,14] used a different notation: where α is the parameter characterizing the strength, q is the quantum number, and µ is the atomic mass in atomic units such that µ = 12 for carbon-12.For example, in the U (1) B case, q is the baryon number, and α = g 2 /(4πGu 2 ), where g is the gauge coupling, and u ≃ 1.66 × 10 −24 g is the atomic mass units.

MICROSCOPE
The MICROSCOPE experiment, which orbits at the altitude of 710 km, uses the Titanium and Platinum alloys as the two test masses and the Earth as the source of the external field [11].The recent MICROSCOPE constraints are η Ti,Pt = −1.5 +3.8 −3.8 × 10 −15 [13], where Ti and Pt denote the two test masses.
To compute the charges of the two test masses under the isospin mediator, we first note that the sizes of the test masses are small compared to the mediator wavelength of interest, and we thus neglect the effects due to their sizes.Table (1) of [14] provides the B/µ and (B − L)/µ values for Ti and Pt.Because B − L = Z and B = Z + N , we find that for the Ti and Pt alloys: (Z/µ) Ti = 0.40358, (Z/µ) Pt = 0.46061, (N/µ) Ti = 0.59668, and (N/µ) Pt = 0.54044.The charge-to-mass ratio for the Ti and Pt alloys can then be computed.For example, the charge-to-mass ratio for Ti is We next compute the charge of the source (the Earth), denoted as Q E .Following Ref. [14], we consider a simple Earth model, where the Earth consists of a core and a mantle, each with a constant density.The charge-to-mass ratio of the Earth is where m E is the Earth mass, β(x) ≡ x 3 Φ(x), R c (R E ) is the radius of the core (Earth), is the charge of the core (mantle). 6In our analysis, we use R c = 3480 km and R E = 6371 km [73].We compute ) where m c is the mass of the core, and we have assumed that the Earth core (mantle) consists of only Fe (SiO 2 ).By averaging the isotopes on Earth as given in Table 2, we obtain that (Z/µ) Fe ≃ 0.4656, (N/µ) Fe ≃ 0.5356, (Z/µ) SiO 2 ≃ 0.4993, and (N/µ) SiO 2 ≃ 0.5013.We use the Preliminary reference Earth model (PREM) [73] to obtain m c /m E ≃ 32.32%.
The proton (neutron) number of the Earth can be obtained by setting (5.4, 5.5, 5.6).The proton and neutron numbers of the Earth as a function of the mediator mass are shown Fig. 1: For the mediator mass ≳ 10 −14 eV (corresponding to λ ≲ 2 × 10 4 km), the proton and neutron numbers of the Earth start to deviate from the long-wavelength limit.Table 3.The isotopes of hydrogen and helium, the natural abundance (mole-fraction) of the isotopes in the solar system [76], and the atomic mass [75].

LLR experiments
The LLR experiments test both weak and strong equivalence principles, with the measurements of round trip travel times of short laser pulses between observatories on the Earth and retro-reflectors on the Moon [22].LLR provides constraints on the Earth-Moon differential acceleration in the field of the Sun: [22]: η Earth,Moon = −3 +5 −5 × 10 −14 .To analyze the constraints on the isospin-violating mediator, we need to compute the charges of the Earth, the Moon, and the Sun.Because the Earth and the Moon orbit the Sun at a distance of ≃ 1.5 × 10 8 km (1 AU), which is very large compared to the radius of the Sun R ⊙ ≃ 7 × 10 5 km (as well as the radius of the Earth R ⊕ and the radius of the Moon R Moon ), the effects due to their sizes are negligible.The calculation of the charge of the Earth is the same as the MICROSCOPE experiment.
We next compute the charge of the Sun.The Sun is mainly composed of hydrogen and helium, with mass fractions of 73.81% and 24.85% [76], and the abundance of the various isotopes shown in Table 3.For the heavy elements, which are small in the mass fraction of the Sun, we take (Z/µ) = (N/µ) ≃ 0.5.We thus obtain (Z/µ) ⊙ ≃ 0.8632 and (N/µ) ⊙ ≃ 0.1309 for the Sun.
We next compute the charge of the Moon.We adopt the Lunar Primitive Upper Mantle (LPUM) model [77] for the Moon, in which the major compositions include 46.1% SiO 2 , 38.3% MgO, 7.62% FeO, 3.93% Al 2 O 3 , and 3.18% CaO.We assume that the natural abundance of the various isotopes are the same as the Earth, as given in given in Table 3.By averaging the various isotopes for the five major compositions and taking (Z/µ) = (N/µ) ≃ 0.5 for all other minor compositions, we obtain (Z/µ) Moon = 0.4958 and (N/µ) Moon = 0.5048 for the moon.

EW experiment
The EW experiment utilizes a remarkable sensitive torsion balance to test the WEP violation, which has achieved a precision of 10 −13 [9,10].The EW experiment compares different pairs of test bodies, including the beryllium-aluminum pair and the beryllium-titanium pair; it also uses different attractors, including geocenter, the Sun, and the galactic center [10].This enables it to probe the composition-dependent fifth force [6,10].The EW constraints on the baryon number case and the B − L case are given in Refs.[9,10].

Asteroids and Planets
Trajectories of asteroids and planets (hereafter AP) around the Sun offer an ideal testing ground for the laws of physics, where deviations from general relativity predictions can be used to probe new physics.In the presence of a Yukawa type interaction, the equation of motion of the asteroids and planets in the gravitational field of the Sun is given by [55][56][57]: where u = 1/r with r being the distance between these celestial objects, φ is the azimuthal angle of the motion, L denotes the orbital angular momentum per unit mass, λ = 1/m V is the Compton wavelength of the mediator, and α is given by Eq. (2.3) with the two objects being the Sun and the asteroid/planet.The first, second and third terms on the right-hand side of Eq. (5.7) describe the effects of Newtonian physics, GR corrections, and the fifth force, respectively.The data of asteroids [55] [56] and planets [18,57] in the solar system have been used to place constraints on new physics, including the fifth force.In our analysis we adopt the upper bound on α in figure 1 of [56] which combines the results in Refs.[55][56][57].To compute the constraints on the isospin-violating mediator, we calculate the charge of the Sun as in section 5; for the asteroids and planets, which mostly consist of heavy elements, we use (Z/µ) AP = (N/µ) AP ≃ 0.5.Because trajectories of asteroids and planets are sensitive to the fifth force mediator in the wavelength ≃ (0.1, 10) AU, we do not consider the effects due to the finite size of the Sun, asteroids, and planets.For heavy mediator mass ≳ 10 −17 eV, the dominant AP constraints come from Ref. [55], which, however, only provides upper bound on α for m V ≲ 2 × 10 −16 eV.To compare with the GW170917 limits, we extrapolate the AP constraints with an exponential factor of e −a m V [57], and fit the coefficient a using the constraints for m V < 2 × 10 −16 eV.The results is shown in Figs.4-6.

BHSR
The energy and angular momentum of a black hole can be transferred into its surrounding "cloud" that consists of ultralight bosonic fields, through the process known as the BHSR [78].The BHSR relies on only the gravitational interaction, thereby making it applicable to a variety of ultralight bosonic fields.Moreover, because the ultralight bosonic fields can be spontaneously produced, the BHSR does not require a preexisting abundance of the bosonic fields [58].The BHSR constraints on ultralight weakly-coupled spin-1 particles have been analyzed in Ref. [58]: the measurements of rapidly spinning BHs in X-ray binaries exclude vector particles in the mass range of (5 × 10 −14 , 2 × 10 −11 ) eV, and the spin measurements of supermassive BH exclude vector particles in the mass range of (6 × 10 −20 , 2 × 10 −17 ) eV (with a lower confidence level).Note that the BHSR constraints are analyzed under the assumption that the ultralight bosonic field does not possess a significant particle interaction with itself or with other particles.The presence of such interactions can potentially invalidate the BHSR constraints; see e.g., Refs.[79,80] and also Refs.[81][82][83] for BHSR constraints on self-interacting axions.

Results
In this section, we derive constraints on the isospin-violating mediator based on GW170817.There are three NP parameters in the model: the mediator mass m V , the coupling to neutrons f n , and the coupling to protons f p .We carry out MCMC scans for fixed m V values, resulting in two NP parameters for each mediator mass.Because the 2D parameter space defined by f p and f n for each m V is uniquely determined by the 2D space spanned by α and γ, we perform MCMC analysis with α and γ as the free parameters.This approach offers advantages: Firstly, α and γ hold more direct physical significance in the context of GW phenomenology.Secondly, the obtained limits on α and γ can be readily applied to other new light mediators.
We next discuss our method to analyze the GW170817 constraints in the parameter space of the isospin-violating mediators.We discuss the bound related to the α parameter first.In our analysis we choose a fixed ratio of f p /f n and use Eq.(2.3) to compute the f n value from the α value (by using the m 1 and m 2 values) for each MCMC sampling point.Thus we obtain a model point in the new 11D parameter space where α is replaced by f n .
We then magnilize all the other 10 parameters to obtain the marginalized posterior on f n .The bound related to the γ parameter is done in a similar way.Fig. 2 shows the GW170817 constraints on the parameter space spanned by f n and f p /f n for two mediator masses: m V = 10 −14 eV (left panel figure), and m V = 10 −12.1 eV (right panel figure); for both cases, the constraints due to α are stronger than γ.We also compare the GW170817 constraints with the MICROSCOPE constraints: for both panel figures in Fig. 2, the MICROSCOPE constraints [13,15] are several orders of magnitude stronger than the GW170817 constraints, except in the vicinity of f p /f n ≃ −1, where the charge of the Earth becomes zero.We note that the m V = 10 −12.1 eV case is also constrained by black hole superradiance [58].In addition to the default EoS, BSk24, we also compute the upper bound on f n in Fig. 2 by using the α and γ constraints, with an alternative EoS, BSk22 [59].We find that the differences between BSk24 and BSk22 are small.Constraints from the α and γ bounds are shown as the red and blue shaded regions, respectively, with the default EoS, BSk24 [59]; in comparison, we also show the upper bounds with the BSk22 EoS (dashed) [59].Also shown are constraints from the MICROSCOPE experiment (gray) [15].The m V = 10 −12.1 eV case is also constrained by the black hole superradiance constraints [58].
Moreover, we provide in Fig. 3 the bounds on the α and γ parameters to show the detectability of the fifth force in the GW170817 event.To obtain bounds on α, we compute the marginalized posterior on α, by marginalizing over all other 10 parameters in the posteriors from the MCMC sampling.The left panel figure of Fig. 3 shows the 3 σ (99.7%) credible regions on α.Unlike α, which can take either positive or negative value, γ is always positive.Thus, the 3 σ bound on γ only has an upper bound, which is defined such that the total probability in the parameter range from zero to the upper bound is 99.7%.The right panel figure of Fig. 3 shows the 3 σ (99.7%) upper bounds on γ.We do not show in Fig. 3 the constraints on α (γ) for m V ≳ 10 −10.5 (10 −12 ) eV, where the constraints on α (γ) become of order one so that perturbative calculations can no longer be trusted.
We note that the bound on α can be used to place bounds in some other ultralight mediator models.We find that the bound on α in Fig. (3) leads to a limit on f n (with m 1 = 1.46 M ⊙ and m 2 = 1.27 M ⊙ ) that is only ≲ 2% different from that in Fig. 2.However, the bound on γ should be used in caution.We find that the bound on γ in Fig. 3 leads to a limit on f n that is about two orders of magnitudes stronger than in Fig. 2.This is largely due to the fact that γ, computed via Eq.(2.6), is proportional to the difference of the charge-to-mass ratios of the two NSs in the GW170817 event, which, however, have similar masses.Thus, the interpretation of γ is very sensitive to the masses of the NSs., where the dots indicate the medians.We compute the limits, by using the GW170817 data from Ref. [52].
We next consider two special cases in the parameter space of the iso-spin violating model: (1) the baryon number case in which f p = f n ; (2) the case where the Earth charge is zero, namely Q E = 0.

The baryon number case
Here we consider the baryon number case, namely f p = f n .For each mediator mass m V , we compute the constraints on f n with the same method as in Fig. 2. The obtained GW170817 constraints are shown in Fig. 4.
In our analysis we have carefully taken into account the effects due to the finite sizes of the NSs.We find that for heavy mediator mass in the range of m V ≃ (2 − 3) × 10 −11 eV, the NS charge is a factor of ≃ (13 − 34)% larger than the naive calculation in which the NS is treated as a point charge.This then leads to an upper bound on the coupling that is smaller by the same factor.We note that because the effects due to the finite sizes of the NSs increase significantly with the mediator mass, as shown in the left panel figure of the Earth.Because of the Q E = 0 condition, the 3D parameter space (m V , f p , f n ) is now reduced to a 2D parameter space, for which we use (m V , α).Thus, we perform MCMC runs with m V and α as free parameter, but with γ given by7 Fig. 6 shows the GW170817 constraints on the Q E = 0 case.Because the Earth charge is zero, the leading constraints in Fig. 4 from the EW experiment [9] and from the MICROSCOPE experiment [13,15], both of which use Earth as the source of the external fields, are no longer present.The other dominant experimental constraints now include LLR [22], the asteroids and planets (AP) [55][56][57], and the BHSR constraints [58].We find that in the Q E = 0 case, the GW170817 data provide the most stringent constraints in the mediator mass range ≃ (3 × 10 −16 , 5 × 10 −14 ) eV.

Summary
We analyze the constraints on ultralight isospin-violating mediators from the GW170817 event, the first and most robust BNS merger event detected by the LIGO/Virgo collaboration.Although fifth-force experiments, such as MICROSCOPE and EW, usually impose more stringent constraints on ultralight isospin-violating mediators (see e.g., the constraints on the U (1) B mediator, a special case of the isospin-violating mediator, in Fig. 4), we find that the GW170817 event can offer leading constraints on the parameter space where the isospin-violating force is screened by the Earth (namely the charge of the Earth is zero).In such cases, the GW170817 event excludes the existence of ultralight isospin-violating . The 3σ constraints on the coupling f n from the BNS merger event GW170817 (denoted as GW) (blue), for the Q E = 0 case.The blue dots indicate the actual model points that are given in Fig. 3.We extrapolate the GW constraint from m V = 10 −14 eV to the massless limit, as the GW170817 data cannot distinguish the fifth force signals for such a low mass mediator.Other constraints are also shown here as shaded regions: asteroids and planets (denoted as AP) data (green) [55][56][57], lunar laser ranging (denoted as LLR) results (orange) [22], and the black hole superradiance (denoted as BHSR) (gray) [58].Note that the AP limits for m V ≳ 2 × 10 −16 eV are extrapolated from Ref. [55].
To analyze the fifth-force-induced GW signal, we compute the gravitational waveform by taking into account the new physics effects, including the long-range Yukawa force and the dipole radiation.We note that these two new physics effects have different dependences on the charge-to-mass ratios of the NSs: the Yukawa force is proportional to the product of the charge-to-mass ratios, while the dipole radiation is proportional to the quadratic power of their difference.Moreover, we note that the GW170817 constraints given by the dipole radiation of the light mediator are suppressed due to the small charge difference between the NSs; we anticipate that future detection of NS-BH mergers can potentially strengthen these constraints by orders of magnitude.
To correctly interpret the experimental data, we have computed the charges of the NSs by taking into account the effects due to their finite sizes, which are particularly important when the wavelength of the mediator becomes small compared to the NS radius.In such cases, the charge is much larger than in the long-wavelength limit where the NS is treated as a point charge.This enlargement of the NS charge at short wavelength (large mediator mass) offsets somewhat the exponential decline of the constraints on the mass, making the large mediator mass region more promising than naively expected.For the interpretation of other constraints, we have provided a simple analytic expression of the Earth charge, correcting a previous erroneous formula in Ref. [14].
In our analysis we have used perturbation calculations to compute the new physics effects and restricted our attention to the mediator mass range of ≲ 10 −11 eV, beyond which the perturbation calculations begin to fail.We note that the mass range can be extended to larger mediator masses, if one uses other calculations instead of the perturbative ones.
We leave this to our future studies.The NS mass as a function of the central mass-energy density, calculated with the BSk24 EoS [59].The primary mass m 1 = 1.46 M ⊙ and the secondary mass m 2 = 1.27 M ⊙ in the GW170817 event [65] are shown.Right: The proton fraction as a function of the mass of the neutron star, calculated with both the BSk24 EoS (solid, black) and BSk22 EoS (dashed, gray).
To compute ρ(r), we solve simultaneously the TOV equation [61,62], the mass balance equation [84], and the EoS.The TOV equation is where G is the gravitational constant, P is the pressure, r is the radial coordinate, and m(r) is the mass within r.The mass balance equation for a spherically symmetric system is dm(r) dr = 4πr 2 ρ(r). (A.4)

Figure 1 .
Figure1.The proton and neutron charges of the two NSs (denoted as m 1 and m 2 ) in the GW170817 event and of the Earth as a function of the mediator mass.We normalize the charge to the longwavelength limit Q pt (left panel), and to Q u = Q pt Φ(R/λ) which is given in Eq. (2.11) (right panel).The masses of the two NSs in GW170817 are m 1 = 1.46 M ⊙ and m 2 = 1.27 M ⊙ .

Figure 2 .
Figure2.GW170817 constraints (3σ) on the parameter space spanned by f p /f n and f n , for m V = 10 −14 eV (left panel) and m V = 10 −12.1 eV (right panel).Constraints from the α and γ bounds are shown as the red and blue shaded regions, respectively, with the default EoS, BSk24[59]; in comparison, we also show the upper bounds with the BSk22 EoS (dashed)[59].Also shown are constraints from the MICROSCOPE experiment (gray)[15].The m V = 10 −12.1 eV case is also constrained by the black hole superradiance constraints[58].

Figure 3 .
Figure 3.The 3σ bounds on α (left panel) and γ (right panel) for different mediator masses from the BNS merger event GW170817, where the dots indicate the medians.We compute the limits, by using the GW170817 data from Ref.[52].

Figure 5 .
Figure 5.The Q E = 0 curve in the parameter space spanned by f p /f n and λ = 1/m V .

Figure 7 .
Figure 7. Left: The proton fraction as a function of the nucleon number density for both BSk22 and BSk24 EoS, which are adopted from figure30of[59].Right: The nucleon number density as a function of the radial distance for neutron stars with different masses, where only the BSk24 EoS[59] is used.

Figure 8 .
Figure 8. Left: The NS mass as a function of the central mass-energy density, calculated with the BSk24 EoS[59].The primary mass m 1 = 1.46 M ⊙ and the secondary mass m 2 = 1.27 M ⊙ in the GW170817 event[65] are shown.Right: The proton fraction as a function of the mass of the neutron star, calculated with both the BSk24 EoS (solid, black) and BSk22 EoS (dashed, gray).

Table 2 .
[74]isotopes of various elements, the natural abundance (mole-fraction) of the isotopes in the Earth[74], and the atomic mass[75].