Non-standard neutrino interactions and low energy experiments

We formulate an Effective Field Theory (EFT) for Non Standard neutrino Interactions (NSI) in elastic scattering with light quarks, leptons, gluons and photons, including all possible operators of dimension 5, 6 and 7. We provide the expressions for the cross sections in coherent neutrino-nucleus scattering and in deep inelastic scattering. Assuming single operator dominance we constrain the respective Wilson coefficient using the measurements by the COHERENT and CHARM collaborations. We also point out the constraining power of future elastic neutrino-nucleus scattering experiments. Finally, we explore the implications of the bounds for SMEFT operators above the electroweak breaking scale.


Introduction
In the SM the neutrinos interact with matter through exchanges of W and Z bosons. In addition, in the presence of new physics, the neutrinos could interact with matter via new mediators. Such Non-Standard neutrino Interactions (NSI) were contemplated already 40 years ago by L. Wolfenstein in the seminal paper on neutrino oscillations in matter [1]. Since then the NSI were studied extensively, but with a strong focus on neutrino oscillations, see, e.g., [2][3][4][5][6][7][8][9][10][11]. The bounds from neutrino oscillations are limited in scope, since they are JHEP09(2019)083 sensitive only to a subset of possible NSI. The common NSI effective Lagrangian relevant for neutrino oscillations contains only dimension 6 operators, see, e.g., [12],

JHEP09(2019)083 2 Operator basis for NSI
We are interested in the experiments where momenta exchanges are q O(100MeV), and thus well below the electroweak scale. The interactions of neutrinos with matter, i.e., with quarks, gluons, photons, electrons and muons, are described by an effective Lagrangian, obtained by integrating out the heavy degrees of freedom. These are the heavy SM particles: t, b, c quarks, τ lepton, W and Z bosons and the Higgs, as well as any heavy new physics particles. The interaction Lagrangian for ν α → ν β transition is given by a sum of nonrenormalizable operators, Here the C (d) a are dimensionless Wilson coefficients, while Λ can be identified, for O(1) couplings, with the mass of the new physics mediators. We consider a complete basis of EFT operators up to and including dimension seven. The sum in (2.1) runs over operator dimensions, d = 5, 6, 7, and operator labels, a, while in the notation we suppress the dependence on neutrino flavors α, β. The renormalization scale is fixed to µ = 2 GeV, unless specified otherwise.
For Dirac neutrinos the dimension 5 and dimension 7 operators, eq. (2.2) and eqs. (2.4)-(2.9), have a chirality flipping neutrino current. An incoming left-handed neutrino of flavor ν α is converted to a right-handed neutrino of flavor ν β . In contrast, the dimension 6 operators, eq. (2.3), preserve the chirality of the incoming neutrino. For Dirac neutrinos there are then two additional dimension 6 operators, Q 1,q , Q 2,q , obtained from (2.3) through P L → P R replacements. These operators cannot be well tested in neutrino experiments, since the production of right-handed neutrinos through SM weak interactions is neutrino mass suppressed. We therefore do not consider the operators Q (6) 1,q , Q 2,q further in our analysis.
In the case of Majorana neutrinos the dimension 5 and dimension 7 operators, eq. (2.2) and eqs. (2.4)-(2.7), violate lepton number by two units (note that we use the conventions of ref. [38] also for Majorana neutrinos). Furthermore, for a Majorana neutrino the operators Q (5) 1 in (2.2), Q 7,f in (2.7), and Q (7) 10,f , Q 11,f in (2.9) vanish identically for α = β, and thus only mediate flavor changing transitions. Finally, for α = β we include in the definitions of the operators an extra factor of 1/2 to compensate for the additional Wick contraction so that our results for cross sections and the bounds on Wilson coefficients can be used without change (cf. App A of ref. [39] for explicit normalization of such operators, albeit for DM interactions).
Note that in general the above operators are not Hermitian, and thus can have complex Wilson coefficients, C (d) a . The exception are dimension 6 operators with α = β, in which case the operators are Hermitian, and thus the corresponding Wilson coefficients are real (for these operators the "h.c." in eq. (2.1) should be dropped).
Note that the SM neutrino interactions with quarks are also described by the effective Lagrangian (2.1), though not all the operators are generated. The SM neutral currents (NC), i.e., due to the tree level Z exchanges, and the SM charged currents (CC), due to the tree level W exchanges, generate the operators Q 1,u(d,s) SM = ∓ 2,u(d,s) SM = ± 1,e SM = 2,e SM = − where s 2 W ≡ sin 2 θ W 0.2223 with θ W the weak mixing angle. The second terms in eq. (2.11) are due to CC, cf. figure 1. In the presence of NSI the above SM Wilson coefficients are modified tô C (6) 1 (2),f =Ĉ (6) 1 (2),f SM +Ĉ (6) 1 (2),f NSI , withĈ , (2.12) where in the last equality we used the ε notation of the NSI Lagrangian, eq. (1.1).
In the SM the dimension 5 and dimension 7 EFT operators, eq. (2.2) and eqs. (2.4)-(2.9), are suppressed by the neutrino masses and thus negligible for all practical purposes. In this case an appreciable Wilson coefficient would immediately signal the existence of NSI.

NSI and elastic scattering
This section describes the nuclear response to the elastic neutrino scattering on nucleus A at low energies, νA → νA, due to either the SM and/or NSI interactions. The calculation is done in several steps. In section 3.1 we first match onto an EFT describing neutrino interactions with non-relativistic protons and neutrons. The corresponding nuclear response to elastic neutrino scattering (CEνNS ) is given in section 3.2. For ease of comparison we also give the naive dimensional analysis (NDA) estimates for CEνNS cross sections induced by each of the EFT operators, while leaving the detailed numerical analysis for section 6.

