The $\nu_{R}$-philic scalar: its loop-induced interactions and Yukawa forces in LIGO observations

Right-handed neutrinos ($\nu_{R}$) are often considered as a portal to new hidden physics. It is tempting to consider a gauge singlet scalar $(\phi)$ that exclusively couples to $\nu_{R}$ via a $\nu_{R}\nu_{R}\phi$ term. Such a $\nu_{R}$-philic scalar does not interact with charged fermions at tree level but loop-induced effective interactions are inevitable, which are systematically investigated in this work. The magnitude of the loop-induced couplings coincidentally meets the current sensitivity of fifth-force searches. In particular, the loop-induced coupling to muons could be tested in the recent LIGO observations of neutron star mergers as there might be a sizable Yukawa force in the binary system mediated by the $\nu_{R}$-philic scalar.


I. INTRODUCTION
Right-handed neutrinos (ν R ) are one of the most intriguing pieces to be added to the Standard Model (SM). Not only can they resolve several problems of the SM including neutrinos masses, dark mater, and baryon asymmetry of the universe, 1 their singlet nature under the SM gauge symmetry also allows for couplings to hidden or dark sectors, a feature known as the neutrino portal to physics beyond the SM.
Among various new physics extensions built on ν R , a gauge singlet scalar φ coupled exclusively to ν R , referred to as the ν R -philic scalar, is arguably the simplest. At tree level, the ν R -philic scalar does not interact directly with normal matter that consists of electrons and quarks, which implies that it might have been well hidden from low-energy laboratory searches. At the one-loop level, there are loop-induced couplings of φ to charged fermions, which are suppressed by neutrino masses (m ν ) in the framework of Type I seesaw [3][4][5][6][7]. The suppression can be understood from that in the zero limit of neutrino masses, which corresponds to vanishing couplings of the SM Higgs to ν R and left-handed neutrinos (ν L ), the ν R sector would be entirely decoupled from the SM content. As we will show, for electrons, the loop-induced effective Yukawa coupling is of the order of where G F is the Fermi constant and m e is the electron mass. Despite the small value of the loop-induced coupling, the magnitude coincides with the sensitivity of current precision tests of gravity. For long-range forces mediated by ultra-light bosons coupled to electrons or quarks, experimental tests of the strong (based on the lunar laser-ranging technology [8]) and weak equivalence principles (e.g., torsion-balance experiments [9,10]) are sensitive to Yukawa/gauge couplings spanning from 10 −20 to 10 −24 . Very recently, gravitational waves from black hole (BH) and neutron star (NS) binary mergers have been detected by the LIGO/VIRGO collaboration [11,12], providing novel methods to test theories of gravity as well as other long-range forces [13][14][15][16][17][18][19][20][21][22][23]. For instance, the process of BH superradiance can be used to exclude a wide range of ultra-light boson masses [14]. The sizable abundance of muons in NS binary systems allows us to probe muonic forces as they could modify the orbital dynamics. It is expected that [23] current and future observations of NS binaries are sensitive to muonic Yukawa/gauge couplings ranging from 10 −18 to 10 −22 which, again, coincidentally covers the theoretical expectation of the loop-induced coupling for muons, G F m µ m ν /(16π 2 ) ∼ 10 −19 .
In light of the frontiers of precision and novel tests of gravity and gravity-like forces, it is important to perform an in-depth study on the loop-induced interactions of the ν R -philic scalar, which is the main goal of this work. We note here that in the seminal work on majorons [24], similar loop-induced interactions have been computed and confronted with experimental limits in the 1980s. While the result in Ref. [24] only applies to a pseudo-scalar boson, in this work we compute loop-induced interactions for a generic scalar and take three lepton flavors into account, with loop calculation details presented. The loop-induced interactions computed in this work could be of importance in phenomenological studies of long-range forces [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42].
The paper is organized as follows. In Sec. II, we briefly review the Type I seesaw extended by a gauge singlet scalar, and derive the tree-level interactions for later use. In Sec. III, we compute the loop-induced interactions of φ with charged fermions, and discuss in detail the UV divergence cancellation, ξ-dependence in the R ξ gauge, and a GIM-like cancellation. The calculation, for simplicity, is first performed assuming only one generation of leptons and then generalized to three flavors in Sec. IV. In Sec. V, we confront the theoretical predictions to experimental limits including searches for long-range forces of normal matter and the LIGO observations of NS events which are sensitive to muonic couplings. We conclude in Sec. VI.

