Sub-GeV $U(1)_{R}$ gauge boson to address the proton radius discrepancy

We propose a Standard Model extension by a $U(1)_{R}$ gauge symmetry where only right-handed chiral fermions can carry a non-trivial charge. Here we show that the simplest anomaly-free solution to accommodate the proton charge radius discrepancy takes right-handed muons $\mu_R$ and first generation quarks, $u_R$ and $d_R$. Consistency with the latest muon's $(g-2)$ measurements is achieved through an extra light scalar, which itself must lie in the tens of MeV mass range to be viable.


I. INTRODUCTION
The Standard Model (SM) of particle physics successfully describes three of the four known fundamental interactions of nature: electromagnetic, weak, and strong. The tremendous accuracy of the SM has been corroborated by contemporary collider experiments at high energies such as the Large Hadron Collider (LHC). Nevertheless, this and other experiments have found several subtle deviations from the SM predictions that remain to be accounted for. For some of these anomalies, if they persist and get confirmed, the presence of new physics is required.
Some of those deviations have received particular attention during the last decade with revitalized interest provided by new measurements. It is peculiar that they all somehow involve the muon, suggesting that subtle effects in its interactions might a be key to new physics. Examples are the B-meson anomalies, which are deviations from the SM predictions in angular observables and the ratios R K(K * ) = Br(B → K ( * ) µ − µ + )/Br(B → K ( * ) e − e + ) that might point to violation of lepton-flavor universality. These have captured a lot of theoretical interest with plenty of effort channeled towards setups that extend the SM with additional scalars and/or gauge bosons coupled to quarks and leptons in a non-universal manner [1]. As of today, the deviations persist at a level above 3σ, based on the latest results from the LHCb collaboration at the LHC [2,3]. Another example of anomalies is the latest measurements of the muon's (g − 2) [4], for which attempts to accommodate it within SM extensions containing muonic mediators already comprises a vast literature.
Our focus is in yet another instance of tentative new physics involving the muon. For decades, the study of the structure of the hadrons found in nuclei -protons and neutrons-has improved as a result of increasingly refined techniques such as lepton-nucleon scattering and the spectroscopy of hydrogen-like atoms. Among the long-standing studied observables there is the (root-mean-square) proton charge radius, r p . Despite the good agreement between measurements of r p by experiments involving electrons prior to 2010, in that year the CREMA collaboration found a highly discrepant r p value through spectroscopy of hydrogen-like atoms where the muon replaces the electron (muonic hydrogen) [5]. The muonic result improved in 2013 [6] and continued to be discrepant with the electronic-measurement averages compiled and published by the CODATA committee in 2014 [7] r (e) p = 0.8751 ± 0.0061 fm r (µ) p = 0.84087 ± 0.00039 fm . (1) Fast forward to 2021, several experiments with electrons and form-factor analyses of the data have found a degree of consistency with the muonic r p value. In particular the CODATA 2018 average [8] and posterior reports [9] tend to favor the muonic r p value. Yet, since a few other recent experimental collaborations [10,11] have reported r p values that persistently agree with the CODATA 2014 r p value, the proton charge radius puzzle is certainly ameliorated but not conclusively gone 1 [10,11]. Similarly, given that the value found by Grinin et. al. [14] is offset the muonic value by about two standard deviations, the puzzle is kept alive.
With the ongoing progress in higher order theoretical corrections in scattering and spectroscopy and the projected sensitivity of near-future experiments, the community has also entertained the possibility that these r p discrepancies are a hint of new physics consisting on muonphilic interactions of the proton [15][16][17][18]. This is the approach taken in this paper: we follow the hypothesis that the ∆r 2 p discrepancy between the two values in Eq. (1) will persist, and is addressed by new degrees of freedom of a complete, anomaly-free gauge U (1) model. Our setup is also consistent with other relevant constraints, including LHC bounds on new scalars and those from the latest muon (g − 2) results by Fermilab [4,19] and the Lattice QCD group [20].
Our work is organized as follows: the new gauge model and its field content is described in in Sec. II, with the proton charge radius discrepancy reviewed and addresssed in Sec. III. Contributions to the muon's (g − 2) are discussed in Sec. IV, where we also study numerically the parameter space where r p can be accomodated while consistent with collider limits on scalars and the muon's (g − 2). Closing remarks follow afterwards.