Interactions of neutrinos with nonrelativistic nucleons
The neutrons and protons inside nuclei are non-relativistic and their interactions are well described by a chiral EFT with nonrelativistic nucleons. The momentum exchange, q, in CEνNS scattering is small so that nuclei remain intact, while neutrons and protons are non-relativistic throughout the scattering event. For instance, in the COHERENT experiment [13], the typical momentum exchange is q ∼ 30 − 70 MeV. This is well below the cut-off of chiral EFT, Λ ChEFT ∼ O(1 GeV), so that the effective neutrino interactions in eq. (2.1) can be included in the chiral EFT framework. We work at leading order in the chiral expansion for each of the EFT operators in (2.1), counting the light pseudoscalar masses to be parametrically of the order m π ∼ O(q). At leading chiral order the neutrino interacts only with a single nucleon, while interactions of a neutrino with two nucleons are JHEP09(2019)083 suppressed by powers of q/Λ ChEFT . The exception to this rule are the dimension seven Rayleigh operators, eq. (2.4), which we discuss separately in section 3.3. The effective Lagrangian describing neutrino interactions with non-relativistic nucleons is given by i,N (+h.c.), (3.1) where N = n, p, while d counts the number of derivatives in the operator, which gives the suppression of the operator in terms of soft momenta, O(q d ). The momentum exchanged, q µ = (q 0 , q), is given by, with k 1(2) , p 1(2) the incoming(outgoing) nucleon and neutrino momenta, respectively, cf. figure 2. The nuclear recoil energy, E R = q 2 /2m A , can, for fixed neutrino energy, be anywhere between E R,min = 0 for forward scattering, to a maximal value of E R,max 2E 2 ν /m A obtained in the case of neutrino back-scattering. The matching of quark and gluon currents onto nonrelativistic nucleon currents is performed using heavy baryon chiral perturbation (HBChPT) theory [40,41], while neutrino currents maintain their relativistic form. In this way one can explicitly show that the chiraly leading interactions of neutrinos are with a single nucleon current. We write the non-relativistic operators in the Lagrangian (3.1) using the heavy nucleon formalism of HBChPT, where the nucleon mass is effectively integrated out. To the order we are JHEP09(2019)083 working, the heavy nucleon field, N v , is given by where v µ is the nucleon four-velocity, which we may take to coincide with the lab frame, The momentum due to the heavy nucleon mass, m N v µ , has been factored out from the definition of N v by the exponential prefactor.
The nonrelativistic operators in (3.1) are, for proton and with d = 0, and a similar set of operators for neutrons with p → n. Here we have defined γ µ and the spin operator S µ N = γ 5 γ µ ⊥ /2. There are two relevant operators with a single derivative, d = 1, and one relevant d = 2 operator, where p ν 12 = p ν 1 + p ν 2 . We work in the isospin limit in which the proton and neutron masses are equal, so that m N = m p = m n 939 MeV. Above, the nucleon operators with σ µν ⊥ are related to the nucleon spin through where µναβ is the totally antisymmetric Levi-Civita tensor, with 0123 = 1. In (3.1) the Hermitian conjugation is present in the sum for almost all the operators -the exception are O 2,p for α = β, in which case the two operators are already Hermitian. The coefficients of these operators are thus real, while for the other operators they can be complex in general.
In COHERENT experiment the incoming neutrinos have energy ∼16-53 MeV, so thatĒ ν ∼ O(q) O(m π ). The coefficients for neutrons are obtained through p → n replacement. The form factors, F i , describe the hadronization of quark and gluon currents. They are functions of q 2 = − q 2 only, and their definitions are given in appendix A.
In order to derive the above expressions for the nonrelativistic coefficients in eqs. (3.9)-(3.15) we used the non-relativistic reduction of the nucleon currents summarized in appendix A. We keep only the leading terms in the q/m N expansion for each of the nonstandard neutrino interaction operators, eqs. (2.2)-(2.9). The leading contributions start at different orders in q expansion, depending on the structure of the NSI operators. For instance, the operators Q 4 , Q 6,q , Q 9,q , and Q 11,q , all match onto non-relativistic operators with one derivative, and thus their contributions to the scattering amplitude start only at O(q). All the other operators have contributions already at O(q 0 ).
There are two specific exceptions, where these leading contributions naturally vanish. For the QED dipole operator, Q 1 , the leading scattering on neutrons comes from operators with two derivatives. This is despite the Q 1 also contributing to the O(q 0 ) nonrelativistic operator, see eq. (3.11). The reason is that for the neutron q Q q F q/n 1 (0) = 0, so that in this case the contribution to (3.11) vanishes. Similarly, the absence of valence strange quarks in nucleons gives F s/N 1 (0) = 0, so that we need to keep the form factor F s/N 1 (q 2 ) expanded to quadratic order. Note as well, that the contributions from the axial current form factor F P (q 2 ) are proportional to the neutrino masses and thus neglected, even though F P (q 2 ) is 1/m 2 π enhanced due to the pion pole, corresponding to the middle diagram in figure 3.
In summary, in eqs. (3.9)-(3.15) most of the form factors are to be evaluated at q 2 → 0, 16) since this gives the chirally leading contribution, and we neglect the q 2 /m 2 N suppressed contributions, denoted above with the ellipsis. The three exceptions are the form factors while for the vector form factor of the neutron we need to go to O(q 2 ), (0) = 0. For simplicity we keep quadratic orders in F q/N 1 (q 2 ) also for q = u, d. 3 The input values of the nonperturbative parameters are given in [42].

Nuclear response to nonstandard neutrino interactions
The NSI coupling neutrinos to nonrelativistic nucleon currents, eq. (3.1), are of two types -the neutrinos either couple to the nucleon number operator,N v N v , or to the nuclear spin,N v S µ N N v . In the notation of ref. [43] the effective Lagrangian is where the ellipsis denote terms with P R ν, irrelevant for our case where the incoming flux is due to left-handed neutrinos. 4 The two Dirac structures are, Note that only the spatial three-vector components of l µ 5 enter the leading order nonrelativistic Lagrangian, (3.20). The EFT counting is such that all components of neutrino momenta count as the neutrino energy, E ν , while the nucleon currents are expanded in q/m N , as discussed in the previous subsection. The results for the c (d) i,N coefficients, that are in general functions of q 2 , are given in (3.9)-(3.15).
The cross sections for the neutrino-nucleus scattering is where E R is the recoil energy of the nucleus and the averaged amplitude square is given by [44], where J A is the spin of the target nucleus, the W M,Σ ,Σ (q) are the nuclear response functions andq ≡ q/q. The sum is over the isospin, τ, τ = 0, 1, with the kinematic factors 3 Incidentally, in this way we also capture the first nonzero term in c for scattering on neutrons, cf. eq. (3.11). In that case the leading term in c  for neutrino scattering on neutrons is described by c

JHEP09(2019)083
where the summations over spatial indices, i, j = 1, 2, 3, are implied. In the evaluation of the kinematic factors we only keep the leading terms in E R /E ν and E ν /m N , q/m N , which gives (note that E R = q 2 /(2m A )) Here τ, τ = 0, 1 denote the isospin so that, The non-relativistic coefficients describing neutrino interactions with protons and neutrons, c (3.32) Note that ∆σ NSI can be negative, if NSI interfere with the SM. In the last equality in (3.32) we show the parametric dependence on kinematical factors and nuclear response functions, The subscript is any of α = M, Σ , Σ , depending on the Wilson coefficientĈ (d) 1,q . In the long wavelength limit, q → 0, the nuclear response functions, W α , have the following parametric sizes, The response functions W τ τ Σ and W τ τ Σ encode the response of nucleus to the transverse and longitudinal axial operators, and thus measure the spin content of the nucleus. The values of W τ τ Σ ,Σ depend critically on the details of the nuclear wave function and can be much smaller than (3.33) for nuclei with all protons and neutrons paired. The W τ τ M response JHEP09(2019)083 functions count, in the long wavelength limit, the number of nucleons inside nucleus. This leads to coherent enhancement, also present for neutrino scattering through SM interaction -the tree level Z exchange, eq. (2.10). The Z boson couples most strongly to neutrons, so that in the SM case the enhancement is O(N 2 ), where N is the number of neutrons in the nucleus. Depending on the flavor structure the NSI can be due to couplings to proton or neutrons or both. The NSI operators, eqs. (2.2)-(2.9), fall into three categories: the operators that interfere with the SM contribution, the operators that do not interfere with the SM but still lead to coherently enhanced scattering, and the operators that are not coherently enhanced. The NDA estimates of the scattering cross sections for each of the three sets of operators are as follows.
The operators that interfere with the SM contribution. These are the operators with quark vector currents, Q 1,q in (2.3). The SM contribution to the corresponding Wilson coefficient is given in eq. (2.10). The NDA estimate of the NSI correction to the scattering cross section is where we used that the interference with the SM dominates over the purely NP contribution.
Coherently enhanced but no interference with the SM. The operators that lead to coherently enhanced scattering, but do not interfere with the SM contribution, are the ones that contribute to c operator in (2.5) that couples the neutrino current to gluons, the operator Q 5,q in (2.6) that couples neutrino and quark scalar currents, and the operators Q 10,q in eqs. (2.8), (2.9) that involve derivatives on the neutrino currents. The corresponding modification of the scattering cross section is parametrically, 3 + O(0.05)C 5,q + E ν Λ C 8,q + C 10,q + α EM 4π No coherent enhancement. The remaining operators do not receive coherent enhancement. The correction to the neutrino scattering cross section is then parametrically given by To shorten the expression we did not include additional numerical suppressions present for the case of strange quarks.