A. Notations
Throughout this paper, Weyl spinors are frequently used in our discussions for simplicity. On the other hand, for Feynman diagram calculations, Dirac or Majorana spinors are more convenient due to a variety of techniques and especially many modern computation packages that have been developed. As both will used in this paper, it is necessary to clarify our notations regarding Weyl spinors versus Dirac/Majorana spinors.
All four-component Dirac/Majorana spinors in this paper are denoted by ψ X with some interpretative subscripts X. Otherwise, they are Weyl spinors. For instance, ν L and R are Weyl spinors of a left-handed neutrino and a right-handed charged lepton, respectively. In contrast to that, ψ is a Dirac spinor of a charged lepton containing both left-and right-handed components.
For Weyl spinors, our notation follows the convention in Ref. [43]. For example, the mass and kinetic terms of ν R are Here and henceforth, the Weyl spinor indices α,α, β will be suppressed. Dirac and Majorana spinors can be built from Weyl spinors. Hence the Dirac spinors of charged leptons and neutrinos can be written as The Majorana spinor of a neutrino mass eigenstate ν i (where i = 1, 2, 3, · · · ) is defined as Note that it is self-conjugate: ψ c i = ψ i . For later convenience, some identities are listed below to convert Weyl spinors into Dirac/Majorana spinors : where P L/R ≡ (1 ∓ γ 5 )/2 and γ µ L ≡ γ µ P L .

B. Lagrangian
We consider the SM extended by several right-handed neutrinos ν R and a singlet scalar φ. In Type I seesaw, the number of ν R needs to be ≥ 2 in order to accommodate the observed neutrino oscillation data. Let us start with one generation of leptons and ignore the flavor structure (for the realistic case including three generations, see Sec., IV). The Lagrangian of ν R and φ reads: The Dirac masses of leptons are generated by where H is the SM Higgs doublet ( H ≡ iσ 2 H * ), L = (ν L , L ) T is a left-handed lepton doublet, and R is a right-handed charged lepton. After electroweak symmetry breaking, H = (0, v) T / √ 2, Eq. (8) leads to the following mass terms: where The Dirac and Majorana mass terms of neutrinos can be formulated as which then can be diagonalized by Here ν 1 and ν 4 are the light and heavy mass eigenstates with their masses determined by The unitary matrix U is parametrized as where c θ ≡ cos θ, s θ ≡ sin θ, and Eq. (14) has been parametrized in such a way that m D , M R , m 1 and m 4 are all positive numbers.

C. Interactions in the mass basis
Since ν L and ν R are not mass eigenstates, we need to reformulate neutrino interactions in the mass basis, i.e., the basis of ν 1 and ν 4 . The two bases are related by Neutrino interactions in the original basis (chiral basis) include gauge interactions and Yukawa interactions, summarized as follows: where g is the gauge coupling of SU (2) L in the SM, c W is the cosine of the Weinberg angle, and H ± is the charged component of H, i.e. the Goldstone boson associated to W ± . Now applying the basis transformation in Eqs. (16) and (17) to Eq. (18), we get Here i and j take either 1 or 4. The couplings g ij Z , g i W , y i ν , y i , y ij R are given by the following matrices or vectors: Eq. (19) can be straightforwardly expressed in terms of Dirac and Majorana spinors according to Eqs. (5) and (6): Note that in the mass basis, φ couples to both heavy and light neutrinos but the coupling of the latter is suppressed by s θ .

