Non-standard neutrino interactions in $U(1)'$ model after COHERENT data

We explore the potential to prove light extra gauge $Z^\prime$ boson inducing non-standard neutrino interactions (NSIs) in the coherent-elastic neutrino-nucleus scattering (CE$ \nu $NS) experiments. We intend to examine how the latest COHERENT-CsI data can constrain this model. A detailed investigation for the upcoming Ge, LAr, and NaI detectors of COHERENT collaboration has also been made. Depending on numerous other constraints coming from oscillation experiments, muon $ (g-2) $, beam-dump experiments, LHCb, and reactor experiment CONUS, we explore the parameter space in $Z^\prime$ boson mass vs coupling constant plane.


I. INTRODUCTION
Oscillations of neutrinos among different flavors are now a well established phenomenon from various experimental searches, which implies that the neutrinos carry non-zero masses and their different flavor are substantially mixed [1]. Currently, we have fairly good understanding of all the neutrino oscillation parameters in three-flavor paradigm, except the Dirac CP violating phase [2][3][4]. This led us into an era of precision measurements in the leptonic sector, where it is possible to observe sub-leading effects originating from physics beyond the Standard Model (SM). Furthermore, this may affect the propagation of neutrinos and eventually it may impact the measurements of three-flavor neutrino oscillation parameters. Among various new physics scenarios beyond the standard three-flavor neutrino oscillations, non-standard neutrino interactions (NSIs) can be induced by the new physics beyond the SM (BSM). In literature, they are traditionally described by the dimension-6 four-fermion operators of the form [5], where f αβ represent NSI parameters and α, β = e, µ, τ , f = e, u, d. The importance of NSIs were discussed well before the establishment of neutrino oscillation phenomena by a number of authors where f C αβ are NSI parameters, α, β = e, µ, τ , C = L, R, denotes the chirality, f = e, u, d, and G F is the Fermi constant 3 . The Hamiltonian in presence of matter NSI, in the flavor basis, can be written as, where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [53], ∆m 2 ij = m 2 i − m 2 j (i < j = 1, 2, 3), and A ≡ 2 √ 2G F N e E represents the potential due to the standard matter interactions of neutrinos and αβ can be written as  [54] where αβ = | αβ |e iφ αβ for α = β. In general, the elements of αβ are complex for α = β, whereas diagonal elements are real due to the Hermiticity of the Hamiltonian as given by Eq. (3). For the matter NSI, αβ can be defined as, where N f is the number density of fermion f and f C αβ = f L αβ + f R αβ . In case of the Earth matter, one can assume that the number densities of electrons, protons, and neutrons are equal (i.e. N p N n = N e ), in such a case N u N d 3N e and one can write, In Table I, we give recent constraints for NSIs obtained from a combined analysis of oscillation experiments and COHERENT measurements [54] at 2σ C.L.
Having introduced general descriptions of NSI and its bounds, in next section we describe our model in great details.

III. THE SETUP
In this work, we extend the SM gauge group to an anomaly free where α and β can be e, µ and τ . In our framework, two different lepton flavors are coupled to the new gauge interaction. The relevant charge assigmnents for the lepton fields as well as the scalar fields that trigger the U (1) gauge symmetry breaking are listed in Table II. In this prescription, we include two scalar fields φ 1 and φ 2 transforming as 1 and 2 under U (1) , respectively. These scalar fields are responsible for the U (1) breaking and therefore to give mass to the Z gauge boson.
The scalar potential for the fields in our framework (Table II) can be split in three parts, The first part is the SM Higgs potential, the second is the coupling of the Higgs doublet with the the singlet fields, and the third part is the potential for the two singlet fields, Now once new scalar fields φ i attain their vev (v i / √ 2), we get mass for the Z gauge boson as where we have used charges for the φ i as mentioned in Table II. Now to give an order of estimation about the breaking scale, we take mass of Z gauge boson M Z = 0.1 GeV, whereas coupling strength g is taken as ≈ 2.8 × 10 −5 . Note that we consider these numerical values in such a way that these can be probed in future COHERENT experiments (for a detail discussion see Sec. IV and Fig. 3). Using these numerical values in Eq. (10), one finds the vevs of φ 1 and φ 2 as v 1 ≈ 3 TeV and v 2 ≈ 1 TeV, respectively. It is worth to mention that the Higgs vacuum stability can be obtained when coupled to singlet scalar field whose vev is at the TeV scale [55].
There are several anomaly-free solutions to the U (1) involving baryon number, for scenarios where CEνNS and NSI have been explored, see for instance [14][15][16]46]. In our approach, we choose the anomaly-free solution for the U (1) = U (1) B−2Lα−L β , see appendix A for details. In this case, lets take one of the solutions, namely (x e , x µ , x τ ) = (0, −1, −2) for instance, the righthanded (RH) neutrino lagrangian is given by The six possible assignments under U (1) for the leptons are depicted in Table III and each of the charge assignments give rise to a different model, namely different neutrino masses and mixings as well as different NSI.