Scattering from Rayleigh operators
Finally, we include the estimates for the CEνNS induced by the Rayleigh operators, Q and Q 2 in (2.4). The CP even Rayleigh operator Q 1 leads to a coherently enhanced cross section, given by eq. (3.23) with the following matrix element squared [45] (for earlier work see [46,47]) (3.37) The first term is due to a contribution from two-body currents, calculated in ref. [45], where the two photon lines attach to two different protons in the nucleus, while the second contribution is due to both photon lines attaching to the same proton. The two contributions to the cross section scale as σ ∝ O(Z 4 Q 2 0 /m 2 N ) and σ ∝ O(Z 2 ), respectively (not showing the common factors and interference terms). For Q 0 /m N ∼ 0.1 the two contributions are parametrically of the same size for light nuclei Z ∼ O(10), while the first term dominates for heavy nuclei, Z ∼ O (50).
For the 2-proton form factor we use the phenomenological ansatz from ref. [45], The single nucleon matrix elements of the di-photon operators are not well known. We parametrize them as 1,N (q 2 ) the form factors. The NDA estimates for their values are, for With these definitions the contributions to neutrino scattering due to two photon exchanges with a single nucleon are obtained by setting the coefficient of the O 1 operator, while for the CP-odd Rayleigh operator, Q The CP odd Rayleigh operator leads to spin-dependent interactions, both from the single nucleon matrix element, (3.40), as well as from the 2 nucleon contributions. The two JHEP09(2019)083 nucleon contributions arise from one photon interacting with the proton charge, and the second photon with the magnetic moment of the other nucleon, be it proton or neutron. The single nucleon and two nucleon contributions to the cross section are parametrically σ ∝ O((q/m N ) 2 ) and σ ∝ O((ZqQ 0 /m 2 N ) 2 ) (not showing the common factors and interference terms). The two-body current contribution is thus expected to dominate for heavy nuclei, Z ∼ 50, while for light nuclei, Z ∼ 10 the single current contributions are important. The formalism for the two-body current contribution was worked out in [45], however, without deriving estimates for the resulting form factor. In the phenomenological analysis we thus conservatively ignore the 2-nucleon term and take as nonzero only the a (1) F ,N coefficient, using the NDA estimate in (3.41). While this estimate for the cross section does not capture the largest contribution for heavy nuclei, where the NDA suggest the cross section to be ∼ (ZQ 0 /m N ) 2 ∼ 20 times bigger for Z ∼ 50, the resulting error on bounds on Λ will be only ∼ (ZQ 0 /m N ) 2/3 ∼ 3, which suffices in view of other uncertainties in our estimates for this particular operator.

NSI and neutrino oscillations in matter
Neutrino oscillation data bound a subset of NSI -the ones that result in a nonzero forward scattering amplitude. The forward scattering of neutrinos on electrons and nuclei gives rise to matter effects in neutrino oscillations, described by an effective potential [48]. In this section we review the simplest case -electrically neutral unpolarized medium at rest. For extension to a polarized medium see [6], while for extensions to sterile neutrinos see, e.g., [49][50][51].
In the SM the effective potential receives contributions from both CC and NC. The NC contribution is neutrino flavor universal. It induces an overall phase shift in neutrino oscillation that is not observable (though it needs to be considered for oscillations into sterile neutrinos). The CC contributes to forward scattering of electron neutrinos on electrons. Electron neutrino scattering on an isotropic, homogeneous gas of unpolarized electrons is therefore described by the following effective Hamiltonian, see, e.g., [48], where n e is the number density of electrons. The resulting potential energy, leads to a change in the oscillation frequency for electron neutrinos. Here p is the neutrino momentum, α its flavor and h its helicity. The integral is performed over a finite volume V which is also included in the normalization of the states, , and thus drops out in the final result. Since weak interactions couple to left-handed fields, the h = −1 ultrarelativistic neutrinos obtain the effective potential energy while the positive helicity, h = +1, neutrinos are exposed to a negligible effective potential, V

JHEP09(2019)083
The above SM results are readily extended to NSI that couple neutrinos to the charged fermion vector current, i.e., the operators Q (6) 1,f in eq. (2.3). These lead to nonzero forward scattering amplitudes and thus induce an effective potential (the dependence on neutrino flavors is implicit) 1,e NSI n e − Ĉ (6) 1,u NSI + 2Ĉ 1,u NSI +Ĉ Here n p = n e is the number density of protons, equal to the number density of electrons in an electrically neutral medium, and n n the number density of neutrons. The global fits to neutrino oscillation data then lead to severe bounds on the parameters αβ (see section 6). The other NSI operators are poorly constrained from neutrino oscillations. This is most easily seen by analyzing the effects of NSI using the nonrelativistic basis, eqs. (3.4)-(3.7). The matrix elements of spin-dependent operators, O where n N is the nucleon number density. This gives the effective potential that is suppressed by the neutrino mass matrix, V and thus gives only extremely weak constraints on NSI.

NSI and deep inelastic scattering
For completeness we include the bounds on NSI that arise from deep inelastic neutrinonucleon scattering (DIS). While the DIS data were obtained at much higher momenta exchanges, q ∼ O(10 GeV), the constraints are severe enough that the EFT description may still be valid at least in parts of the parameter space. Throughout this section we thus assume that the EFT Lagrangian in eq. (2.1) is valid also for DIS. We comment on the validity of this assumption in section 6 where we confront predictions with data.
In neutrino-nucleus DIS the typical momentum exchange, q, is much larger than the inverse radius of the nucleus. The total cross section is therefore an incoherent sum of contributions from neutrino scattering on protons and neutrons inside the nucleus, The ellipses indicate power suppressed corrections to the factorization [52]. The sum runs over the partons, j = u, d, s, c, g, γ, with f j/N (x) the corresponding parton distribution functions (PDF) for the nucleon N . The kinematics of the process are shown in figure 4, with p 1 (p 3 ) the incoming neutrino (outgoing lepton) momentum, p 2 the incoming nucleon momentum a fraction x of which is carried by the parton, and we sum over the hadronic final state, X. We work in the limit of large neutrino energy, E ν , and large momentum exchanged, q 2 , i.e., E ν m N and q 2 m 2 N . The usual DIS variables, the partonic center of mass energy squared,ŝ, and the fraction, y, of incident neutrino energy transferred to the hadronic system, are, where E ν( ) is the out-going neutrino (lepton) energy. The double-differential neutrino-nucleon DIS cross section is, whereŝ = xs. We neglect the intrinsic charm content of the proton, take the off-diagonal CKM matrix elements to zero, while V ud , V cs → 1. In νq andνq (νq and νq) scattering, the total spin of the initial state is J = 0(1), so that the cross section is y independent (has (1 − y) 2 angular dependence). We will be interested in neutrino DIS measurements by the CHARM collaboration [53], which used the target material that is to a very good approximation isoscalar, i.e., composed of equal number of protons and neutrons. The average neutrino-nucleon scattering cross section for an isoscalar target is given by, For an isoscalar target the SM predictions for the CC neutrino-nucleon scattering, νN → JHEP09(2019)083 where we assumed isospin symmetry, (5.8) and similarly for antiquarks,q(x) ≡ fq /p (x). In (5.6), (5.7) we integrated out the W , and traded the m W dependence for the Fermi constant, G F .
In order to calculate DIS cross sections, we use the ManeParse package [54] to get the quarks and gluon PDFs (we use the CT10 NLO pdf set), while we take the photon PDF from [55].