III. LOOP-INDUCED INTERACTIONS OF φ WITH CHARGED LEPTONS
As shown in the previous section, at tree level the scalar singlet φ only couples to neutrinos, including light and heavy ones in the mass basis. It does not interact with other fermions directly. In this section, we show that one-loop corrections lead to effective interactions of φ with charged leptons.
From Eq. (19), it is straightforward to check that at the one-loop level, in the unitarity gauge (which means Goldstone boson interactions can be ignored), there are only two possible diagrams that can connect φ to charged leptons or quarks, as shown in Fig. 1. The second diagram involving the Z boson actually leads to a pseudo-scalar coupling (see calculations later on). In unpolarized matter, pseudo-scalar interactions cannot cause significant long-range forces [44,45] because the Yukawa potential between two fermions are spin dependent. When taking an average over the spins,  (35), and the right diagram leads to a pseudo-scalar coupling (with γ 5 ), the effect of which however is suppressed in unpolarized matter. The diagrams are presented in the mass basis (ν i and ν j are mass eigenstates). For an equivalent description in the chiral basis, see Fig. 2. Figure 2. The W ± -mediated loop diagram in the chiral basis, which is equivalent to the left diagram in Fig. 1 in the mass basis. It shows explicitly how chirality changes in the process. Since in the chiral basis W ± only couples to left-handed leptons and φ only to ν R , we need two mass insertions of m D to connect ν L and ν R . Other two mass insertions, M R and m , are also necessary due to additional requirements-see discussions in the text.
the effect of pseudo-scalar interactions vanishes. Therefore, we will focus our discussions on the first diagram where the external fermion lines have to be charged leptons. The diagrams in Fig. 1 are in the mass basis which is technically convenient for evaluation. Nonetheless it is illuminating to show Fig. 2, another diagram in the chiral basis which explicitly shows how chirality changes in the process. The physical results should be basis independent. Fig. 2 follows directly from Eq. (18), which suggests that φ only couples to ν R while W ± interacts with ν L . Therefore, two Dirac mass insertions (m D ν L ν R and m D ν † L ν † R ) are necessarily introduced to connect ν R and ν L , or ν † R and ν † L . Note that the two W ± vertices have to be conjugate to each other, which implies that from the W ± side, a pair of ν L and ν † L is provided. On the other hand, the Yukawa vertex couples φ to two ν R 's rather than a pair of ν R and ν † R . So a Majorana mass insertion is required to flip the lepton number and convert one of them to ν † R . The direction of lepton-number flow in this diagram are represented by the arrows. Note that according to the conventions in Sec. II A, ν L and ν R have opposite lepton numbers. So for ν R ν R φ, the arrow of ν R should be outgoing. In contrast to that, the arrow of ν L in the W − µ † L σ µ ν L vertex goes inwardly. Finally, there should be a mass insertion of m L R in one of the external legs because it is impossible to write down an effective operator that consists of φ and two L 's -the operator φ L L is not allowed due to electric charge conservation.
The chirality analysis in Fig. 2 indicates that the diagram would be proportional to m D 2 M R m if all these masses are sufficiently small. This will be verified in the loop calculation presented below. Now let us compute the loop diagrams explicitly. Using the Dirac/Majorana spinor representation in Eq. (22), we can write down the amplitudes of the two diagrams in Fig. 1: where (i) 3 comes from three vertices; p 1 and p 2 are the momenta of the upper and lower external fermion lines; p i and p j are the momenta of ν i and ν j ; q = p 2 − p 1 = p j − p i ; and k is the momentum of W propagator. The symbol ∆ denotes propagators. For Majorana spinors in the mass basis, their propagators have the same form as Dirac propagators: The propagators of W ± and Z are gauge dependent. Most generally, in R ξ gauges, they are: The unitarity gauge corresponds to ξ → ∞. Except for the unitarity gauge, other gauges with finite ξ, e.g., the Feynman-'t Hooft gauge (ξ = 1) and the Lorentz gauge (ξ = 0), require the inclusion of Goldstone boson diagrams. The unitarity gauge, albeit involving fewer diagrams by virtue of infinitely large masses of Goldstone boson propagators, has a disadvantage in that the cancellation of UV divergences is less obvious-see discussions in Sec. III A. Nonetheless, it is straightforward to compute iM W and iM Z for general values of ξ. First, we inspect the iM Z amplitude. The loop integral of the trace part gives rise to a quantity proportional to q ν :ˆd which can be expected from Lorentz invariance, explained as follows. On the left-hand side of Eq. (28) there are only two independent momenta p j = p i + q and p i . After p i being integrated out, the only quantity that can carry a Lorentz index is q so the result is proportional to q ν . Now plugging this in Eq. (24), we can immediately get a γ 5 sandwiched between u(p 2 ) and u(p 1 ): Therefore, due to the aforementioned suppression of pseudo-scalar effect in unpolarized matter, we will not investigate more into the iM Z amplitude. The iM W amplitude can be computed by splitting the W ± propagator in Eq. (26) to two parts: where the first part does not contain ξ and the second part is important for cancellation of UV divergences. Note that when computing Eq. (23), because of the chiral projectors in y ji R P L +y ji * R P R , the product of Dirac matrices gives It implies that if m i → 0 and m j → 0, the result would be zero, which agrees with our analysis in the chiral basis.
With the above details being noted, we compute 2 Eq. (23) in the soft scattering limit (q → 0) and in the approximation of m m W : where and F 1 and F 2 correspond to the contributions of the first and second parts of the W ± propagator in Eq. (30), respectively. The explicit expressions F 1 and F 2 read as follows: Here we have used dimensional regularization which means the integrals are computed in a d = 4 − 2 dimensional spacetime. And the generalization of integration measure´d introduces a dimensional constant µ which, together with 1/ , should be canceled out in physical results.
We have verified that the above expressions are symmetric under i ↔ j: In addition, though m 2 i − m 2 j appears in some of the denominators, it does not cause additional divergences when m i → m j : lim m j →m i To obtain the final result of iM W , one needs both Eqs. (34)- (35) and Eqs. (37)- (38) to sum over i and j as it involves cases of i = j and i = j. The latter requires the above limit.
There are several noteworthy cancellations in the summation which we would like to discuss as follows.