II. MODEL SELECTION
In this scenario the SM gauge group is enlarged by an Abelian gauge U (1) R symmetry (and its corresponding In addition, the model contains the following U (1) R -charged scalars, needed to break U (1) R spontaneously and give a mass to the Z : a complex SM singlet s, and a SU (2) L doublet Φ 2 (same hypercharge as the SM Higgs field, denoted as Φ 1 ) which will also be responsible for generating the masses of the muon, up, and down quarks. As it is shown in the next section, due to its relative size to the vacuum expectation value (vev) of Φ 2 , the vev of s generates most of the Z mass while accomodating the charge radius discrepancy. Meanwhile, a range of sizes for the Φ 2 vev is discussed at the end of this section. The rest of the SM particles are uncharged under the new symmetry, including the Higgs-like doublet.
The fact that SM fermions of different generations carry distinct U (1) R quantum numbers implies that gauge anomalies do not cancel automatically per generation, thus extra (chiral) fermion content is needed to satisfy the anomaly-free conditions. For this purpose one SM singlet RH neutrino ν R is introduced into the model (with a yet unspecified charge q ν = 0 that does the job when exactly one quark and one lepton generations are charged 3 ).
Given the introduction of the SM singlet ν R , there are 27 different anomaly-free scenarios that one could generate.
The diagram in Fig. 1 exemplifies the 9 gauge anomaly-free combinations involving the RH u-quark, with all three RH d-type quarks and RH charged leptons. The other 18 solutions involve c R and t R . Seeking to establish a connection between nucleons and muons that tackles the proton radius puzzle, the valence RH-quarks of the first generation and the second generation lepton are the ones given a non-trivial charge under U (1) R . This choice is illustrated by the solid red line in Fig. 1. For reference, Table I  QL (1,2,3) uR We choose to work out the q ν = +1 solution for the sake of generality. In this scenario, a Dirac neutrino mass is generated through the Yukawa term ∼ L L ν RΦ2 , whose smallness might be associated to the small induced vev of the additional SU (2) scalar doublet as in the type-II seesaw for Dirac neutrinos [23]. In order to account for all neutrino oscillation observables the model can be further extended with extra RHNs neutral or not under U (1) R . The actual construction of the neutrino sector is out of the scope of this paper.

A. New interactions
The relevant operators and scalar potential are listed next. The fermion interactions with the Z are given by Then, the overall right-and left-handed couplings C R,L to the muon become C R = −g and C L = 0 for q ν = 1. We normalize the vector (V ) and axial vector (A) couplings as C V,A = 1 2 (±C L + C R ) . The gauge-invariant Yukawas of the SM fermions are with i = 1, 2, 3 denoting the fermion generations and where Regarding the generation of the muon mass, we consider the Yukawa matrix for charged lepton as diagonal, so that, A natural muon Yukawa occurs if one selects v Φ2 ≈ m µ within one order of magnitude. Then 0.01 v Φ2 /GeV 1 is equivalent to 0.1 y 2µ 10, with the possibility of saturating y 2µ = 4π.
The sizability of v Φ2 (or equivalently y 2µ ) will prove to be important in Sec. IV to accomodate the (g − 2) µ bounds.
The scalar potential of our Two-Higgs doublet model (2HDM) plus complex singlet reads where the U (1) R charge of the scalar singlet is taken as q s = −1/2. For this charge choice, the allowed four-scalar operator with the Higgs doublets is s 2 Φ † 1 Φ 2 , with dimensionless coefficient κ as in Ref. [23]. An alternative, not pursued here, is choosing q s = −1 under which the allowed four-scalar operator would be sΦ † 1 Φ 2 with a dimensionful mass parameter as coefficient.

III. THE CHARGE RADIUS DISCREPANCY
In our model the contributions to the proton charge radius come from the exchange of the mediators Z , Φ 2 and s between muons and protons. We later discuss why the scalar contributions are subleading/negligible. Furthermore, since left-handed chiral neutrinos do not interact with the nucleus through the Z , they automatically avoid neutrinonucleon scattering constraints [16].
The energy splitting in muonic hydrogen by exchange of the vectorial (V ) component of a spin-1 mediator is given by [15] with the analog ∆E (S) by CP-even scalar exchange 4 sharing the same expression. Above, the C µ(p) V are the vectorial muon (proton) couplings to the gauge boson, m red ≡ m p m µ /(m µ + m p ) is the muon-proton reduced mass, m V is the mediator mass, and α em is the QED constant. The Lamb energy shift between levels 2S 1/2 and 2P 1/2 depends on the radius r (in fm) through the following relation [12], So we are using Eq. (7) to fit the proton charge radius discrepancy shown in Eq. (1). This will correspond to an energy shift of ∆E ≈ 306.4 µeV.
For purposes of simplicity, we take the Z coupling to the u and d quarks as an estimate of C p V , the vectorial Z coupling to the proton 5 . Then from Eq. The square of this value is the diagonal entry of the new gauge boson in the full 3 × 3 interaction-basis mass matrix of the neutral vector bosons after the Weinberg rotation. In the limit of vanishing kinetic mixing, this value is a good approximation for m Z because the off-diagonal mass mixing of the U (1) Y and U (1) R gauge bosons goes as and this is much smaller than |m 2 Z − m 2 Z | (see for example the paramametrization of Ref. [24]). On the (m Z , g ) plane in Fig. 2, a solid black curve satisfying |∆E (V ) | ≈ 306.4 µeV is shown, together with a few contours of the corresponding v s singlet vev (magenta dashed), which we will assume dominates the Z mass, i.e m Z ≈ g v s /2. These dashed contours indicate that 5 v s /GeV 10 works well for ∆r 2 p , except when m Z gets as low as 1 MeV 6 .
With the previous information it is suggested that the proton radius puzzle can be addressed in this model thanks only to existence of a light Z that interacts only with one family of RH-quarks and -leptons. However, the assumption on the vevs (v s ∼ few GeVs v Φ2 ∼ m µ ) has an impact on the scalar mass spectrum which is constrained by current experimental data. Moreover, given the features of the model, we will also consider the latest muon (g − 2) results, which provide important parameter space restrictions and predictions.

IV. MODEL CONSTRAINTS AND RESULTS
In what follows we will describe the experimental limits applied to our setup. As discussed earlier, from Fig. 2, and under our choice q s = −1/2, a vev v s ≈ 6.5 GeV fits ∆r 2 p for m Z above 1 MeV. This and other parameters choices will be collected in Table II, where subindices will be dropped for the muon Yukawa, y ≡ y 2µ , to ease notation. Now we proceed to discuss the scalar sector and its experimental restrictions.

A. Limits on the scalar mass spectrum
The singly-charged Higgs from the Φ 2 doublet is subject to severe LHC constraints. In the context of a muonphilic 2HDM, a bound m H + 640 GeV has been reported [25]. We choose two representative benchmarks, one with the m H3 , m A , m H + sitting right at this bound, and another, more conservative one where they sit at 1 TeV. In order to have a grip on the s-like and ϕ 2 -like scalar masses, we look closer at the mass matrices in a regime of small mixing.
For a O(1), sizable κ coupling and under the vev hierarchy v Φ1 v s v Φ2 , the M 2 R mixing matrix in (A3) can always be brought to an approximate block-diagonal form under a suitable choice of quartics. Above, ε, ε entries are no larger than the (1, 1) entry, but much smaller than the (3, 3) and (2, 3) entries. Then the diagonalization of the lower-right 2 × 2 block in (σ R , ϕ 0 2R ) approximates the light s-like and ϕ 2 -like masses, and it does it in terms of λ s , κ , v Φ2 only (recall that v 2 Φ1 = (246 GeV) 2 − v 2 Φ2 , and that v s is fixed by ∆r 2 p ). The approximate eigenvalues of this block are At next-to-leading order in v Φ2 , the light mass is whereas the heavy one is well approximated by the diagonal entry in ϕ 2R at leading order in v Φ2 The charged scalar mass is also approximated by (11), and this is also true for the pseudoscalar A too. Therefore, the H 3 , A, H + in this regime are highly degenerate. Given these masses, m heavy is freed from the 640 GeV bound by LHC on charged scalars if Regarding the s-like singlet scalar, its mass m light must be rather light in order to help the total ∆a µ reach the current limit. Since its coupling to muons is supressed by the small s-Φ 2 mixing, y must be large enough to compensate. At light mass values, however, a tachyonic m light must be avoided. This is achieved if Contours of the light singlet (solid green) and charged scalar (dashed orange) masses are shown in the panels of We list the mass spectrum at two benchmarks in Table II, and mark them with red crosses in the panels of Fig. 3.
Notice that the charge scalar mass m H + sits at either 640 GeV or 1 TeV, with m A ≈ m H + ≈ m H3 in both benchmarks.    Table II.
We follow the mass ordering m H1 ≤ m H2 ≤ m H3 , hence H 1 is a sub-GeV s-like light state, H 2 the SM-like Higgs, and H 3 is a heavier scalar (with the physical states mostly unmixed).
Having established the scalar mass spectrum one can also determine the size of the scalar-fermion couplings, such as, the scalar-muon and -quark couplings which are needed to compute the contribution of the scalars to the charge radius ∆r 2 p . From Eq.(3), the scalar-muon couplings are given by where O R is the CP-even mixing matrix and is formally determined from Appendix A. Meanwhile, the scalar couplings to the proton are given by where y q is the quark Yukawa. Using Eqs. (14) and (15), one can write one coupling in terms of the other as The value of these couplings stays almost identical in either benchmark, since y = 4π in both.
Plugging the resulting couplings and scalar masses into ∆E (S) , as in Eq. (6), one gets that the Higgs-like H 2 and H 3 provide a negligible contribution to ∆r 2 p . The sub-GeV scalar H 1 adds to the required energy shift at the level of 10%, representing a subleading contribution to the charge radius. This makes plausible to stick to the limit in which the ∆r 2 p deviation is addressed by Z alone 7 .

B. The muon anomalous magnetic moment
As we have seen, our model is characterized by its mounphilic interactions, which give additional corrections to the muon anomalous magnetic moment, commonly parametrized as ∆a µ ≡ (g − 2) µ /2. Therefore, in what follows we show the allowed parameter space where our model simultaneously satisfies ∆r 2 p , the limits on the new scalars, and the constraints coming from latest (g − 2) µ results.
In light of the latest (g − 2) µ Fermilab measurement [4,19], we interpret these results either as a bound on BSM states under the assumption that the SM deviation eventually fades away or as a hint of new physics. This is done by substracting the (g − 2) µ measurement to two theory predictions, namely ∆a latt µ /10 −11 = 109 ± 71 , ∆a SM,th µ /10 −11 = 251 ± 59 .
∆a latt is the current Fermilab+BNL average after substracting the latest QCD theoretical prediction by lattice including vacuum hadron polarization [20]. A summary of more lattice results by other groups is discussed in Refs. [27,28]. On the other hand, ∆a SM,th µ is the current Fermilab+BNL average after substracting the theoretical SM contribution [29].
In this scenario, ∆a µ receives contributions from all physical scalars and the new gauge boson, all of them are listed in Appendix B. Since the Z couples exclusively to RH fermions, the contribution to (g − 2) µ is negative, given that the axial piece of ∆a µ is approximately 10 times larger that the vectorial one at equal couplings [30]. Fig. 5 shows ∆a (Z ) µ for different Z masses and gauge couplings, displayed as solid blue contours in the (g , m Z )-plane. According to the lattice bound (solid red) in Eq. (17), ∆a (Z ) µ < 0 is not excluded as long as its magnitude lies within its 2σ range, this would correspond to ∆a µ contours outside the red shaded region. However, notice that values allowed from the lattice constraint are incompatible with ∆r 2 p (solid black line). In other words, the safe region from the lattice constraint has smaller r p than the required one to fit the charge radius discrepancy, as noticed in Refs. [15,31] . For the scalar contributions to ∆a µ one needs to extract their couplings to muons. The H 1,2,3 couplings are already listed in Eq. (14) and those of the physical pseudoscalar and the charged scalar are expressed in terms of the corresponding mixing matrices given in Appendix A. Since y ν ∼ m ν /v Φ2 ∼ 10 −8 y, the couplings C µ S [H + ] and C µ S [H + ] are practically the same which leads to a negligible contribution from the charged scalar to ∆a µ .
As we have seen, the parameter κ controls the masses of the heavy scalars and tends to be of order O(1), look at Fig. 3 for reference. This leads to a high degree of degeneracy of the ϕ 2 -like scalars. Then, the contributions to (g − 2) µ by near-degenerate CP-even and CP-odd scalars partially cancel among each other, for comparable couplings. This is similar to what occurs with the contribution from the charged scalar [30]. As a result, the only sizable scalar contribution to the anomalous magnetic moment of the muon comes from the s-like state H 1 with sub-GeV mass. We then look for regions of the parameter space where its contribution to (g − 2) µ together with that of the Z satisfies the experimental constraints as well as explain the charge radius.
The left panel of Fig. 6 plots ∆a µ by the Z alone (in blue) as a function of its mass, where at each m Z the g coupling shown in the top axis fits ∆r 2 p according to Fig. 2. Even though g is varied, the required v s for the corresponding m Z in the bottom axis stays constant at the benchmark value as long as m Z is above a few MeV.
From the graph it is clear that ∆a (Z ) µ alone is too large and negative to hit either the (g − 2) µ lattice bound on new physics or the 2σ band of the tentative signal. Here is where the light s-like scalar remedies this situation: the black curves in the same panel indicate the combined ∆a µ by Z and s, evaluated at the scalar-muon couplings (16). One sees that for tens of MeV it is possible for the black curves to hit the region relevant for (g − 2) µ bounds. The right panel of Fig. 6 zooms in into the lattice bound (red dashed) and signal hint (green band) ballpark. Therefore, the model is capable of simultaneously fit ∆r 2 p and be compatible with (g − 2) µ measurements, while generating the muon mass and evading current collider constraints on heavy scalars.

C. Further restrictions on Z
Before concluding, we briefly mention constraints that may apply to our specific kind of Z . First, our Z expects to face looser constraints compared to other types of Z s that couple to electrons at tree-level 8 . For instance, this Z automatically avoids the stringent limits from rare pion decays π 0 → γZ → γ(e + e − ) and ν − e scattering.
µ-flavored Z s are however subject to a variety of limits [32]. For m Z exceeding the µ + µ − threshold, there are CMS and BaBar searches on Z Z → 4µ as well as B + → K + Z → K + (µ + µ − ) in LHCb. Below the µ + µ − threshold there are π 0 → γZ → γ(νν) decay constraints. Due to possible kinetic mixing and its modification to m Z , there are also bounds from atomic parity violation and the electroweak T -parameter.
More importantly, there are current stringent constraints on the (m Z , ε) plane, where ε = eg /(16π 2 ) is the kinetic mixing parameter. See for instance the limits by the NA64 and other collaborations in Ref. [33]. A simple analysis shows that the model is consistent in a wide region of the allowed parameter space. In addition, there are neutrinorelated bounds from non-standard induced ν-interactions, or coherent elastic neutrino nucleus scattering. At masses near 10 MeV there are also ∆N eff constraints on dark radiation degrees of freedom.
It is worth mentioning that a variety of nuclear experiments impose bounds on the ratio C n V /C p V of neutron-Z to proton-Z coupling, where C n For example, measurements on another muonnucleus bound state, the Muonic Deuterium (µD), put important exclusion limits on this ratio. From Ref. [26], one can see that results on µD allow the window −1 C n V /C p V 0 for m Z < 100 MeV. In our setup, see Sec. II, given that the U (1) R charge assignments for u and d are fully determined by q ν = 1, the ratio of neutron-and proton-Z couplings satisfies C n V /C p V = −1 for m Z in the 10-100 MeV range. This small tension with the µD constraint can be easily diluted if the u and d are placed in distinct anomaly-free sets. That is, for instance, taking another anomaly free solution {u R , b R , µ R , ν R } with U (1) R charge parametrized by q ν and {c R , d R , τ R , ν R } parametrized by q ν = q ν , then C n V /C p V = (2q ν − q ν )/(−2q ν + q ν ) can be varied and brought to safety from the µD limits while leaving our general conclusions for ∆r 2 p and (g − 2) µ unaffected 9 . A detailed analysis of the multitude of bounds on Z is left for an upcoming work [34]. The general Z constraints mentioned above must be implemented by taking into account 1) that Z is coupled to RH fermions only, and 2) that we have at our disposal some freedom in the quark sector to lock the RH quark rotations as diagonal matrices, and to turn off some u and d Yukawas to avoid dangerous H + -mediated meson processes. Before the seemingly surviving discrepancy in the proton charge radius succumbs to the improvement in precision techniques or better yet, a discovery is made, models accommodating it can be implemented and tested. This work proposes a SM extension based on a U (1) R gauge boson with exclusive coupling to right-handed muons, up and down quarks at tree-level, that addresses the ∆r p discrepancy for sub-GeV Z . Spontaneous breaking of the U (1) R , the generation of the µ, u, and d masses, and agreement with the recent (g − 2) measurements demand extra scalar degrees of freedom. The new scalar states can be brought to consistency with LHC limits albeit in a tight region of the parameter space, with a light singlet-like scalar at the tens of MeV, and doublet-like, near-degenerate neutral, charged and pseudoscalar states no lighter than ∼ 640 GeV. Further constraints arise for the gauge boson that deserve a study on its own, focusing on the preferential coupling of the Z to right-handed chiral fermions. so the extremum conditions are Upon minimization, the scalar mass matrices read and Define the rotation matrices O R and O I according to diag(m 2 H1 , m 2 H2 , m 2 Unless κ is tiny, neither the CP-even (squared) masses m 2 H k nor O R have a simple analytical form, they are numerically evaluated instead. Turning now to the charged sector we have, in the basis (ϕ + 1 , ϕ + 2 ), the following mass squared matrix whose eigenstates are the longitudinal W ± boson and the physical H ± . The (squared) masses and rotation matrix satisfying diag(0, m 2

Charged fermion sector
From Eq. (3), the textures of the charged fermion mass matrices are Given the number of Yukawas in the up and down quark mass matrices, the Cabbibo-Kobayashi-Maskawa mixing matrix can be straightforwardly fitted.
Appendix B: ∆aµ contributions The notation and normalization for the muon couplings in Secs. III and IV is adopted from Ref. [35], with S, P denoting true scalar and pseudoscalar, V and A polar and axial vector, and S, P refering to charged scalar couplings With vanishing charged lepton mixing, for i = j = µ ∆a (Z ) For the CP -even scalars H k , The analog pseudoscalar and charged scalar contributions are listed below,