Constraints on new neutrino interactions
Utilizing the results from sections 3-5 we now derive the bounds on the NSI Wilson coefficients, eq. (2.1), from neutrino oscillations, CEνNS [13], DIS [56] and from searches for neutrino dipole moment [57]. We also explore the reach of future CEνNS measurements in reactor experiments [58]. We restrict the analysis to interactions of ν e and ν µ since these are the NSI probed in CEνNS . The results are summarized in figures 6 and 7.

Constraints on NSI from oscillations
The oscillation data constrain the NSI contributions to the operators Q 1,f , eq. (2.3). The global fits to oscillation data allow at 95% C.L. [9] (in the notation of eq. (2.12)) −0.182 < ε eV ee < 0.264 , −0.120 < ε eV µµ < 0.120 , There are also bounds on operators that change neutrino flavor, or involve ν τ (for details see ref. [9]). The oscillations do not constrain NSI couplings to strange quarks, because the corresponding forward scattering matrix elements vanish. The above ranges on ε f V αβ imply the following lower bounds on the NP scale, for f = e(u, d), setting C

Constraints on NSI from CEνNS
Roughly a year ago the COHERENT collaboration measured for the first time the cross section for coherent neutrino-nucleus scattering [13] (see also data release in [59]). The target was 14.6 kg of CsI[Na], while stopped pion decays, π + → ν µ (µ + → e + ν eνµ ), acted as a source of neutrinos. The resulting time integrated neutrino fluxes per energy interval, φ ν i , are well known [19], Here m π (m µ ) is the charged pion (muon) mass, E ν the energy of the neutrino, and N = rN POT /(4πL 2 ) the time integrated neutrino flux, for each flavor, reaching the COHERENT detector. It depends on N POT = 1.73 × 10 23 , the delivered number of protons on target (POT), on r = 0.08, the number of neutrinos per flavor produced for each POT, and on L = 19.3m, the distance between the neutrino source and the detector. The expected number of CEνNS events for each neutrino flavor, α = ν e , ν µ ,ν µ , is where E R is the nuclear recoil energy, E ν the energy of the incoming neutrino, and n N the number of nucleons of type N = n, p, in the detector. The lower boundary in the integration over E ν is given by E ν,min ≈ M A E R /2, the minimal energy neutrinos need to have in order to induce nuclear recoil energy E R . The upper integration boundary, E ν,max , is given by the highest energy in the incoming neutrino flux. The ν e andν µ are produced in muon decay and thus have the maximal energy E ν,max = m µ /2, while ν µ is produced in pion decay, and has E ν,max ≈ 30 MeV. The maximal nuclear recoil energy deposited by ν e andν µ in the detector is thus E R,max 47 keV, while for ν µ it is E R,max 15 keV. The differential elastic neutrino-nucleus scattering cross section, dσ A /dE R , is given in eq. (3.23).
The prediction for the total number of CEνNS events expected in the COHERENT experiment is obtained by integrating eq. (6.9) over E R ∈ [0, 47] keV, convoluted with the signal acceptance fraction for COHERENT, given in figure S9 of [13] (which has an onset at about 4.25 keV). The experimentally allowed difference from the SM prediction then translates into bounds on the Wilson coefficients for the NSI operators,Ĉ In the numerical analysis we take only a single NSI Wilson coefficient at a time to be nonzero (apart from the SM contributions, eq. (2.10)). For simplicity we assume that the NSI affect either only ν e or only ν µ . To derive the 90 % C.L. allowed ranges forĈ 1 ,Ĉ 1,u/d ,Ĉ 2,u/d ,Ĉ 2,s ,Ĉ  follow the COHERENT collaboration [13,60] and define, 1,N . However, for heavy elements the latter is negligible, giving the estimate for σ α in table 1. The central values for W i are taken from [44]. Our estimates for the theoretical errors on W i are educated guesses. While this suffices at present, since these are subleading to the other uncertainties, a dedicated study would be desired in the future.
Our prediction for the SM rate in the COHERENT detector is N th Ĉ (d) a | SM = 188 ± 53 events. Comparison of this prediction with the COHERENT measurement gives the following 90% C.L. bounds on the NSI due to dimension 6 operators, −0.11 < ε uV ee < 0.49 , −0.10 < ε dV ee < 0.44 , (6.11) −0.06 < ε uV µµ < 0.12 , −0.06 < ε dV µµ < 0.11. For NSI couplings to the strange quark we obtain a relatively weak bound, |ε sV ee,µµ | 10 3 , because the sensitivity comes only from the O(q 2 ) term in the expansion of F s/N 1 (q 2 ), see eq. (3.19). The axial couplings ε qA ee and ε qA µµ are also poorly constrained, since they lead to spin-dependent interactions that are not coherently enhanced.
We collect the 90% C.L. bounds on the NSI coefficients in table 2 for ν e → ν X transitions and in table 3 for ν µ → ν X transitions. The bounds onĈ The bound onĈ (7) 7,s is controlled by the s-quark tensor form factor F s/N T,0 , which is not well known. At present it is even consistent with zero, which implies that there is no reliable bound on theĈ 1,2 where we (i) use the phenomenological form factor for two body currents forĈ (7) 1 , eq. (3.38), (ii) neglect the two-body current contributions forĈ (7) 2 , and (iii) and