A. Cancellation of UV divergences
As can be seen from Eqs. (34) and (35), both F 1 and F 2 contain UV divergences 1/ in their first terms. When combined together in Eq. (32), there is obviously a cancellation between the two divergences: Further cancellations of log ξ will be discussed in the next subsection.
Here we would like to address a subtlety concerning UV divergences in the unitarity gauge. If we had naively taken the ξ → ∞ limit of Eq. (26) at the beginning of the above calculations, we would get a divergent result because which corresponds to exactly the F 1 contribution according to Eq. (30). And our calculation has shown that the F 1 contribution itself is UV divergent. We also know that the divergence is actually canceled out by the F 2 contribution, which however would vanish if ξ → ∞ had been taken in the naive way. That implies that taking ξ → ∞ should be after the loop integration. Actually from the second term of (30), one can see that when the loop integral contains , the ξ → ∞ limit does not commute with k → ∞ in the integral. Taking ξ → ∞ after the integration can make the large momentum contribution with k 2 > ξm 2 W be included, which is crucial for the UV cancellation.
In other gauges, it is more straightforward to see that iM W is finite. Taking the Feynman-'t Hooft gauge for example, when it is applied to Eq. (23), using Eq. (31), the loop integral becomeŝ where λ j ≡ m i y ji * R and λ i ≡ y ji R m j . Is is now obvious to see that the loop integral converges because the integrand is proportional to k −5 as k → ∞.

B. Cancellation of ξ dependence in R ξ gauges
The F 1 contribution is ξ independent. So we are only concerned with F 2 . Let us make a series expansion of F 2 with respect to ξ −1 : Usually in R ξ gauges, it is expected that the ξ dependence of a W ± diagram is canceled by a similar diagram with W ± replaced by its Goldstone boson H ± . In our case, it would be the diagram in Fig. 3. However, a straightforward calculation shows that the amplitude of this diagram is which is impossible to cancel the log ξ term in Eq. (43) when ξ increases to sufficiently large values. Figure 3. The Goldstone boson diagram that complements the W ± diagram to cancel the ξ dependence in R ξ gauges.
This problem is essentially related to the completeness of the model. For an arbitrary matrix of y ij R , indeed the result would be gauge dependent and the log ξ term remains for each case of (i, j) = (1, 1), (1, 4), (4, 1), and (4, 4). However, in Sec. II we have shown that the elements in y ij R are correlated by active-sterile neutrino mixing-see Eq. (20). Besides, g i W also depends on the mixing-see Eq. (21). As a consequence, when summing up the contributions of both light and heavy neutrinos in Eq. (32), the log ξ term cancel out because That implies that the log ξ term can be safely ignored when computing the full amplitude. Actually if we inspect G ij in the chiral basis, the cancellation is more manifest. From Eq. (33), we can express i, j G ij in the matrix form: where w and Y R are the vector and matrix of g i W and y ij R in Eqs. (21) and (20) respectively. They are transformed from the chiral basis by: Therefore, in the chiral basis, we have Here M 2ν is the neutrino mass matrix in Eq. (11). It shows that the vanishing product of w and Y R (more specifically, w † Y R = 0 and Y † R w = 0) leads to i, j G ij = 0, which is due to the absence of W ± -ν R and φ-ν L couplings.