IV. COHERENT ELASTIC NEUTRINO-NUCLEUS SCATTERING
Coherent elastic neutrino-nucleus scattering has already been measured by the COHERENT experiment [20], using a scintillator detector of made of CsI. The low energy neutrino beam was generated from the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory.
The SM differential cross section for CEνNS is given by [12,56,57] where T is the nuclear recoil energy, E ν is the incoming neutrino energy, M N is the nuclear mass, and Q 2 w is the weak nuclear charge, given by where Z(N ) is the proton(neutron) number, F Z(N ) (q 2 ) its nuclear form factor, and g V p = 1/2 − 2 sin 2 θ W , g V n = −1/2, are the SM weak couplings. It is important to notice that the cross section depends highly on the mass of the detector and type of material, especially on the number of neutrons N , since the dependence on Z is almost negligible due to the smallness of g V p (∼ 0.02). For the SNS energy regime, another important feature of the CEνNS cross section are the nuclear form factors. From here on, we will adopt the Helm form factors [58] with an average neutron rms radius of R n = 5.5 [29].
The differential recoil spectrum can be computed as Here, N A is Avogadro's number, M det is the detector mass, M m is the molar mass of the material, and φ α (E ν ) is the neutrino flux for each flavor. The SNS neutrino flux consists of monochromatic ν µ coming from π + decays, along with delayed ν e andν µ from subsequent µ + decays. Each of these flux components are given by to the experimental measurements [20].
where r = 0.08 is the fraction of neutrinos produced for each proton on target, N POT represents the total number of protons on target, and L is the distance to the detector.
From Eq. (14) we can compute the expected number of neutrinos per energy bin: where A(T ) is the acceptance function, taken from the COHERENT data release [59]. In Fig. 1 we show the measured number of events from the COHERENT collaboration as a function of the nuclear recoil energy T, along the expected number of events in the SM framework.
In presence of NSI, the cross section for CEνNS is affected through the weak nuclear charge (see Eq. (13)) in the following way: where α = (e, µ, τ ). Notice that with this new contribution, the differential cross section from Eq. (12) is now flavor dependent.
It is possible to write an effective low-energy Lagrangian for the neutrino-fermion interactions with the Z boson: where Q 2 is the transferred momentum. Therefore, by comparing this effective Lagrangian with the NSI Lagrangian in Eq. (2), we can relate the NSI parameters with the Z interaction parameters as ε qV αα =    [20], along with the future setups using other type of detectors [60].
In Fig. 2, we plotted the number of events versus the nuclear recoil energy, for different detectors (Ge, LAr and NaI) from the upgrade plan of the COHERENT collaboration. The features of the future detectors that are used in our numerical simulations, along with the current CsI detector are presented in Table IV. We show the expected events in the SM framework, and compare with the case of NSI terms in the cross section. For this particular example, we considered the A 2 model, and the values M Z = 0.1 GeV and g = 2 × 10 −5 . We give the details about the values of M Z and g that are considered here in Fig. 3. As expected, the number of events in the presence of NSI increases with respect of those in the SM, but this increase is higher for smaller values of the nuclear recoil energy T .
Given the relation in Eq. (21), it is now clear how the NSIs can be generated from the interactions of a new vector boson Z . By computing the number of events including NSI contributions, we are now able to compare with the COHERENT measurements in order to set boundaries to the coupling and mass of the Z boson.
As mentioned before, the first part of the analysis consists in comparing with the first measurements of CEνNS, provided by the COHERENT collaboration [20]. A CsI detector of 14.6 kg was used at a distance of 19.3 m from the source. The cross section for this type of detector has to be computed separately for cesium (Cs) and iodine (I) in the following way: We performed a fit of the COHERENT data by means of a least-squares function where N i meas (N i th ) is the measured (expected) number of events per energy bin, B i is the background, and σ stat = N th + B i is the statistical uncertainty. We marginalized over the nuisance parameters α and β, which quantify to the signal and background normalization uncertainties σ α and σ β , respectively. Following the COHERENT analysis, we choose σ α = 0.28 and σ β = 0.25 [20].
Since the fit to the quenching factor was done for the bins from i = 4 to 15, we follow our analysis only for those energy bins.
In order to extract information about the Z boson, we compute the expected number of events N th including NSI effects, according to Eq. (18) and considering the weak nuclear charge in Eq. (19).
It must be pointed out that in the NSI scenario, the differential cross section is now flavor dependent.
As we have mention before (see section III for details), the proposed model has six possibilities on the U (1) charges of the charged leptons. Since only four of these cases are allowed by oscillations data (A 1 , A 2 , B 3 and B 4 ), we will perform the χ 2 analysis only for these cases.
Since all the quarks have same U (1) charge, we get ε uV αα = ε dV αα , reducing the number of free parameters. Also, the neutrino source does not produce tau neutrinos, and hence, we can not extract any information about ε τ τ .
In Fig. 3, we show the exclusion regions at 95% C.L. in the (M Z , g ) plane. Each panel corresponds to one of the four possible models, where the resulting light neutrino mass matrix is of type A 1 , A 2 , B 3 and B 4 . The constraint coming from the current COHERENT data has been presented using the light-green shaded region. In order to have a more complete study, we also corresponds to the constraint set by the current COHERENT measurement using a CsI detector [20]. The solid purple line shows the limit from oscillation experiments [54]. The limits set by the future detectors setup from the COHERENT collaboration [60], namely, Ge, NaI, and LAr are shown using red dash-dotted, yellow dotted, and blue dashed lines, respectively. The limits from the CONUS reactor experiment [22] are shown by the magenta (long dashed) lines. The exclusion regions set by beam dump experiments [61][62][63][64][65][66][67][68][69][70], BBN and CMB [71], and LHCb dark photon searches [72] are presented using color code yellow, gray and sky-blue regions, respectively. The pink shaded band corresponds to the region where the muon (g − 2) anomaly is explained [73] (see text for more details).
included exclusion regions coming from the future upgrades of the COHERENT collaboration: Ge, NaI, and LAr detectors, considering an exposure of four years. We show them using red dashdotted, yellow dotted, and blue dashed lines, respectively. We can see how these future setups can improve the current COHERENT limit for the coupling g by almost one order of magnitude.
In Fig. 3, we also included limits coming from oscillation experiments (see purple solid line) using the relation For models A 1 and A 2 , we take the smallest value of ε µµ from the first column of Table I, when setting ε ee = ε τ τ = 0. Then we use Eq. (24) to get a limit for g as a function of M Z . For B 3 , we have directly a value for ε ee because ε µµ = 0 5 . The limit from BBN +CMB [71] is also presented using gray band in Fig. 3. For the cases where ε ee = 0 (i.e., for B 3 and B 4 ), we have also included boundaries for a light Z boson, obtained by different electron beam dump experiments as shown by yellow region. We have used the Darkcast [74] code to translate the beam dump limits to our specific model. In the cases where ε µµ is present, we also considered limits set by dark photon searches for LHCb limits [72] shown using sky-blue region. We also used the Darkcast [74] code to translate these limits to the different cases of our model, which has been shown using sky-blue regions in Fig. 3.
Since the existence of new light vector bosons can explain the inconsistency in the anomalous magnetic moment of the muon [75,76], we have incorporated in our plots from Fig. 3 the region of the (M Z , g ) plane where our model can explain the discrepancy ∆a µ = (29 ± 9) × 10 −10 [73] (pink region). The interaction of the Z boson with muons leads to an additional contribution to the anomalous magnetic moment: where Notice that only in the B 3 this region is absent, since there is no interaction between muons and the Z boson (ε µµ = 0).
Furthermore, there are several proposals aiming to measure CEνNS using nuclear reactors, such as CONNIE [77], CONUS [22], MINER [78], RED100 [79], TEXONO [80], etc. For example, the CONUS experiment will consist of a 4 kg Germanium detector, located at 17 m from the nuclear power plant at Brokdorf, Germany [22]. They expect ∼ 10 5 events over a 5 year run, assuming the SM signal. 5 For model B4 we considered the smallest possible value between εee and εµµ from the other models.
Furthermore, we present limits for the Z boson considering the CONUS experiment. Since a nuclear reactor produces only electron antineutrinos, we give an exclusion regions only for the cases where ε ee = 0 (i.e. for B 3 and B 4 ). These regions are shown in the lower panels of Fig. 3   The first panel of Fig. 3, i.e., A 1 (U (1) B−Lµ−2Lτ ) has µµ and τ τ with τ τ > µµ . In this scenario, it can be seen that the future COHERENT experiment with LAr detector will explore a parameter space for masses between 7 MeV to 3 GeV and couplings as small as g ∼ 10 −5 .
For masses between 11 MeV and 3 GeV the future COHERENT bounds will be competitive with the current LHCb exclusion limits. However, we notice that above 3 GeV bounds coming from LHCb drak-photon searches will give the strongest constraints, where g can be ∼ 10 −3 (see skyblue region). Bounds arising from the calculation of ∆N ef f of BBN will rule out M Z < 7 MeV as shown by the gray band. We now proceed to discuss our results for A 2 (U (1) B−2Lµ−Lτ ) as shown by the second panel of the first row of Fig. 3. It has µµ and τ τ as in A 1 but in this case τ τ < µµ . Here, we have found that the future COHERENT experiment will explore a parameter space for masses between 7 MeV to 0.55 GeV and couplings up to g ∼ 10 −5 . For masses between 11 MeV and 15 MeV the future COHERENT bounds will be comparable as exclusion coming from LHCb. Unlike A 1 , LHCb can explore more parameter space for this scenario, i.e., M Z ≥ 0.55 GeV, compared to COHERENT-LAr bounds. ∆N ef f also shows similar bounds as A 1 .
Unlike the scenarios A 1 and A 2 , we also have contributions coming from beam dump experiments and reactor experiment CONUS that is because of non-zero ee for B 3 and B 4 , which we show at the second row of Fig. 3, respectively. The model B 3 (U (1) B−Le−2Lτ ) predicts NSI parameters like ee and τ τ (see Table III  By investigating all the four scenarios, it has been seen that the bounds arising from (g−2) µ (see pink band ) is ruled out by current COHERENT-CsI data, while limits from oscillation experiments (as shown by the solid purple line) will be ruled out by the future COHERENT data. Finally, we present a set of benchmark values that can be explored by different experiments in the Table IV.

V. TWO-ZERO TEXTURES
Here we revisit the phenomenology of two-zero textures that are allowed in this model, as given in Table III, viz A 1 , A 2 , B 3 , and B 4 , in light of the latest global-fit data. Two-zero textures are phenomenologically very appealing in the sense that they guarantee the calculability of the neutrino mass matrix M ν from which both the neutrino mass spectrum and the flavor mixing pattern can be determined [81][82][83]. In what follows, we first parameterize M ν in terms of three neutrino mass eigenvalues (m 1 , m 2 , m 3 ) and three neutrino mixing angles (θ 12 , θ 23 , θ 13 ) together with three CP violating phases (δ, α, β). Note that here δ is the Dirac type CP-phase, whereas α, and β are the Majorana type CP-phases. Therefore, the mass matrix M ν can be diagonalized by a complex unitary matrix U as where m diag ν = diag{m 1 , m 2 , m 3 }. In the standard PDG formalism, the neutrino mixing matrix U , also known as PMNS matrix is given by where s ij = sin θ ij and c ij = cos θ ij . Given the parameterization of U , it is now straight forward to write down the elements of neutrino mass matrix M ν with the help of Eq. (27).
The two-zero textures of the neutrino mass matrix M ν (see Eq. (27)) satisfies two complex equations as where a, b, p and q can take values e, µ and τ . Above equations can also be written as where V has been defined in Eq. (28). We notice that these two equations involve nine physical parameters m 1 , m 2 , m 3 , θ 12 , θ 23 , θ 13 and CP-violating phases α, β, and δ. The three mixing angles (θ 13 , θ 12 , θ 23 ) and two mass-squared differences (∆m 2 12 , ∆m 2 23 ) are known from the neutrino oscillation data. Note here that from the latest global-fit results, we have some predictions about the CP-violating phase δ, however at 3σ, full range i.e., 0 • − 360 • is still allowed. Therefore, in this study we kept δ as a free parameter. The masses m 2 and m 3 can be calculated from the known mass-squared differences ∆m 2 12 and ∆m 2 23 using the relations m 2 = m 2 1 + ∆m 2 12 , and m 3 = m 2 2 + ∆m 2 23 . Thus, we have two complex equations relating four unknown parameters viz. m 1 , α, β and δ. Therefore, one can have the predictavility of all these four parameters within the formalism of two-zero textures.
We numerically solve Eq. (30) for the concerned types of two-zero textures, see Table III. It has been known from the latest global analysis of neutrino oscillation results [2][3][4] that the least unknown parameter among the three mixing angles is the atmospheric mixing angle θ 23 . Therefore, considering two benchmark values of θ 23 , we calculate remaining unknown parameters, which we present in Table VI and VII. We take the latest best-fit value of θ 23 from [4] as one of our benchmark value, whereas maximal value of θ 23 i.e., θ 23 = 45 • is taken as the second benchmark value. Notice that the seed point θ 23 = 45 • has the great importance in perspective of flavor symmetries as well as flavor models building. Among numerous theoretical frameworks, µ − τ symmetry that explains θ 23 = 45 • has received great attention in the neutrino community, for the latest review see Ref. [84].
From Table VI, we notice that the textures A 1 , A 2 can explain both the latest best-fit as well as the maximal value of θ 23 . Further, given these benchmark values we calculate unknown parameters m 1 , α, β and δ. It is to be noted from the fourth column that the predicted values of δ for all the cases lies within 1σ of the latest best-fit value [4], which is 237.6 +37.8 • −27.0 • . We also calculate m 2 , m 3 (see second column of the  For the textures B 3 and B 4 , one can have non-zero |m ee | (see Table II are dedicated to look for the signature of 0νββ-decay are namely, GERDA Phase II [86], CUORE [87], SuperNEMO [88], KamLAND-Zen [89] and EXO [90]. It is to be noted here that, this process violate lepton number by two-units and the half-life of such decay process can be read as [91,92], where G 0ν is the two-body phase-space factor, and M 0ν represents the nuclear matrix element 10 -4. 10 -3. 10 -2. 10 -1.
10 -1. We present our predictions for the effective Majorana neutrino mass |m ee | for both the textures in Fig. 4. The 3σ allowed parameter space of |m ee | considering the latest global-fit data [4] for the normal neutrino mass hierarchy is shown by the light-orange band 6 . The magenta band shows the latest bounds on |m ee |, arises from KamLAND-Zen 400 experiment [89] which is read as |m ee | < (61 − 165) meV at 90% C.L. by taking into account the uncertainty in the estimation of the nuclear matrix elements. We also show the the first results of KamLAND-Zen 800 collaboration using the lighter-green band, which was presented in the latest meeting TAUP 2019 [93]. Besides this, the predictions for |m ee | for the textures B 3 and B 4 are shown by the blue (cyan) patch at 3σ (1σ) significance level. We notice from both the panel of Fig. 4  predictions can be probed by the KamLAND-Zen 800 data. 6 Note that the present oscillation data tends to favor normal mass hierarchy (i.e., ∆m 2 31 > 0) over inverted mass hierarchy (i.e., ∆m 2 31 < 0) at more than 3σ [2][3][4], therefor, we focus only on the first scenario.
It is to be noted here that the latest bound on the sum of neutrino masses m ν come from   Fig. 3, whereas other neutrino phenomenology are given in Fig. 4 and in Table VI, VII. Depending on our analysis, we make our final remarks as follows: • Texture A 1 : in this case, we notice that the future COHERENT experiments with NaI or LAr detectors will explore a parameter space for masses 7 MeV ≤ M Z ≤ 3 GeV within the coupling limits 0.8 × 10 −5 ≤ g ≤ 10 −3 . Also, parameter space below 5.3 MeV can be ruled using the measurement of ∆N ef f coming from the observation of Big Bang nucleosynthesis.
Notice here that this observation holds true for remaining cases. Furthermore, it can be seen that above 3 GeV LHCb can put the strongest bound. Also, in this scenario, one does not have prediction 0νββ-decay as it gives m ee = 0 .
• Texture A 2 : findings of A 2 is similar as A 1 . However, we notice that the future COHERENT experiments will show tightest constrain upto the mass limit ∼ 550 MeV and above this LHCb will give the stringent bound. It is to be noted here that LHCb can exclude more parameter space for A 2 compared to A 1 , which is simply because µ− field carry 2-units U (1) charge than of A 1 (in case of A 1 , U (1) charge of µ− field is 1).
• Texture B 3 : outputs of B 3 is very different compared to A 1 and A 2 . Here we notice CEνNS experiment like, CONUS can explore the most of the parameter space for the masses of M Z above ∼ 25 MeV and coupling constant g ≥ 5 × 10 −6 . On the other hand, below this has been ruled out by beam dump experiments.
Moreover, one also have predictions for 0νββ-decay which can be explored by KamLAND-Zen collaboration (see left panel of Fig. 4).
• Texture B 4 : in this case CONUS can rule out the parameter space for the mass range, 25 ≤ M Z ≤ 500 MeV corresponding to copuling strength 5 × 10 −6 ≤ g ≤ 1.5 × 10 −4 . Above this mass limit and coupling strength LHCb can put the tightest constraint. Moreover, beam dump experiments can exclude the parameter space below 25 MeV. We also have predictions for 0νββ-decay and the parameter space are marginally consistent with the present limit of both KamLAND-Zen and Planck bound as given in the right panel of Fig. 4. In this section we will show an example of how to compute the light neutrino mass matrix, for a specific choice of U (1) charges. By considering a type-I seesaw scenario [95], the neutrino mass matrix is given by where M D and M R are the Dirac and Majorana mass matrices, respectively.
In our prescription, the Yukawa Lagrangian invariant under SM ⊗ U (1) for charged-leptons and neutrinos is given by −L Y ⊃ y e L e e H + y µ L µ µ H + y τ L τ τ H + y ν 1 L eH N 1 + y ν 2 L µH N 2 + y ν 3 L τH N 3 .
This leads to the Dirac neutrino mass matrix of the form For the Majorana mass matrix, we need to specify the U (1) fermion charges. For example, with the choice (x e , x µ , x τ ) = (0, −1, −2), the RH neutrino Lagrangian is Therefore, the Majorana mass matrix takes the form Plugin M D and M R in Eq. (B1) leads to a light neutrino mass matrix of the form which corresponds to the type A 1 neutrino mass matrix. One can follow the same procedure for the other charge assignments to get different mass matrices (see Table III for details).