JHEP09(2019)083
Lower bounds on Λ in GeV, ν e → ν X transitionŝ   use the NDA estimates (3.41) for the non perturbative single nucleon matrix elements, and do not assign any associated errors to these approximations. The bounds shown are thus just giving a rough potential reach of CEνNS experiments once theoretical errors will be under control (with probably a better guesstimate forĈ (7) 1 thanĈ (7) 2 ). Tables 2, 3 show in the case of dimension 6 operators the bounds for flavor non-diagonal processes ν e → ν X with X = e, and ν µ → ν X with X = µ, respectively. For flavor diagonal transitions, ν e → ν e or ν µ → ν µ , the bounds on NSI from the COHERENT measurement are instead 2,u(d,s) : Λ > 22.9(11.7, 4.5) GeV (ν e ), Λ >23.1(11.9, 4.6) GeV (ν µ ), (6.14) The relative sizes of bounds in tables 2, 3 are easy to understand using approximate scalings in eqs.
1,q , Q 3 , Q  Table 3. The 90 % C.L. lower bounds on Λ from COHERENT, CHARM and NaI 2T for ν µ → ν X (X = µ, for X = µ see main text) NSI Wilson coefficients,Ĉ 10,q , eq. (3.35), translate into more stringent bounds for Λ. The bounds are significantly weaker for the remaining non-enhanced operators. The bounds on operators with strange quarks are also weaker, since the corresponding form factors are smaller.
The bounds on NSI from CEνNS experiments are set to improve in the future with a number of new detectors either already taking data or being planned. The COHERENT collaboration is operating a 10kg Ge detector, a 22kg single-phase liquid Ar detector, and a 185 kg NaI[Tl] scintillating crystal detector [61]. The liquid Ar may increase to 1ton, and NaI[Tl] to 2 tons, in the future [62], cf. table 4. To take full advantage of these experimental progress an increased precision in the predictions of nuclear response functions and nuclear form factors will be called for. The sensitivity reach of the proposed upgrades at COHERENT is studied in details in ref. [63].
In tables 2 and 3 we show the expected improvements in the sensitivity to NSI due to the 2 ton NaI detector proposed by the COHERENT Collaboration [61], with the same JHEP09(2019)083 COHERENT CHARM NaI 2T Borexino C 1 (5) C 1 (7) C 2 (7) C 3 (7) C 4 (7) C 5,u (7) C 5,d (7) C 5,s (7) C 6,u (7) C 6,d (7) C 6,s (7) C 7,u (7) C 7,d (7) C 7,s (7) C 8,u (7) C 8,d (7) C 8,s (7) C 9,u (7) C 9,d (7) C 9,s (7) C 10,u (7) C 10,d (7) C 10,s (7) C 11,u (7) C 11,d (7) C 11,s ν e COHERENT CHARM NaI 2T Borexino C 1 (5) C 1 (7) C 2 (7) C 3 (7) C 4 (7) C 5,u (7) C 5,d (7) C 5,s (7) C 6,u (7) C 6,d (7) C 6,s (7) C 7,u (7) C 7,d (7) C 7,s (7) C 8,u (7) C 8,d (7) C 8,s (7) C 9,u (7) C 9,d (7) C 9,s (7) C 10,u (7) C 10,d (7) C 10,s (7) C 11,u (7) C 11,d (7) C 11,s neutrino source but with a baseline of 28m and a lower energy threshold of ∼ 13 keV. Furthermore, in the projections we assume that the total theoretical uncertainty is reduced 10-fold compared to the present ones, quoted in table 1. This would give the projected total theoretical uncertainties σ F i ∼ 3 − 5%. This will require more precise determinations of the neutrino flux, which is already planned, as well as a much better knowledge of the quenching factors, and major advances in the purely theoretical inputs -the form factors and nuclear response functions entering the SM prediction. While such a decrease of uncertainties may be aggressive, they do give us a useful gauge of the potential reach of CEνNS experiments. In NaI detector the neutrino recoils on both the iodine and sodium nuclei. For the SM neutrino interactions the scattering on iodine completely dominates. The coherently enhanced cross section is ∼ 40 times larger for neutrino scattering on iodine as it is for sodium. This is the case also for coherently enhanced NSI interactions, where scattering on iodine similarly dominates. Spin-dependent interactions, on the other hand, can be comparable, depending on the operator. We find that scattering on iodine dominates except for the operators Q 2,q , Q 7,q and Q 9,q , for which the main contribution to the scattering rate is from sodium, while the two cross sections are of the same order for Q (5) 1 . The expected bounds from ν e,µ → ν X scattering are shown in tables 2 and 3.
For flavor diagonal transitions, ν e,µ → ν e,µ , the expected bounds on dimension 6 operators from the NaI 2T detector are,   Table 4. A list of proposed experiments to detect CEνNS using (anti)neutrinos from Stopped Pion Decay (SPD) or from reactors, with recoil energy threshold, T th , the distance from the source, the target material and its mass given in 2nd to 4th columns.
While COHERENT uses stopped pions as a source of neutrinos, there are also a number of planned or already operating CEνNS experiments that use reactor antineutrinos, see table 4 as well as, e.g., refs. [15,71,72]. Reactors produce large quantities of low energy electronic antineutrinos. On average about ∼6 antineutrinos are produced per fission, for a total of ∼ 2×10 20ν e per second per GW of thermal reactor power [72][73][74], with a maximum energy of ∼ 8 MeV.
As two representative examples of reactor CEνNS experiments we chose the proposed RED100 [58] and MINER [64] experiments, and checked their respective sensitivities to different NSI, assuming in both cases a total uncertainty of σ F i ∼ 10%. RED100 has a proposed 100kg target of liquid Xenon, with a baseline of 19m from a 1GW reactor and an energy threshold of 500 eV. While this energy threshold is lower than for stopped pion decay experiments, it is the highest among the reactor experiments. MINER [64] has a proposed 30kg detector composed of 72 Ge and 28 Si in 2:1 ratio, with a baseline of 1m from a 1MW detector and a very low energy threshold of 10 eV.
Because of lower thresholds, both RED100 and MINER would have one to two orders of magnitude better sensitivity to the neutrino dipole moment,Ĉ 1 , compared to NaI 2T. The reach for the vector NSI current operators,Ĉ 1,u(d) could exceed the ones from global oscillation fits, while there would be also an appreciable improvement on the derivative couplings,Ĉ 8,q ,Ĉ 10,q . The operators inducing spin-dependent interactions, on the other hand, cannot be probed better in these experiments, since xenon has a smaller nuclear spin than iodine, while 72 Ge and 28 Si have nuclear spin 0.
Future experiments can improve their sensitivity to specific NSI operators by changing experimental conditions. As just stated, lowering the energy threshold improves the sensitivity to the magnetic monopole operator Q (5) 1 , whose contribution is enhanced by the 1/ q 2 photon pole, cf. eqs. (3.11), (3.35). For instance, the CONNIE collaboration proposes a 1kg Si detector with the energy threshold of 28 eV that would be situated 30m away from JHEP09(2019)083 the reactor [65]. Taking all the other parameters as in our projection for the NaI 2T limit, leads to a projected limit for Q (5) 1 of Λ 50 TeV, to be compared with 8.3 TeV at NaI 2T, table 2.
Varying the neutrino energy, E ν , as well as increasing the range of nuclear recoil energies, E R , can also be beneficial. In this way one could uncover q 2 dependence of the neutrino scattering cross section that differs from the SM one. We illustrate this in figure 8, taking as the example the Q  4 (the NP scale was taken very low in figure 8 to exaggerate the effect).

Constraints on NSI from deep inelastic scattering
The CHARM neutrino detector [75] was composed of 78 plates of marble (CaCO 3 ) with a total fiducial mass of 87.4 tons. The target is to a very good approximation an isoscalarthe correction to the cross section from the isotriplet component is O(0.2%) [56]. Data were recorderd exposing the detector to neutrinos and antineutrinos from the CERN 400 GeV SPS proton beam dump.