C. GIM-like cancellation
If all the neutrino masses (including heavy ones) are much smaller than m W , when summing over i and j in Eq. (32), the leading-order contribution vanishes in a way similar to the Glashow-Iliopoulos-Maiani (GIM) mechanism [47]. At the next-to-leading order (NLO), a nonzero result can be obtained. But in case of zero mass splitting of neutrinos, the NLO contribution would vanish again. This is also similar to the GIM cancellation, where if u and c quarks were of equal mass, the K 0 → K 0 amplitude would be zero.
Let us compute iM W in the unitarity gauge. The previous discussion in Sec. III B concludes that the log ξ term can be safely ignored in this complete model. Hence we define which is finite. Then iM W in the unitarity gauge can be computed by: Now if we assume m W m i and m j , F 12 can be expanded as follows: As we have shown i, j G ij in Sec. III B, the constant term in F 12 does not contribute to Eq. (49), which means that the leading-order contribution to iM W vanishes. Only the second or higher-order terms in Eq. (50), which is further suppressed by m 2 i,j /m 2 W , can lead to nonzero iM W . Plugging in the second term in Eq. (50) to Eq. (49) and using the explicit form of G ij in Eq. (33), we obtain Note that the expression in the square bracket is antisymmetric under the interchange of m 1 and m 4 . Therefore if the mass splitting m 1 − m 4 is zero, the NLO contribution vanishes as well, similar to the GIM cancellation. In Type I seesaw, the scale of heavy neutrino masses is often assumed to be much higher than the electroweak scale. Hence a more likely scenario is m 4 m W m 1 . For such hierarchy, there is no GIM-like cancellation, as we shall show below.
First, we need to expand F 12 in other regimes. If the diagram contains heavy neutrinos running in the loop, we expand it with respect to m W : If the diagram contains one light and one heavy neutrinos, we have For m j m W m i , the result can be obtained by an interchange of i and j in Eq. (53). Combining the results in Eqs. (50), (52), and (53), we sum over i and j in Eq. (49), which gives: where m ≡ m 2 1 + m 2 4 . Now taking G F = √ 2g 2 /(8m 2 W ), ms 2 θ = m 1 , and m m W , we obtain with It implies that the loop diagram generates the following effective interaction where the effective coupling y φ , given in Eq. (55), is suppressed by the neutrino mass m ν and the charged lepton mass m .

IV. GENERALIZATION TO THREE FLAVORS
So far we have only considered leptons of a single flavor for which we have computed the loopinduced coupling y φ , as given in Eq. (55). Now we would like to generalize it to the realistic scenario with three flavors.
Assuming there are three generations of ν L and ν R , we can express the neutrino mass terms in a similar way to Eq. (11) except that now the mass matrix is interpreted as a 6 × 6 matrix: where m D and M R are 3 × 3 Dirac and Majorana mass matrices respectively. In principle, the number of right-handed neutrinos does not have to be three. It can be two or more. But to make it concrete, let us concentrate on the case with three ν L plus three ν R . The neutrino mass terms and Yukawa terms of ν R are formulated as: where Y 0 R is the 3 × 3 Yukawa coupling matrix of φ to ν R . The 6 × 6 symmetric mass matrix can be diagonalized by a 6 × 6 unitary matrix U : Although in general analytical diagonalization of the 6 × 6 mass matrix may be not difficult or impossible, numerically solutions can always be obtained. Then we transform neutrinos to the mass basis: where both ν L and ν R are 3 × 1 vectors. First, we convert the gauge interaction g √ 2 W − µ † L σ µ ν L from the flavor basis to the mass basis: This generalizes g i W in Eq. (21) from a 1 × 2 vector to a 3 × 6 matrix: Next, we perform a similar transformation for the Yukawa interactions ν R , i.e., the second term in Eq. (58): where Y R is a 6 × 6 generalization of the y ij R matrix in Eq. (20). The generalization of G ij is quite straightforward. By replacing g i W and y ij R in Eq. (33) with g i W and Y ij R , we get: As for F 1 and F 2 , the expressions in Eqs. (34) and (35) can be used directly except that now i and j run from 1 to 6 instead of 1 and 4. With the generalized G ij , F 1 and F 2 , the loop-induced effective Yukawa coupling is where M d has been defined in Eq. (59) and F 12 takes the same definition as Eq. (48). Note that any constant terms in F 12 can be ignored because Eq. (64) is our final result for the most general 3ν L + 3ν R scenario. Although its dependence on the PMNS matrix and light neutrino masses is not manifest, each quantity in Eq. (64) can be readily evaluated using numerical methods.
Below we would like to discuss a special case in which Eq. (64) can be further simplified and expressed in terms of the PMNS matrix and light neutrino masses. Assuming that m D , M R and Y 0 R can be simultaneously diagonalized 3 as follows: we can transform neutrinos to an intermediate basis where In this basis, each 3 × 3 sub-blocks of the neutrino mass matrix are diagonal: And the Yukawa coupling matrix in this basis is diagonal: Now the mass matrix on the right-hand side of Eq. (69) can be diagonalized by: where U 14 , U 25 and U 36 are generalizations of the matrix in Eq. (14) to perform small rotations between ν 1 and ν 4 , ν 2 and ν 5 , and ν 3 and ν 6 respectively. Combining the three matrices together, we get the following explicit form of U : where (c i , s i ) ≡ (cos θ i , sin θ i ) and θ i = arctan m 1 /m i+3 .
Therefore, the U matrix which is defined to diagonalized M ν should be And Y R , defined in Eq. (63) as the Yukawa coupling matrix in the mass basis, now reads: Only diagonal elements of each 3 × 3 sub-block of Y R are nonzero: This greatly simplifies the calculation in Eq. (64). It implies that for F 12 (m i , m j ) only i = j and i = j ± 3 need to be computed: where f (m i ) ≡ 12 log m i m W and elements denoted by "." are irrelevant to our calculation. In Eq. (76) we have assumed m 4,5,6 m W m 1,2,3 , which also implies that s 1,2,3 1. Now supplying all the matrices required by Eq. (64), we obtain Note that U L is approximately identical to the PMNS matrix: if we ignore the small difference caused by U . On the other hand, the light neutrino mass matrix can be diagonalized by U T PMNS M ν U PMNS = diag(m 1 , m 2 , m 3 ). Therefore, if we define Eq. (77) cab be reformulated as Note that Eq. (80) is derived under the assumption that m D , M R and Y 0 R can be simultaneously diagonalized. For more general cases, we refer the readers to Eq. (64) to compute the effective couplings.