JHEP09(2019)083
The ratios R e and R νµ are predicted in the SM to be [77,78] R SM e = 0.3221 ± 0.0006, R SM νµ = 0.3156 ± 0.0006. (6.19) The dominant theoretical uncertainty in the two predictions is due to the approximation that the target was taken to be an isospin singlet, which is correct within O(0.2%).
Since the neutrino CC cross section is strongly constrained, we can assume that NSI only affect NC transitions. The ratio R e in eq. (6.19) then receives the NSI correction as where σ CC is the total neutrino and antineutrino CC cross section and ∆σ NSI is the NSI contribution to the NC cross section. In the analysis we take only a single NSI operator at a time to be nonzero. Similarly, R νµ is modified to The relative values of bounds are easily understood from the matrix elements in eqs. (5.13)-(5.17). In particular, the two scalar operators Q (7) 5,q and Q (7) 6,q have exactly the same matrix elements and thus the same bounds for given flavor q. The difference between the bounds on Λ for three light quark flavors comes predominantly from the factor m q that is part of the definition of the operators, leading to ∼ (m s /m u,d ) 1/3 larger Λ exclusion for the strange quark. The matrix element of the tensor operator, Q 7,q , is bigger, cf. eq. (5.17), which translates to roughly factor of 2 more stringent bounds on Λ.
In the case of flavor diagonal transitions, the bounds on dimension 6 operators are For the reader's convenience we translate the CHARM bounds on dimension 6 operators also into the bounds on ε i parameters, eq. (1.1). For the electron neutrino, one obtains

JHEP09(2019)083
For the case of the electron neutrino these bounds are comparable, yet somewhat stronger, than the bounds from oscillations for vector currents involving u and d quarks, eqs. (6.2) and (6.3), and are significantly stronger for the case of the muon neutrino. Note that for a number of operators the bounds on Λ from the CHARM experiment are comparable to the momentum exchange q ∼ O(10 GeV). This means that the EFT analysis may be applicable only for strongly coupled mediators, with couplings larger than O(1). Another general comment regarding CHARM constraints on NSI is that for light mediators these are comparatively less effective than CEνNS constraints where the momentum exchange is smaller.

Other constraints
The contribution to the neutrino scattering rates due to the neutrino magnetic moment has a pole at q 2 = 0, cf. figure 9 in appendix B. This means that experiments with lower E R thresholds will have better sensitivity to the magnetic moment. Furthermore, scattering on electrons will in general lead to lower q 2 . Measurements of solar neutrinos scattering on electrons in Borexino [57] give the current most stringent limits on the magnetic moment µ να , or, equivalently, on the Wilson coefficient C 1 : Λ > 2.7 · 10 6 GeV (ν e ); Λ > 1.8 · 10 6 GeV (ν µ ). (6.30) The measured neutrino scattering rates in Borexino can also be translated in a bound on Rayleigh operators C 1,2 . The neutrino interactions mediated by the two Rayleigh operators result either in νA → νA scattering through 1-loop matrix elements or in νA → νAγ, i.e., with an emission of an extra photon. In the Borexino experiment both processes lead to the same signal. Using the results from section 3.3 for the 1-loop contribution, with NDA estimates for the single nuclear matrix elements and neglecting two-body currents for Q A different set of bounds on NSI operators (2.3)-(2.9) comes from collider experiments. For a proper analysis we need to extend the EFT analysis to above the electroweak symmetry breaking scale, which we do in the next section. However, some of the bounds directly apply to the operators (2.3)-(2.9). For instance, the searches for dark matter can be reinterpreted as bounds on NSI operators with two neutrinos replacing the two DM particles in the final state. Finally, Borexino Phase II results also bound dimension 6 NSI for couplings to electrons [79] − 0.14 < ε eV ee < 0.21, (6.33) where we conservatively quote the combined allowed ranges marginalized over the highand low-metallicity input choices for the standard solar model. The resulting bounds are JHEP09(2019)083 somewhat more stringent than from neutrino oscillations, eq. (6.1). In terms of the NP scale the bound in (6.33) gives Λ > 0.76 TeV.
The monojet searches at CMS and ATLAS [80][81][82] thus result in lower limits of Λ 1 TeV for (axial-)vector current operators, eq. (2.3), Λ 200 GeV for gluon-gluon operators, eq. (2.5), Λ 40 GeV for scalar operators, eq. (2.6), and Λ 20 GeV for tensor current operators, eq. (2.7), after converting to our normalization of the operators. The monophoton searches bound di-photon operators, eq. (2.4), giving Λ 30 GeV [83]. These bounds, apart from (axial-)vector current operators, are quite weak, with values of Λ allowed that are even below the kinematical cuts on p T and/or MET. One may thus question the applicability of the EFT with the actual bounds from colliders in reality even weaker, unless the couplings are large. We do not attempt to correct for these effects since this is beyond the scope of present manuscript, see, however, [84] on how to properly obtain EFT bounds from LHC searches.
Once the EFT is uplifted above the electroweak scale, the bounds from collider searches can become more severe. For instance, searches for charged fermion contact interactions at LEP [85] give constraints on the SU(2) L ×U(1) Y symmetric operators, defined in eqs. (7.4)-(7.7) below. These operators result in dimension 6 NSI operators (2.3) below the electroweak scale, but with bounds on Λ from LEP of about O(1 TeV) even when coupling only to leptons (see section 7 for details).