V. PHENOMENOLOGY
The loop-induced interaction of φ with electrons leads to a Yukawa potential between two objects containing N 1 and N 2 electrons, The effective Yukawa coupling y φee is of order G F m e m ν /(16π 2 ) ∼ O(10 −21 ), which reaches the current sensitivity of long-range force searches. If we replace electrons with muons, the effective coupling is generally two orders of magnitude larger because m µ /m e ≈ 200. The muonic long-range force can be tested in binary systems of neutron stars (NS) which contain O(0.1 ∼ 1%) muons of the total mass [51]. In particular, the recent gravitational wave observation of NS binary mergers by the LIGO collaboration [11,12] are able to test the muonic force with unprecedented sensitivity. As indicated by Eqs. (80) and (79), the value of y φ depends on neutrino masses and the Yukawa couplings of φ to ν R . Since there are many free parameters in diag(Y R1 , Y R2 , Y R3 ) and M ν (where Majorana phases, the Dirac CP phase, the lightest neutrino mass are still unknown), we would like to simply parametrize it as follows: where Y ν is not suppressed or vice versa, analogous to the well-known fact that m ee for neutrinoless double beta decay can vanish in the normal mass ordering.
Next, we shall confront the theoretical predictions with experimental limits, as shown in Figs. 4 and 5 for y φee and y φµµ respectively.
For y φee , current limits come from long-range force searches of normal matter, which have long been investigated in precision tests of gravity, including tests of the inverse-square law (ISL) and tests of the weak equivalence principle (WEP). The Yukawa force mediated by φ can affect the former by contributing an exponential term to the total force and affect the latter due to its leptophilic coupling, which causes differential free-fall accelerations for different materials. So far, the Eöt-Wash torsion-balance experiment has performed WEP tests with the highest precision [9,10], leading to the most stringent constraint on y φee in the regime of very small m φ . In addition, the lunar laser-ranging (LLR) technology which is able to measure the varying distance between the moon and the earth to high precision using laser pulses is also sensitive to new long-range forces [8]. These two bounds, reviewed in Ref. [52], are presented in Fig. 4 and overlap with the theoretically most favored region (red lines).
In Fig. 4 (also Fig. 5), we add hatched bands to represent the constraint from black hole superradiance [14], which is independent of the Yukawa couplings because the effect is cause by φ coupling to the spacetime.
For y φµµ , the aforementioned laboratory constraints do not apply since normal matter does not contain muons. Neutron stars (NS), however, can be a powerful probe of muonic forces due to a significant abundance muons in them, which is expected when the Fermi energy exceeds the muon mass. According to the calculations in Refs. [23,51], the number density of muons is typically of O(1 ∼ 10)% of the total number density, which is lower than but still comparable to the electron number density 4 .
In fact, since the electron and muon number densities are of the same order of magnitude while y φ ∝ m , for NS binaries we have where F (µ) and F (e) are the forces caused by muons and electrons respectively. The recent observations of NS-NS and NS-BH mergers by the LIGO collaboration provide very promising data to probe the muonic force in this model. For a NS-NS merger, the effect of φ is twofold [18]. First, the attractive force affects the motion of the binary system in a classical way, i.e., modifying the Kepler's law when r ∼ m −1 φ . Second, since φ is a ultra-light boson, there is radiation The muonic force could be probed in binary systems of neutron stars (NS) due to the considerable abundance of muons. The blue and green curves represent current sensitivity of the LIGO observations of GW170817 (NS-NS merger) and GW190814 (NS-BH merger) events, respectively. Solid (dashed) curves take conservative (optimistic) estimates of the muon abundance [23]. In addition, precision measurements of binary pulsar systems are also sensitive to the muonic force (orange curves) [23].
of φ from the rotating dipole, which causes extra energy loss. As for a NS-BH merger, only the effect of φ radiation is relevant. An in-depth analysis of the sensitivity to muonic forces based on the recent two events GW170817 (NS-NS merger) and GW190814 (NS-BH merger) has been performed in Ref. [23]. Their results have been incorporated in Fig. 5, where solid (dashed) curves are derived using a conservative (optimistic) estimate of the muon abundance. For GW170817, the sensitivity curves of the two effects are evaluated and presented separately. The one sensitive to m φ < 10 −15 eV is derived from the second effect (φ radiation), and the other comes from the first effect. In addition to binary mergers, precision measurements of binary pulsars can also be sensitive to muonic forces [23,76].
As shown in Fig. 5, the LIGO curves cross the red lines of Y (µ) R = 10 −3 ∼ 1, which implies that the loop-induced muonic force in this model could be probed in the theoretically most favored regime. Future experiments such as the Einstein Telescope 5 and Cosmic Explorer [77] can substantially improve the sensitivity to muonic forces and thus have great potential of probing this scenario.

VI. CONCLUSIONS AND DISCUSSIONS
The ν R -philic scalar model naturally gives rise to extremely small couplings of charged leptons to a long-range force mediator via loop-level processes. The small values of the loop-induced couplings coincidentally meet the current sensitivity of long-range force searches in laboratories and in astrophysical observations such as the recent detection of GW from NS mergers by LIGO, as we have shown in Figs. 4 and 5. In this model, loop-induced couplings to quarks also exist, due to the Z-mediated diagram in Fig. 1. However, our calculation shows that only pseudo-scalar couplings are generated in this case, the effect of which is suppressed in unpolarized matter.
Our loop calculation result for the most general three-flavor case is given by Eq. (64) which, though involving diagonalization of the full 6×6 mass matrix, can be numerically evaluated. For the special case where m D , M R and Y 0 R can be simultaneously diagonalized, the result can be further simplified to Eq. (80), where the dependence on the PMNS matrix is manifestly extracted.
Our results can also be used to obtain loop-induced interactions for other similar models that contain the diagrams in Fig. 1, via properly replacements of the couplings in vertices and masses in propagators. However, one caveat should be noted here that incomplete models where the treelevel couplings of φ to light neutrino mass eigenstates are not governed by the active-sterile neutrino mixing would lead to gauge dependent results.