NSI above the electroweak scale
In this section we turn to the question of how the NSI interactions (2.1) are generated. Since the bounds on many of the operators are relatively mild, cf. tables 2 and 3 , it is possible that light NP could be responsible for their generation. This interesting direction was pursued, e.g., in refs. [34,35,86,87].
The other option is that the NP responsible for NSI is heavy, heavier than the electroweak scale. If this is the case, the NP states can be integrated out leading to an EFT that is valid between the scale of NP, Λ, and the electroweak scale, v EW , with the effective Lagrangian (see also [88][89][90][91][92]) where Q d a are SU(2) L × U(1) Y invariant operators. In constructing the electroweak (EW) EFT operators we add to the SM field content the right-handed neutrino, ν R , which is a SM singlet. This allows for SM neutrinos to be either Dirac or Majorana. In the rest of this section we discuss the bounds on operators Q d a and their matchings onto the low energy EFT for NSI, eq. (2.1). We only consider operators that do not violate lepton number, since the lepton number violating operators are severely constrained by bounds on neutrinoless double beta decay.
The two dimension 6 operators in EW EFT that lead to an neutrino magnetic dipole moment are (throughout this section we do not display generational indices on leptonic fields, and assume flavor conservation for quark currents) 2) JHEP09(2019)083 whereH = iσ 2 H * , with H the Higgs doublet, L L the lepton doublet, B µν and W a µν the hypercharge and SU(2) L field strengths, and τ a the SU(2) L generator in the fundamental representation. Both of the above operators contribute below the EW scale to the dimension five magnetic dipole operator, eq. (2.2), where v EW 246 GeV is the Higgs vev. If either Q 1,B or Q 1,W are generated in the UV, the neutrinos will have magnetic moments. The exception is, if the two contributions cancel against each other, i.e., for 2Ĉ 3F,ii and should be dropped), which will then give for the Wilson coefficients of dimension 6 operators below electroweak scale (u 1 = u, d 1 = d, d 2 = s and α the lepton flavor, not displayed on l.h.s.) 1F,α1 + 1 4 C 2F,α1 + C 6F,α1 , (7.8) 1F,αi + 1 4 C 2F,αi + C 7F,αi , (7.9) 3F,α1 + 1 4 C 4F,α1 + C 5F,α1 . (7.10) The problem with generating large contributions toĈ (6) 1(2),f Wilson coefficients in this way is that the dimension 6 operators in eq. (7.4) are extremely well bounded. Translating the results from [78,93] the bounds for Q   21 , but a precise determination would require a correlation matrix to properly account for the change of basis). For electron-quark couplings the most stringent bounds come from measurements of d → ueν transitions and atomic parity violation, and are much more severe than the bounds from neutrino scattering. For muon-quark couplings the most stringent bound for Q (6) 2F,21 is from d → uµν transitions, while for the other operators, Q 1F,2i , Q 6F,2i , Q 2F,i , with all the other Wilson coefficients zero, since then the quarks only couple to neutrinos.
Another option is that the leading contributions arise from dimension eight operators with two Higgs insertions (see also, e.g., [94]), After the Higgs obtains a vev only the neutrino is projected out of the lepton doublet, The Wilson coefficients for low energy dimension 6 operators are in this case, , (7.14) Since these operators only lead to couplings of quarks to the neutrinos and not to charged leptons, this relaxes some of the bounds. The remaining bounds are "inevitable", as they come from processes that involve neutrinos -these are the bounds discussed at the end of section 6.4. The gauge-gauge operators in (2.4) and (2.5) can arise from dimension eight operators giving the Wilson coefficientŝ The reinterpretation of the 8 TeV ATLAS W +MET search, ref. [95], bounds Q 4W operators of roughly comparable strength, but with the details dependening on the relative sizes of photon and Z exchange contributions [98].
The operators (2.6) and (2.7) can arise from the following dimension six operators

JHEP09(2019)083
The resulting Wilson coefficients arê 5,e = Re C 2m e Λ 2 ,Ĉ 6,e = Im C 2m e Λ 2 ,Ĉ 7,e = C (6) 6R 2m e Λ 2 . (7.25) The bounds on these chirality flipping operators are quite stringent. For instance, from pp → + / E T + X searches at the LHC we can expect bounds on the NP scale of the operators (7.21)-(7.23) at the order of Λ O(5TeV) and at the similar level from semileptonic decays of light pseudoscalar mesons. (This estimate is based on the bounds for tensor current SM-EFT operators with left-handed neutrinos obtained in [99], which only give quadratic corrections to the SM rates, as do the operators (7.21)-(7.23). A dedicated analysis of operators with right-handed neutrinos is called for, which, however, is beyond the scope of our work.) No such bounds exist for dimension eight operators, 12H = ν R σ µνH † τ a L L L L τ a Hσ µν e R , (7.31) which, after the Higgs obtains the vev, have the form (ν R . . . ν L )(f . . . f ). The relevant constraints on these operators thus come only from neutrino scattering experiments, discussed in the previous section. In principle there are also constraints from Higgs decaying to four body final states, h → 2jνν. However, these are at present much less constraining. The Wilson coefficients of the low energy EFT operators (2.6) and (2.7), generated from the above operators, arê 8H,i , (7.32) 10H,i , (7.33) 11H ,Ĉ 12H . (7.34) Finally, the dimension 7 low energy EFT operators with derivatives on the neutrino current, eqs. (2.8) and (2.9), can arise from the following electroweak EFT dimension 8 operators, After the Higgs obtains the vev, the operators have the form of (neutrino current)×(charged fermion current). The most relevant bounds on these operators are, again, due to neutrino scattering experiments discussed in the previous sections. The matching gives for the low energy Wilson coefficientŝ 10{11}, 10{11}, 10D ,Ĉ 10{11},e = v EW Λ 4 Re{Im} C 9D . (7.42)

Conclusions
In this manuscript we obtained predictions for CEνNS in the presence of nonstandard neutrino interactions described by an EFT at 2 GeV. Our analysis covers the complete basis of EFT operators up to and including dimension 7, and thus covers most of the viable models as long as the mediators are heavier than about 100MeV. We recast the recent measurement of CEνNS by the COHERENT collaboration using a CsI detector to obtain bounds on the EFT operators, assuming that only one NSI operator at the time contributes appreciably. The main results are collected in figures 6 and 7, where they are compared with the bounds on NSI from neutrino oscillations, the solar neutrino flux measurements at Borexino, and from deep inelastic scattering. The obtained bounds apply to incoming electron or muon neutrinos scattering either through flavor diagonal interaction, or even, if the neutrino flavor changes (including scattering to sterile neutrinos). We see that already now the CEνNS measurements lead to the most stringent limits for some of the NSI operators, for instance for scalar currents. The NSI reach of CEνNS experiments is set to significantly improve in the future, with a number of new experiments either already running or being planned. In figures 6 and 7 we also show the projected limits for the NaI 2 ton detector planned by the COHERENT collaboration, also assuming that systematic errors can be decreased by an order of magnitude. The new experiments, as we show in the paper, can also increase their sensitivity to NSI by modifying running conditions, such as the incoming neutrino energy, the nuclear recoil energy thresholds, but also by trying to perform measurements at higher recoil energies. For this it will be important to investigate how well one can distinguish between NSI and the subleading corrections to our predictions (for instance from q 2 dependence of nuclear form factors away from zero recoil point, see, e.g., [100]).
When using our results it is important to note their validity. The DIS bounds assume that the mediators are heavier than a few 10s of GeV, the COHERENT bounds that they are heavier than about 100 MeV, while Borexino and oscillation bounds apply also to very JHEP09(2019)083 light mediators (a typical momentum exchange in solar neutrinos scattering on electrons in Borexino is q ∼ few 100 keV−1 MeV, which sets the lower bound on mediator mass, such that use of EFT is justified for interpreting Borexino measurement). For light mediators, with masses below 10 GeV the bounds from DIS would get suppressed, and similarly for COHERENT bounds for mediators lighter than a few 10s MeV. It is easy to recast our bounds also in such cases, but one does need at least simplified models for the mediators in that case (some examples are, e.g., [101][102][103]). More challenging would be to extend our work to the intermediate range of q 2 ∼ (few) GeV 2 , where neither ChPT methods that we used, nor the factorization used for DIS, apply. The benefit, on the other hand, is that many of the neutrino experiments are taking data precisely in this theoretical difficult intermediate regime, so that any theoretical advances would be highly desired.

Acknowledgments
We thank Gil Paz for the discussions regarding Rayleigh operators, and Rex Tayloe

A Nucleon form factors and nonrelativistic limits
The nucleon form factors, F i , in eqs. (3.9)-(3.15), are defined as [42], where we suppressed the dependence of nucleon states on their momenta, N | ≡ N (k 2 )|, |N ≡ |N (k 1 ) , and similarly,ū N ≡ u N (k 2 ), u N ≡ u N (k 1 ), while k µ 12 = k µ 1 + k µ 2 and q µ = k µ 2 − k µ 1 , and γ [µ q ν] = γ µ q ν − q µ γ ν .  Figure 9. The ν e → ν X (left) and ν µ → ν X (right) rates in the COHERENT CsI detector including Q 1 NSI operator contribution, forĈ The form factors, F i , are functions of q 2 only. In the derivations of eqs. (3.9)-(3.15) a non-relativistic reduction of nucleon currents is required. Counting v · ∂ ∼ O(q 2 ), with q the typical soft three-momentum, the leading terms are [42] where Λ a;min is the current lower limit for this operator, as obtained in section 6.2 from the COHERENT measurement, and listed in tables 3 and 2.
The predicted differential scattering rates, dN/dE R , eq. (6.9), are plotted in figures 9-12. For the nuclear response functions, W i , we use the values from [44], while the value of nuclear form factors are taken from [42]. In figures 9-12 we show separately the scattering rates due to the ν e (left panels) and ν µ +ν µ (right panels) incoming neutrinos. The

JHEP09(2019)083
corresponding fluxes are given in eqs. (6.6)-(6.8). Note from eq. (3.11) that the operators Q (7) 8,q and Q (7) 10,q have the same matrix element and will give rise to the same differential scattering rate; same happens for the operators Q (7) 9,q and Q (7) 11,q , see eq. (3.13). The neutrinos are due to stopped muons, which sets the maximal recoil energy, E R,max , to be around 47 keV for ν e → ν X andν ν → ν X , and about 15 keV for ν µ → ν X transitions. The muon neutrino flux is monoenergetic, with E ν ≈ 30 MeV. As a consequence, forν µ → ν X scattering there is an abrupt drop in the predicted differential rate, dN/dE R , at E R ≈ 15 keV, since none of theν µ can contribute to more energetic recoils. This discontinuity is clearly visible for Q (6) 2,q , see figure 10 bottom right, Q 1 , . . . , Q 4 , see figure 11 right, and for Q (7) 5,q , . . . , Q  (3.29). This prefactor goes to zero when when the maximal E R for given value of E ν is reached, i.e., when the incoming neutrinos backscatter. This means that the contribution fromν µ → ν X to the SM dN/de R scattering rate goes to zero at E R ≈ 15 keV.
In order to obtain the number of events predicted in the CsI[Tl] COHERENT experiment the predictions in the left and the right panels of figures 9-12 need to be added up, and then convoluted with the signal acceptance fraction of the detector. At present the acceptance has a lower threshold at around 4.25 keV, denoted as a vertical red line in figures 9-12.
From the figures we see that a number of operators lead to a significantly different E R dependece compared to the SM predictions. The magnetic dipole moment leads to dN/dE R that has a pole at q 2 = 0, clearly showing a significant increase in the rate at small values of E R , see figure 9. Lowering the energy threshold can thus lead to an increase sensitivity to this NSI operator, as long as the background can be kept low. In other case probing larger recoils may be beneficial. For instance, the scattering matrix element due to Q 1,s is proportional to F s/N 1 (0) q 2 and thus grows with the increased E R , see figure 10 (top). A dedicated analysis is required to see to what extend this contribution can be distinguished from the subleading corrections in the SM rate that come from the q 2 dependence of the F u/N 1 and F d/N 1 form factors, and from the uncertainties in the nuclear response function W M .
A very striking difference in the kinematical dependence of dN/dE R arises in the case of monoenergetic neutrino beams, as already mentioned above. This can be seen in figures 10-13 (right panels). Observing experimentally any such discontinuity would be a clear signal for the presence of NSI. The discontinuity if especially pronounced for Q (7) 2 , Q 4 , figure 11, and Q 6,q , figure 12, since these operators contribute to the non-relativistic operator O (1) 1,N , see eqs. (3.13), (3.40). This leads to an additional q 4 dependence in the scattering rate, see eq. (3.29). For NSI generated by Q 2 , Q 4 , or Q 6,q operators the sensitivity increases for higher E R recoils. Q 4 , or Q (7) 6q operators the sensitivity increases for higher E R recoils. This was also illustrated in the main text, in figure 8.  Figure 10. The ν e → ν e (left) and ν µ → ν µ (right) scattering rates in the presence of a single NSI operator, either Q 1,q (top) or Q 2,q (bottom), settingĈ (6) a = 1/Λ 2 a;min , where Λ a;min is taken to be the corresponding lower bound in tables 2 and 3, respectively. The SM event rate is denoted by the black solid line, dashed lines denote NSI only and solid lines the sum of SM and NSI contributions (with blue, brown and purple representing couplings to up, down and strange quark currents, respectively). The energy threshold of the experiment is denoted by the vertical red line.
The operators Q (7) 1 , Q 3 and Q (7) 5,q match onto the nonrelativistic operator O (0) 3,N which leads to a q 2 prefactor in the scattering rate instead of the 4E 2 ν − q 2 one for the SM, see eq. (3.28). This different E R dependence clearly shows in figures 11 and 12 (left panels). Similar comment applies to the Q (7) 9,q operator, which matches onto O (1) 2,N , which leads to a kinematic prefactor E 2 R in the scattering event rate, cf. eq. (3.29), and thus a very different different recoil dependence compared to the SM rate, see figure 13.

C Numerical values of CEνNS scattering cross sections
In this appendix we provide numerical expressions for CEνNS differential cross sections in the presence of a single nonzero NSI operator Q (d) a . We normalize the NSI cross sections to the SM, so that they take the form  Figure 11. The ν e → ν X (left) and ν µ → ν X scattering rates in COHERENT CsI detector in the presence of NSI operators Q 1 , Q 2 , Q 3 , Q 4 (top to bottom) settingĈ  Figure 12. The ν e → ν X (left) and ν µ → ν X scattering rates in COHERENT CsI detector in the presence of NSI operators Q 5,q , Q 6,q , Q 7,q , (top to bottom) settingĈ (7) a = 1/Λ 3 a;min . The notation is as in figure 10.
Below we give the numerical values for the coefficients g for CEνNS on nuclei 23 Na, Ge, 127 I and Xe. For Germanium and Xenon we calculate the average cross sections for natural abundance of stable isotopes, namely 70 Ge, 72 Ge, 73 Ge, 74 Ge and 76 Ge for Germanium and 128 Xe, 129 Xe, 130 Xe, 131 Xe, 132 Xe and 134 Xe for Xenon. For the numerical evaluation we use the nuclear response functions from [43,44]. We only quote central values for the coefficients g a , but comment when these estimates are particularly uncertain.
Since only dimension 6 operators interfere with the SM amplitude for CEνNS , these are the only ones that have both g 8,q (top) and Q 9,q (bottom).
the lowest order in E R in the expressions of g (d) a,q (E ν , E R ) and f (d) a,q (E ν , E R ). These quoted results for these functions are thus reliable for recoil energies up to E R ∼ 10 − 20 keV, while for higher energies one needs to take into account corrections from higher powers of E R .

C.1 Numerical results for CEνNS on 23 Na
To shorten the notation we define the following functions of incoming neutrino energy, E ν , and nuclear recoil, E R , D = (E ν /MeV) 2  We first give the results for NSI parametrized by ε i , eq. (1.1). Using ε i , instead of the Wilson coefficientsĈ