EFT analysis of New Physics at COHERENT

Using an effective field theory approach, we study coherent neutrino scattering on nuclei, in the setup pertinent to the COHERENT experiment. We include non-standard effects both in neutrino production and detection, with an arbitrary flavor structure, with all leading Wilson coefficients simultaneously present, and without assuming factorization in flux times cross section. A concise description of the COHERENT event rate is obtained by introducing three generalized weak charges, which can be associated (in a certain sense) to the production and scattering of $\nu_e$, $\nu_\mu$ and $\bar{\nu}_\mu$ on the nuclear target. Our results are presented in a convenient form that can be trivially applied to specific New Physics scenarios. In particular, we find that existing COHERENT measurements provide percent level constraints on two combinations of Wilson coefficients. These constraints have a visible impact on the global SMEFT fit, even in the constrained flavor-blind setup. The improvement, which affects certain 4-fermion LLQQ operators, is significantly more important in a flavor-general SMEFT. Our work shows that COHERENT data should be included in electroweak precision studies from now on.


Introduction
Precision measurements play an ever-increasing role in particle physics. A broad range of observables has been designed to test various aspects of the Standard Model (SM), such as accidental symmetries, CP violation, flavor structure, electroweak symmetry breaking, etc. Beyond the SM (BSM), precision measurements are sensitive to new particles and interactions, often well beyond the direct reach of high-energy colliders.
Experiments where neutrinos are detected constitute a vital ingredient in the precision program. Neutrino scattering on matter (electrons, nucleons, nuclei) probes the interaction strength and structure of the weak force. In the SM context it is a measure of the weak mixing angle, although currently it cannot compete with the much more precise determination based on the W mass and Z-pole physics. More generally, it is sensitive to new BSM particles interacting with neutrinos.
Coherent elastic neutrino scattering on nuclei (CEνNS) is the newest member of the family of electroweak observables. The process was theorized long ago [1,2], but the experimental challenges were overcome only recently by the COHERENT experiment [3]. The peculiarity of CEνNS is that vector interactions between neutrino and nuclear constituents (nucleons or quarks) add up coherently. Consequently, possible small deviations between the actual neutrino interaction strength and that predicted by the SM are also amplified. This effect, roughly proportional to the number of neutrons in the nucleus, makes CEνNS a powerful probe of physics beyond the SM.
Assuming the BSM particles are heavy, the most convenient description of the nonstandard effects in CEνNS involves the effective field theory (EFT) framework [4,5]. There is in fact a ladder of relevant EFTs, with each consecutive rung involving higher energy scales, different degrees of freedom, and more stringent assumptions about the scale of new physics [6,7]. At the energy scales typical for CEνNS, nucleons (protons and neutrons) are convenient degrees of freedom. In this case, experiments probe the strength of effective 4-fermion contact interactions between neutrinos and nucleons, that is to say, neutron and proton weak charges in the standard parlance [8]. Moreover, the observed CEνNS rate is obviously also sensitive to BSM effects in neutrino production. In the SM, the CEνNS event rate is simply proportional to the square of the weak charges weighted by the number of protons and neutrons in the target nucleus [9]. At higher energies, above the QCD phase transition but below the electroweak scale, quarks become the relevant degrees of freedom. In this language CEνNS experiments probe the Wilson coefficients of the effective 4-fermion neutral-current (NC) interactions between neutrinos and quarks [10], as well as additional charge-current (CC) interactions involved in neutrino production (in pion and muon decay in the case of COHERENT experiment). We refer to this effective theory as the Weak EFT (WEFT) [11]. 1 The Wilson coefficients of operators involving neutrinos can be translated into the language of non-standard interactions (NSI) [14] prevalent in the neutrino literature.

JHEP05(2023)074
Finally, above the mass of the Z boson the relevant effective theory is the so-called SM Effective Field Theory (SMEFT) [16,17], where the degrees of freedom are those of the SM and the full SU(3) C × SU(2) L × U(1) Y gauge symmetry is realized linearly. This is the most used EFT in the particle physics community because it can be directly (and now automatically [18]) matched to a host of popular BSM models with supersymmetric particles, leptoquarks, W'/Z' bosons, etc.
In this paper we interpret the results of the COHERENT experiment in the language of the EFTs. We analyze the entirety of the available COHERENT data to extract the nuclear and nucleon weak charges (which we properly generalize so that production effects are taken into account). We translate these into constraints on the WEFT Wilson coefficients. Our results are provided in a completely general form that is straightforward to apply to any BSM model or EFT analysis. In fact, we analyze the impact of current CEνNS data in a global SMEFT fit to electroweak precision observables (EWPO), which include such emblematic data as the W boson mass measured in hadron colliders, Z boson decays measured in LEP-1, or atomic parity violation in cesium. We will demonstrate that this impact is non-trivial, significantly reducing the allowed parameter space for the SMEFT Wilson coefficients. Moreover, from the point of view of EWPO and SMEFT fits, we will show that COHERENT is the most relevant neutrino-detection experiment to date! There are several ways in which this paper improves on the previous literature regarding the EFT approach to COHERENT data [9,10,[19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]: • Our work, together with the recent refs. [33,34], are the first ones to use the entire COHERENT dataset presently available (2D energy and time distributions with argon and cesium-iodine targets) [35,36] to probe fundamental interactions, and thus represents the current state of the art.
• We perform for the first time a complete analysis in the framework of WEFT, including the nonlinear effects of the various Wilson coefficients in detection and production.
The latter, as well as the non-trivial interplay between the two, have not been correctly discussed before. We remark that correlated effects in production and detection are generic in new physics models, since the SU(2) L gauge symmetry relates CC and NC interactions.
• We identify the linear combinations of EFT Wilson coefficients that are strongly constrained by COHERENT.
• The model-independent SMEFT analysis of COHERENT data and the combination with other EWPO is performed for the first time.
The rest of this article is organized as follows. In section 2 we describe the EFT framework, and in section 3 we obtain the EFT prediction for the event rate measured at COHERENT. In section 4 we discuss the various COHERENT measurements and how they will be implemented in our numerical analysis, which is presented in section 5 using the WEFT setup. The implications for the SMEFT and the combination with Electroweak Precision Observables is discussed in section 6. Finally, section 7 contains our conclusions.

EFT below the electroweak scale
Given the low energies involved in CEνNS, our starting point will be the most general effective Lagrangian below the electroweak scale with the SM particle content, the WEFT Lagrangian [11]. In this effective theory, electroweak gauge bosons, the Higgs boson, and the top quark have been integrated out and the electroweak symmetry is explicitly broken. Let us stress that right-handed neutrinos are not part of the WEFT particle content.
Here we focus on the lepton-number-conserving parts of the WEFT Lagrangian that are relevant for COHERENT physics, namely the NC interactions mediating CEνNS and the CC interactions involved in pion and muon decay. Let us first present the NC interactions between neutrinos and quarks [11,13]: where P L,R = 1 ∓ γ 5 /2 are the chirality projection operators and 1 is the unit matrix. The dimensionful normalization factor is related to the measured value of the Fermi constant G F ≈ 1.166 × 10 −5 GeV −2 [37] via v ≡ ( √ 2G F ) −1/2 ≈ 246. 22 GeV. In the coefficients of the operators we separated the SM values g qq V,A from the new physics contributions qq V,A . The former are given at tree-level by where T 3 q and Q q are the weak isospin and charges of the quark q and θ W is the weak mixing angle. See ref. [8] for the discussion of radiative corrections and numerical values of these SM couplings. The parameters qq V,A are zero in the SM limit and more generally they are Hermitian matrices (in the neutrino indices). As is well known [1], the contribution of the vector Wilson coefficients to the coherent neutrino scattering rate is enhanced by (A − Z) 2 , where A − Z is the number of neutrons in the nucleus. On the other hand, the contribution of the axial Wilson coefficients does not benefit from such an enhancement, and moreover it vanishes completely for scattering on spin-zero nuclei. In view of the current experimental uncertainties, we can thus neglect the axial contributions in the following, leaving us only with the vector Wilson coefficients. 2 To lighten the notation the vector subindex will be omitted and we will work instead with the more compact notation qq αβ ≡ [ qq V ] αβ . Let us now introduce the parts of the WEFT Lagrangian that describe the interactions relevant for neutrino production at COHERENT. Pion decay is described by [11,39] 3)

JHEP05(2023)074
where α = e, µ, τ are charged lepton fields, V ud is the (1,1) element of the Cabibbo-Kobayashi-Maskawa (CKM) matrix. Here, new physics is parametrized by ud L,R,S,P,T , which are general complex matrices in the lepton indices. Finally, the WEFT interactions describing muon decay are [11,39]: Here new physics is parametrized by the [ρ L,R ] aαβb tensors. Hermiticity of the Lagrangian implies [ρ L ] * aαβb = [ρ L ] bβαa , while [ρ R ] aαβb is a general complex tensor. Contrary to eqs. (2.1) and (2.3), in the case of eq. (2.4) we cannot normalize the interactions using v. This is because that parameter is defined via G F , which is extracted from the experimental measurement of the muon decay rate, to which the interactions in eq. (2.4) contribute. Instead, we use the tree-level value of the Higgs vacuum expectation value, v 0 . Using eq. (2.4) it is trivial to calculate the muon decay width and relate v and v 0 , finding v 2 In other words, we are working with an input scheme such that the µ → eνν decay rate is controlled, exactly as in the SM, by v (equivalently, by G F ), without any tree-level dependence on the ρ L,R Wilson coefficients. 3 Likewise, one should also take into account that new physics effects in eq. (2.3) affect the extraction of the CKM factor V ud . See refs. [41,42] for the discussion of this issue in the context of extracting V ud from beta decays. For the present purpose, however, V ud appears in COHERENT observables in combination with the pion decay constant f π ± , and the product of the two is usually extracted from the experimental measurements of π → µν decay width.
In the WEFT Lagrangian above, the quarks u, d and charged leptons α are in the basis where their kinetic and mass terms are diagonal, whereas the neutrino fields ν α are taken in the flavor basis. The latter are connected to neutrino mass eigenstates through the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix: ν α = 3 n=1 U αn ν n . As usual, flavor indices are denoted with Greek letters, whereas mass eigenstate indices are denoted with Roman letters.
The Wilson coefficients f f X and ρ X parametrize the effect of new interactions mediated by non-standard heavy particles (heavy compared with the COHERENT physics scale), such as, e.g. leptoquarks or Z and W vector bosons. It is customary to define the Wilson coefficients in eq. (2.3) at the renormalization scale µ = 2 GeV in the MS scheme, since lattice QCD provides hadronic decay constants at that scale and scheme. We make the analogous scale choice for the Wilson coefficients qq αβ in eq. (2.1). Below 2 GeV, the quarklevel Lagrangians will be matched in the next subsection to the nucleon-level EFT. On the other hand, for the Wilson coefficients in the leptonic interactions in eq. (2.4) it is convenient to choose a lower renormalization scale, µ m µ . In any case, in this paper we will take JHEP05(2023)074 into account only the one-loop QCD running of the Wilson coefficients, which can lead to substantial effects, but only concerns the parameters S,P,T in eq. (2.3) [43]. Although the Wilson coefficients are complex in general we will treat them as real throughout this paper for the sake of simplicity.

Nucleon-level EFT
The EFT Lagrangian in the preceding subsection contains quarks fields. However, quarks are not useful degrees of freedom at the energies relevant for CEνNS experiments. In particular, in the COHERENT experiment some of the neutrinos are produced in pion decay, and they subsequently scatter on heavy nuclei. Thus, we need to connect the quark-level formalism to hadronic-and nuclear-level observables. Regarding pions, it is customary to connect the WEFT Lagrangian in eq. (2.3) directly to the decay amplitude, using the matrix element of the quark bilinears between a pion state and the vacuum, see section 3.1. On the other hand, regarding scattering on nuclei, it is convenient to introduce an intermediate step in the form of an effective Lagrangian where the degrees of freedom are nucleons rather than quarks. To match the nucleon-and quark-level Lagrangians we need the matrix elements of the quark operators between the nucleon states. Let us denote the incoming nucleon momentum k, and the outgoing nucleon momentum k . In the near-zero recoil limit, k − k ≈ 0, we have (2.6) Thus, in the limit where the strange content of the nucleon is ignored 4 one finds that Using these matrix elements, we can write down the effective Lagrangian for vector NC interactions between neutrinos and nucleons: where the leading order matching of the Wilson coefficients of the two EFTs reads 4 See ref. [38] for the corrections due to the strange content. Taking this into account in our analysis would introduce a weak dependence of the COHERENT observables on the WEFT Wilson coefficients ss V . However, currently COHERENT is sensitive only to | ss αβ | 1, and for this reason we ignore the strange content of the nucleon.

JHEP05(2023)074
In the SM limit at tree level we have One can see that the NC interactions between neutrinos and vector nucleon currents are approximately protophobic, due to the accidental fact that sin 2 θ W ≈ 1/4. For this reason, low-energy neutrinos scatter mostly on neutrons in nuclei. At leading order, the nucleon weak charge Q N w can be defined simply as the value of the effective couplings in eq. (2.8) at some fixed renormalization scale, [Q N w ] αβ ≡ g νN αβ (µ). One can generalize this definition so that it includes radiative corrections, which cancels the renormalization scale dependence, but it becomes process dependent [44]. In that approach, the SM values for the weak charges are the following Note that with this definition the proton weak charge depends slightly on the neutrino flavor. The difference between the proton weak charge experienced by electron-and muonneutrinos is however numerically irrelevant at present, given the accuracy of the COHER-ENT experiment. Accordingly, we will just neglect such differences and work with the muonic weak charge also in relation with electron neutrinos (see section 3.1). In order to connect the nucleon-level EFT to the nuclear scattering observables, it is convenient to take the non-relativistic limit of the Lagrangian in eq. (2.8), since that will allow us to calculate observables for nuclei of arbitrary spin. At the zero-recoil level it takes the form Above, we traded the relativistic nucleon Dirac fields for the non-relativistic ones denoted as ψ N , which satisfy the Schrödinger equations of motion. We also dropped all terms containing spatial derivatives ∇, which correspond to recoil effects. 5 Now, coherent neutrino scattering amplitudes will involve matrix elements of ψ † N ψ N between nuclear states. For a nucleus N with momentum k, energy E N , charge Z, mass number A, spin J, and spin projection along the z-axis J z , the rotational and isospin symmetry requires the matrix elements to take the form where the primed variables refer the final nuclear state (E N = m N ≈ E N in the target rest frame). The nucleon form factors F p,n (q 2 ) are equal to 1 at q 2 ≡ (k − k ) 2 = 0 due to JHEP05(2023)074 isospin symmetry, although we will not take that limit in our studies. The factor A − Z in the neutron matrix element is at the origin of the coherent enhancement of the neutrino scattering on heavy nuclei.

COHERENT event rate
Let us consider neutrinos produced by a source S through the process S → X α ν k , where X α is a one-or more-body final state that contains a charged lepton α = e, µ, τ , and ν k is a neutrino-mass eigenstate (k = 1, 2, 3). These neutrinos propagate a distance L -conserving its mass index k -and are detected via the process ν k N → ν j N , where j = 1, 2, 3 is again a mass index and N denotes the target nucleus. Let us consider the differential number of detected events per time t, incident neutrino energy E ν , nuclear recoil energy T and target particle where N T stands for the number of target particles. Previous calculations of the CEνNS rate [1,9,10,38] have been carried out within the SM or under the assumption that the neutrino production is unaffected by New Physics (NP) and that one can thus simply calculate the rate as a flux times cross-section, i.e., where ν α,β are neutrino flavor eigenstates. In this work we present a derivation that is more general in the treatment of new physics contributions. We calculate the event rate in terms of the WEFT Wilson coefficients introduced in section 2. The interactions of neutrinos relevant for their production and detection are assumed to be most general at the leading order in the WEFT expansion, i.e., we allow for the simultaneous presence of all interactions described in eqs. (2.1), (2.3) and (2.4). We allow for arbitrary flavor mixing, both via the PMNS matrix, as well as via the various WEFT Wilson coefficients of 4-fermion interactions. Finally, we take into account that NP affecting neutrino production also affects the muon and pion decay widths, which are used in the COHERENT analysis and to determine some of the input observables (such as G F ).
For this derivation, we need to connect the observable event rate in eq. (3.1) with the production and detection QFT amplitudes, denoted by M P αk ≡ M (S → X α ν k ) and M D jk ≡ M (ν k N → ν j N ), which encode the fundamental physics taking place at production and detection. This connection was obtained in ref. [15] using a QFT approach for the case of charged-current interactions both at production and detection. The main idea behind such derivation was, instead of considering the neutrino production and detection separately, to treat both interactions as a single process [46]. In our CC-NC configuration, that translates into the following process where the ν k neutrino is considered just as an intermediate particle in the amplitude. We can adapt the result found in ref. [15] to describe the CEνNS rate observed at COHERENT

JHEP05(2023)074
taking into account the various differences. In NC neutrino scattering we have no information about the neutrino final mass eigenstate (or flavor) and hence, we should sum over the corresponding mass index j. Secondly, the time variation of the number of source particles N S cannot be neglected in this case (in fact, it produces a time-dependent signal that is measured). Furthermore, at COHERENT one measures the differential number of events per recoil energy T , and thus we will not integrate over that detection kinematic variable.
Finally, the incident neutrino energy is not observed in COHERENT, which can be taken into account trivially integrating R S α over that variable. We note that the assumption of neutrinos emitted isotropically from a source at rest applies to the case of COHERENT.
All in all, the event rate for a source S is given by where complex conjugation is denoted with a bar, m S,N are the masses of the source and target particles respectively and ∆m 2 kl ≡ m 2 k − m 2 l is the mass squared difference between neutrino (mass) eigenstates, which appears in the formula through the usual e −i L∆m 2 kl /(2Eν ) oscillatory factor, which we can simply approximate as one for the case of COHERENT given its very short baseline. The phase space elements for the production and detection processes, dΠ P and dΠ D , are defined as usual: where k i and E i are the 4-momenta and energies of the final states and p n is the total 4-momentum of the initial state. However, in order to obtain the observable of interest at COHERENT, we use a slight modification of the standard production and detection phase spaces denoted by primed subindices and defined by dΠ P ≡ dΠ P dE ν and dΠ D ≡ dΠ D dT , such that R S α then provides the differential number of events per incident neutrino energy E ν and recoil energy T via eq. (3.1). The integral sign involves both integration as well as sum and averaging over all unobserved degrees of freedom such as spin. Finally N S (t) = n POT r S/p τ S g S (t) is the time-dependent number of source particles S, where n POT is the total number of protons on target delivered at the Spallation Neutron Source (SNS), r S/p is the number of S particles produced per proton, τ S is the S lifetime and g S (t) encodes the N S time dependence (normalized to one over each bunch cycle characteristic of the proton beam at the SNS). The precise form of g S (t) depends on several factors such as the S lifetime, time-dependent efficiency and the time-structure of the S pulse.
For antineutrinos the event rate is defined equivalently.

Production and detection amplitudes
Let us first discuss the detection process. The amplitude for neutrino scattering on nuclei as a function of the nucleon-level EFT parameters in eq. (2.8) reads where u ν k,j are the Dirac wave functions of the incoming and outgoing neutrinos. Taking into account the current uncertainties affecting the COHERENT event rate, it is convenient JHEP05(2023)074 to approximate neutron and proton form factors to be equal, F p (q 2 ) = F n (q 2 ) ≡ F (q 2 ), which allows us to write the detection amplitude in the more compact form: where the (dimensionless and Hermitian) nuclear weak charge is defined as To work instead with different neutron and proton form factors would entail working with q 2 -dependent weak charges (or, equivalently, with the g νN coefficients instead of the weak charges), which would make our subsequent discussion and intermediate results more cumbersome. Additionally, the final results for the quark-level WEFT Wilson coefficients are not expected to be affected by this approximation given current uncertainties. On the other hand, the q 2 dependence in the form factor can not be neglected in the studied recoil energy range. To describe it, we will make use of the Helm parametrization [47] (see appendix B for further details). Other recoil effects are numerically unimportant given the current experimental precision.
On the production side we have the 2-body leptonic pion decay π + → + α ν k (α = e, µ) and the 3-body muon decay µ + →ν m e + ν k . Their amplitudes are given by 6 where v α ,ū ν k are the Dirac spinor wave functions of the charged lepton and the neutrino, respectively, and the pion decay constant f π ± is defined by 0|dγ µ γ 5 u(0)|π + (p) = ip µ f π ± . Above we introduced the shorthand notation The transposition in the P L,R matrices in eq. (3.7) is defined such that it only affects the two neutrino indices when applied to the ρ L,R Wilson coefficients.

Amplitudes product and phase space integrations
On the detection side, for neutrinos we find In the muon decay amplitude M P,µ mk we omit the charged lepton subindex (α = e in this case) and we include both the neutrino and antineutrino mass eigenstate indices. This allows us to use the same notation for the neutrino and antineutrino rates.

JHEP05(2023)074
where T = E N − m N is the (kinematic) nuclear recoil energy, and thus q 2 = 2M N T . The same result holds for antineutrinos except for the ordering of the (l, k) indices on the right-hand side.
On the production side, for pion-decay neutrinos we obtain where E ν,π = (m 2 π ± − m 2 µ )/(2m π ± ) is the energy of the neutrino emitted in the 2-body pion decay. On the other hand, for muon-decay neutrinos and muon-decay antineutrinos we find, respectively, the following results 7 We have neglected O(m e /m µ ) corrections, which include the crossed RL and LR terms.
We would like to write these results in terms of the observable pion and muon decay widths, since COHERENT uses those quantities in their flux predictions. They can be obtained from the expressions above integrating over the (anti)neutrino energy and working with equal neutrino indices in M P andM P , which are summed over. This gives Tr P L P † L + P R P † R , (3.12) up to radiative and m e /m µ corrections. The muon decay width can be expressed as Γ µ→eνν = m 5 µ /(384π 3 v 4 ), which represents the definition of v (equivalently, of G F ≡ ( √ 2v 2 ) −1 ), and remains valid in the presence of new physics. 8 In other words, in our input scheme the possible new physics contamination in the determination of the Fermi constant is absorbed into the parameter v, whose value is fixed by experiment. These NP effects are also absorbed in the Wilson coefficients qq V,A and ud X , since v is used in their definitions, cf. eqs. (2.1) and (2.3). This can be seen explicitly matching to the Warsaw-basis SMEFT, cf. section 6 and appendix C. 7 For the sake of clarity we have written explicitly the sum over the (anti)neutrino of mass m k instead of considering it implicit inside the integral sign. 8 It is straightforward to see this implies v 2 µ mµ , which is consistent with the discussion after eq. (2.4).

JHEP05(2023)074
Using these results we can rewrite eq. (3.10) and eq. (3.11) as The NP effects that appear in the numerators (involving the PMNS matrix) are those contributing directly to the production amplitudes in eq. (3.10) and eq. (3.11), whereas those in the denominators enter indirectly because they affect the pion and muon decay widths. We will refer to these contributions as direct and indirect NP effects respectively. 9 We stress that both contributions appear at the same order and are generated by the same EFT operators, so it is not consistent to include only the direct piece.

Event rate
The total rate per recoil energy T and time t detected at COHERENT is given by: where R π µ and R µ e (R µ e ) denote the event rates mediated by neutrinos produced in pion decay and by (anti)neutrinos produced in muon decays respectively. These three event rates are obtained plugging the results of eq. (3.13) in eq. (3.3). The integral over the neutrino energy is trivial in the pion-decay case since the neutrino energy is fixed, whereas for muon decay the lower integration limit is E min ν (T ) = T 2 1 + 1 + 2 m N T (i.e. the minimum energy required to produce CEνNS with a recoil energy T ) and the upper one is simply m µ /2. Working in the L = 0 limit and separating the events in prompt (i.e. produced in pion decays) and delayed (i.e. produced in muon decays), we can rewrite eq. (3.14) as where g π,µ (t) encode the N π,µ time dependences (normalized to one over each bunch cycle), as discussed at the beginning of section 3. The prompt and delayed components are given by

JHEP05(2023)074
The generalized squared chargesQ 2 f are defined as the following positive and targetdependent quantities 10 The explicit form of the f f (T ) functions, which is not very enlightening, can be found in eq. (A.1). Finally the number of (anti)neutrinos produced per proton via pion (muon) decay, f π(µ) ν/p , are given by Finally, the allowed values for the recoil energy are for the delayed ones.
The prompt and delayed event rates in eq. (3.16) represent one of the main results of this work, which is thus worth analyzing in some detail. First, let us note that the PMNS factors are not present anymore (they were removed using the unitarity condition U U † = U † U = 1), which means that COHERENT is not sensitive to the PMNS mixing angles and phases, as expected in an L ≈ 0 experiment. Secondly, let us note that expressions for the rates in eq. (3.16) are equal to the SM expressions except for the fact that the nuclear weak charge has been replaced by a generalized weak charge that (i) is different for muon neutrino, electron neutrino and muon antineutrino; and (ii) contains non-standard effects affecting detection and production. To put this in more explicit terms, we can re-write the prompt and delayed event rates as follows: where the fluxes are the usual ones and the cross sectionsσ f are defined in the usual form but using the generalized chargesQ f , cf. appendix A. Even though we started from very general premises, the final result is similar to the SM formulas and to the usual NSI expressions (with NP only in the detection side), except for the introduction of the generalized weak charges. This unexpected result makes 10 Let us note that the RR term inQ 2 µ (Q 2 e ) is generated by non-standard effects in ν-mediated (νmediated) events. Thus one should only identify the fμQ 2 µ (feQ 2 e ) terms in the delayed event rate with ν-mediated (ν-mediated) events if the RR contributions are zero. The same caveat holds for the flavor indicesμ and e used in those two terms, which only refer to the flavor of the mediating (anti)neutrino in the case of flavor-diagonal interactions in production (and no RR terms).

JHEP05(2023)074
the phenomenological analysis very simple, since it represents a simple modification with respect to the standard approach in the previous literature. The conceptual change is however much deeper and one should keep in mind that in our general analysis the generalized weak charges contain non-standard lepton flavor violating effects affecting neutrino production. For instance, our general expression includes possible contributions through the process π → µ ν τ , despite not having introduced a ν τ flux in eq. (3.19). Thus one should keep in mind that the event rate in eq. (3.19) is just a practical parametrization, but the factorization in fluxes and cross section, as well as the subindices ν µ , ν e andν µ , do not have physical meaning except in the SM case and some simple BSM scenarios.
In general, it is not possible to carry out a naive factorization of the event rate in eq. (3.16) in fluxes and cross sections, since there is a matrix multiplication between production and detection quantities in the generalized squared chargesQ 2 f in eq. (3.17). For simplicity let us consider the case of pion decay production and CEνNS detection, where one can easily see that Before discussing some interesting specific cases, let us mention briefly how the analysis is modified if we consider different q 2 -dependent form factors for neutron and proton, i.e., F n (q 2 ) = F p (q 2 ) = const. In that case, theQ 2 f parameters are not convenient objects to summarize experimental results, because they become q 2 dependent. In the COHER-ENT rate expression, the product of the weak charge matrix squared and the form factor squared, Q 2 (F (q 2 )) 2 , would be replaced by the matrices (g νp ) 2 , (g νn ) 2 , g νp g νn , and g νn g νp , accompanied by the appropriate powers of the ZF p (q 2 ) and (A−Z)F n (q 2 ) functions. Thus, we would go from 3 parameters per target (Q 2 f =e,µ,μ ) to 12 target-independent parameters. They are reduced to 9 parameters if the (production and detection) NP parameters are real, because the generalized squared charges obtained "replacing" Q 2 with g νp g νn and g νn g νp are equal.

Interesting limits
SM limit. If all NP effects are switched off we recover the SM prediction, with a single nuclear weak charge where (3.21) and the nucleon weak charges are given in eq. (2.11). Note that in principle the muon and electron weak charges have slightly different values, however the difference is irrelevant given the COHERENT accuracy, cf. eq. (2.11). In our analysis we will take the muon weak charge as the reference value. The SM scenario has of course been thoroughly studied before [8,44]. For each target nucleus there is a single quantity, Q SM , to be extracted from experiment, which is predicted in the SM in function of the weak mixing angle. Thus, in the SM limit coherent elastic neutrino-nucleus scattering can be regarded as a probe of the weak mixing angle, see e.g. [20,33,34,[48][49][50]. One should remark however that

JHEP05(2023)074
other probes, such as Z-pole physics, atomic parity violation, or parity-violating electron scattering have currently a much better sensitivity to the weak mixing angle.
New physics in production. We move to the case where new physics affects the CO-HERENT observables via neutrino production in pion and muon decay. The matrices P and P L,R , which encode these effects, are allowed to be completely generic. On the other hand, we assume here that we can ignore new physics in detection. This implies that the weak charge Q defined in eq. (3.6) is proportional to the unit matrix and thus it commutes with P and P L,R . It then follows from eq. (3.17) that the new physics production effects completely cancel out in the generalized weak charges, and we recover the SM limit in eq. (3.20) with a single nuclear weak charge. All in all, COHERENT data are completely insensitive to new physics affecting only the CC semileptonic and leptonic WEFT operators due to cancellations between direct and indirect new physics effects. This observation invalidates the bounds found in ref. [51], where the indirect NP effect was not taken into account. A more intuitive way of understanding this null sensitivity is the following. The CC operators in eq. (2.3) certainly affect the pion decay rate to muon and neutrino, but they do not distort the kinematics. Their effect has been fully absorbed into the experimental value of BR(π → µν), which is used to calculate the neutrino flux. Moreover, in the particular case at hand BR ≈ 1, that is, ∼100% of the pions will decay to muon and neutrino for any reasonable values of [ L,R,P ] αβ . Similarly, new leptonic operators in eq. (2.4) affect the muon decay rate, but their effect has been fully absorbed into the experimental value of the Fermi constant.
Our work is the first one that takes into account the direct and indirect effects in production, as well as the possible cancellations. Let us stress however that new physics in production cannot be ignored completely: its effects do not cancel out if there are accompanying new physics effects in detection.
New physics in detection. If we neglect NP in production our expressions reduce to those found previously in the NSI literature [9,10,48], where two free parameters are present (instead of three). Namelỹ where Q SM is given in eq. (3.21), and That is, for a given target, COHERENT is linearly sensitive only to two linear combinations of the four WEFT Wilson coefficients: uu µµ , dd µµ , uu ee , and dd ee , which describe flavor-diagonal NC interactions between neutrinos and quarks. 11

Experimental input
The COHERENT collaboration uses a series of detectors to detect neutrinos produced by the Spallation Neutron Source (Oak Ridge National Laboratory) through CEνNS. At this facility, high-energy protons with E ∼ 1 GeV hit a mercury target to produce π + and π − . The latter are absorbed, whereas the positive pions decay at rest into the prompt neutrinos and positive muons. The latter correspond to the source particles of the delayed (anti)neutrinos.
We will analyze the two available measurements of the CEνNS interaction performed by this experiment: one performed on a liquid argon target (LAr) [35] and another one using a target consisting of a mixture of cesium and iodine (CsI) [36]. The input needed to calculate the number of prompt and delayed events for these two measurements through eq. (3.16) is summarized in table 1. The number of target particles N T is obtained as the ratio of the active mass of the detector m det and the mass of the interacting nuclei. For the CsI measurement, we treat cesium and iodine as a single nucleus with (Z, A) = (54, 130). This will allow us to analyze CsI data in terms of only 3 charges (instead of 6) and we do not expect it to have any impact in the final bounds on the WEFT Wilson coefficients, since the atomic numbers for Cs and I are very similar (namely Z Cs = 55, Z I = 53).
Naively, one only has to integrate the expression in eq. (3.16) over the recoil energy T in each bin to obtain the expected number of prompt/delayed events (for a given value of the generalized weak chargesQ f ), that is This statement depends on the definition of the WEFT coefficients and on the input scheme. In particular, if one uses v0 (instead of v) in eq. (2.1), then we would find that COHERENT is also linearly sensitive to the CC interaction [ρL]eµµe because of its effect on the muon decay, and hence on GF (or v), which is used to calculate the CEνNS cross section. In our approach such effects have been absorbed inside the NC coefficients qq ll . Both approaches are of course equivalent, as can be seen explicitly when they are matched to the Warsaw-basis SMEFT, cf. section 6 and appendix C. Last, we note that this caveat also applies to the previous discussion about NP in production.

JHEP05(2023)074
Parameter CsI [36] LAr [35] (Z, N ) ( Table 1. Input used to calculate the differential number of prompt and delayed events dN/dT for each setup, taken directly from the corresponding experimental article. Following those references we approximate the number of pions and muons per proton to be equal However, this simple step needs to be modified to take into account various experimental effects. The first thing to consider is that COHERENT does not measure its events in nuclear recoil energy (T ), but in electron-equivalent recoil energy (T ee ). These two magnitudes are related as follows where QF(T ) is the so-called quenching factor, which depends on the detector and the recoil energy T , as we indicated explicitly. Moreover, one has to introduce an energy resolution function R(T rec ee , T ee ), which relates the true value of the electron-equivalent recoil energy, T ee , with the reconstructed one, T rec ee , that is registered at the detector. Finally, the efficiency of the detector, (T rec ee ), should also be taken into account. These considerations are collected in the following modified expression for the number of prompt and delayed events in the i-th T rec ee bin [36,52] where a = prompt/delayed. The theoretical prediction for the (differential) number of events, dN a /dT , is given in terms of the nuclear recoil energy T , as provided in eq. (3.16). Note that we have expressed the energy resolution R(T rec ee , T ee ) in terms of T instead of T ee by means of the QF. That way, the integral over T allows us to go from the nuclear recoil energy T to the reconstructed electron-equivalent recoil energy T rec ee . The integration limits, the quenching factor and the efficiency and energy resolution functions that we use for each dataset are discussed in appendix B.
Additionally, the CsI analysis presents the data in terms of the number of photoelectrons (PE) that are recorded for each event instead of using the electron-equivalent recoil energy. These two magnitudes are simply connected by PE = LY × T rec ee , where the light yield LY is the number of PE produced by an electron recoil of one keV. Thus, the general expression in eq. (4.3) holds also for CsI with the replacement T rec ee → PE. Above we presented the expression for the prompt and delayed events, which have to be summed to provide the observed number of CEνNS events. Equivalently this result can JHEP05(2023)074 be obtained integrating t in eq. (3.15) over the entire bunch cycle and taking into account that the g a (t) functions are normalized to one. If instead one is interested in the double distribution in nuclear recoil energy T and time t, then we will have where we introduced a second index j that refers to the j-th time bin. The g a j factors are the prompt/delayed probability distributions for the timing of the events (calculated integrating the g a (t) functions over the j-th time bin), which can be extracted from the COHERENT publications, cf. appendix B.
Finally, one also has to include background events and nuisance parameters to parametrize the most relevant uncertainty sources. Thus, the predicted number of events has the following generic form where we have indicated explicitly the dependence of the expected number of CEνNS events on the generalized squared weak charges for the nucleus N , denoted by The expected number of background events of type a, denoted by N bkg,a ij , is obtained following COHERENT prescriptions, as described in detail in appendix B. The typical background sources are the steady state (SS) background, the neutrino-induced neutron (NIN) background and the beam-related neutron (BRN) background, although the way each of them is characterized differs slightly in every measurement. Finally, the h signal/bkg,a ij ( x) functions are linear in the nuisance parameters x and vanish at their central values x = 0. The specific form of these functions for each experimental analysis is given in appendix B following once again the COHERENT prescription.
In our numerical analysis we use the 2D distributions (in time and recoil) measured in the CsI and Ar works [35,36]. For each of these two datasets we work with a Poissonian chi-squared function with the following generic form where σ n is the uncertainty of the x n nuisance parameter. All in all we have 52x12 bins in CsI and 4x10 bins in LAr, cf. appendix B for further details.

Generalized nuclear weak charges
In this section we present the results of the analysis of LAr and CsI data in terms of the generalized nucleus-dependent weak chargesQ f . Since the event rate depends quadratically on these charges, it is convenient to work with their squared valuesQ 2 f .

Argon charges
We carry out a 2D fit to the nuclear recoil energy and time distributions, as described in the previous section. This fit has 40 experimental inputs(with their associated uncertainties and backgrounds), 9 nuisance parameters (with their uncertainties) and the threẽ Q 2 f charges that contain the UV information. We find that the distribution of these three charges is approximately Gaussian, with the following results: 12 along with the nine nuisance parameters that we do not display. We have normalized the results using the SM value, Q 2 SM,Ar = 461, with an associated small error that can be neglected for the purpose of this work. The results in eq. (5.1) are in perfect agreement with the SM predictionQ 2 f /Q 2 SM = 1. We can disentangle the first charge,Q 2 µ , from the other two thanks to its different time dependence: the former enters the rate via (prompt) pion decay, while the latter do it via (delayed) muon decay. On the other hand the recoil energy distribution only allows for a mild separation ofQ 2 e andQ 2 µ . To get rid of the large correlations, which obscure the strength of the results, let us rewrite them as the following uncorrelated bounds where we have highlighted the most stringent constraints (note that the SM prediction is one by construction). A particularly interesting case is the SM supplemented by the following contributions:Q 2 µ =Q 2 µ = Q 2 µµ andQ 2 e = Q 2 ee , which is the most general setup that we can have when considering NP only at detection or in a linear analysis (see section 3.4). In this case we find:

JHEP05(2023)074
The results of this 2D fit are shown in the left panel of figure 1, where we also present the allowed regions if one only uses the total number of events, the energy distribution or the time distribution.
Finally in the SM scenario there is only one weak charge (Q 2 µ =Q 2 µ =Q 2 e ≡ Q 2 ), and we find Q 2 /Q 2 SM Ar = 1.03 ± 0.45.

CsI charges
We have carried out a fit to the 2D distributions (nuclear recoil energy and time) provided in the CsI analysis. This fit has 624 experimental inputs (with their associated uncertainties and backgrounds), 4 nuisance parameters (with their uncertainties) and the three CsI generalized weak charges. Once again we find that the distribution of theQ 2 f charges is approximately Gaussian, with the following results: along with the nuisance parameters. Here, Q 2 SM,CsI = 5572. As in the LAr case, we can separate much betterQ 2 µ from the other two charges thanks to the use of the time information. The results can be rewritten as the following uncorrelated bounds: where we find again good agreement with the SM predictions (one).
Considering only NP at detection we find which can be re-written as the following uncorrelated bounds: The results of this 2D fit are shown in the right panel of figure 1. We present as well the allowed regions if one uses only the total number of events, the energy distribution or the time distribution. Finally, we also show the result obtained using the full 2D distribution of the first CsI dataset [3], which is in good agreement with figure 6 in ref. [29]. One observes a clear improvement when the entire CsI dataset is used. Finally in the SM scenario there is only one weak charge (Q 2 µ =Q 2 µ =Q 2 e ≡ Q 2 ), and we find Q 2 /Q 2 SM CsI = 1.04 ± 0.16. For each dataset we include information from the total number of events (N), from the recoil energy distribution (E), from the timing distribution (t) and from the 2D distribution (E+t). In the right panel we also show the contour obtained using only the first CsI dataset [3]. The shaded area indicates the unphysical region (negative squared charges).

WEFT coefficients
In this section we move to consider the constraints on the nucleon-and quark-level EFT Wilson coefficients, stemming from our analysis of LAr and CsI CEνNS data.

Linear BSM expansion
As shown in eq. (3.23), at linear order in New Physics there are only 2 free parameters per target: the flavor diagonal muon and electron weak charges, [Q 2 ] N µµ,ee . Using eq. (3.24) we can express our bounds on the four weak charges (with N = Ar, CsI, see eqs. (5.3) and (5.7)) as bounds on the four nucleon-level EFT Wilson coefficients δg νp ee,µµ and δg νn ee,µµ . We find that we can constrain the following orthogonal and uncorrelated linear combinations of couplings: (5.9) At the quark level, using the map in eq. (2.9), we can translate these results into constraints on the following orthogonal and uncorrelated linear combinations of WEFT

JHEP05(2023)074
Wilson coefficients: The last two constraints in eq. (5.10) are not expected to change with the inclusion of quadratic corrections 13 and they represent another central result of our work. Once again, their errors are Gaussian to a good approximation and the application of these EFT constraints to more specific setups is straightforward. We stress that we did not neglect NP affecting production, as usually done in the past literature. Instead, we showed that, with our Lagrangian input choice, they are absent at this order in the EFT expansion.
It is interesting to discuss the results in a more constrained WEFT scenario when lepton-flavor universality of the relevant Wilson coefficients is assumed: uu µµ = uu ee ≡ u and dd µµ = dd ee ≡ d . Then the 4-parameter fit in eq. (5.10) reduces to the two-parameter one: The results above lead to weak and effectively meaningless marginalized constraints on δg νN and qq because only two of the constraints have uncertainties O(1), where the linear expansion makes sense. The situation changes when only one operator is present, which leads to stringent and reliable individual limits: (5.13)

Non-linear effects
In this section we take into account the full expression for the COHERENT event rate, including non-linear terms in the new physics Wilson coefficients, and discuss how this changes the results compared to the linear analysis. For simplicity, we start by discussing the impact of non-linear effects in analyses where only one operator at a time is present. As expected, the linear bounds presented in eq. (5.13) for the neutral-current Wilson coefficients, qq αα , are barely affected by non-linear terms, cf.  = β), which only enter the event rate at quadratic order in NP. The corresponding results, obtained putting one operator at a time, can be found in table 2. 14 They are also compared with the flavor-diagonal ones in figure 2. Let us stress that, with our definition of the WEFT coefficients and input parameters, the effect of charged-current operators cancels in the rate (at all orders) if only one operator at a time is considered, as discussed in section 3.4.

JHEP05(2023)074
Let us now discuss some cases where more than one operator is present at the same time. In figures 3 and 4 we consider scenarios with two free NP parameters, with the remaining ones set to zero. A quick glance reveals that CsI data drive the constraining power of the fit in all cases.
Comparing the upper two panels, we can see that the fit with only electron couplings yields noticeably weaker constraints than the one with muon parameters. This is to be expected since the contribution to the event rate coming from the muonic neutrinos is larger than the one coming from the electronic ones. We see that the CsI data is precise enough to separate the allowed region in two bands: one compatible with the SM and a second one corresponding to the above-mentioned dark solution. The well-known blind directions that these panels display is due to the potential cancellations among the linear combinations of up-and down-quark couplings in the weak charges, namely (A + Z) uu V + (2A − Z) dd V = constant, cf. eq. (3.22). Since the blind directions are almost parallel for CsI and Ar, the combined dataset also shows this feature.
14 In this work we have assumed that all Wilson coefficients are real. However, our one-at-a-time bounds on lepton-flavor off-diagonal coefficients are trivially generalized to bounds on their modulus squared if they are complex, since they do not interfere with the SM contributions. The second row in figure 3 studies the cases where neutrinos are coupled only to down quarks and only to up quarks. In these fits we have NP in both the electron and muon charges and thus there are four solutions, corresponding to ([Q] ee = ±Q SM , [Q] µµ = ±Q SM ). Current COHERENT data is only able to separate the upper two dark solutions, whereas a third dark solution remains connected with the SM one. Adding CEνNS measurements at reactors isolates the SM solution [32]. The scenarios that we have discussed above have been thoroughly discussed in previous works using different COHERENT datasets and projections [19,20,23,26,27,[29][30][31][32][33]. Our work and the recent work of ref. [34] are the first ones to use the entire COHERENT dataset presently available (2D distribution in LAr and CsI) [35,36] to constrain NSI coefficients, and thus represent the current state of the art. Since the CsI detector has been decommissioned, they represent the final results with that target. Previous works have studied these NP scenarios using the SM fluxes and modified cross sections. Our more general approach reduces to this "factorized" NSI description if one neglects NP effects in production, as discussed in section 3.4. Our numerical results for the separated LAr and CsI analyses agree well with previous works, including COHERENT analyses [35,36]. Our combined bounds (LAr+CsI) are also in good agreement with the recent results in ref. [34].

JHEP05(2023)074
Let us now compare our COHERENT results with other NSI probes, which were reviewed and compiled in ref. [53]. For simplicity we focus on bounds obtained switching on one operator at a time. For the muonic couplings, dd µµ and uu µµ , the bounds obtained from JHEP05(2023)074 COHERENT data match the best existing constraints, which come from atmospheric and accelerator neutrino data [54]. For the electronic couplings, dd ee and uu ee , our COHERENT results are much stronger than the limits extracted from CHARM data [55] and comparable to those obtained from Dresden-II reactor data [32,56]. For the flavor violating NSIs, our results for qq µτ are roughly 20× weaker than those obtained from IceCube [57]. Finally, our one-at-a-time bounds on qq eµ coefficients from COHERENT data are roughly two times weaker than those obtained in a global fit to oscillation data [29], whereas for the qq eτ coefficients they are similar. The relatively weak sensitivity from our analysis to off-diagonal NSIs is to be expected since oscillation observables are linearly sensitive to them, whereas CEνNS is only quadratically sensitive. On the other hand, CEνNS is best suited to study flavor-conserving NSIs, with interesting synergies observed in combined analyses with oscillation data [29].

Production and detection effects together
Our general approach allows us to go beyond the well-known cases discussed above, and study situations where NP effects are present both in production and detection.
For instance, we can study a setup where the NC coefficient dd µτ is accompanied by the CC semileptonic Wilson coefficient ud µτ ≡ [ ud L ] µτ . The latter affects neutrino production (it generates π → µ ν τ ), whereas the former affects the detection of muon and tau (anti)neutrinos (it generates ν µ N → ν τ N , ν τ N → ν µ N and likewise for antineutrinos). Figure 5 (left panel) shows the allowed regions when both parameters are present at the same time. Let us stress once again that in our formalism these bounds are obtained without introducing a ν τ flux. As explained in section 3.3, one would obtain different (and thus incorrect) results if one calculates the event rate in a flux-times-cross-section factorized form. Namely, one would lose all sensitivity to the ud µτ parameter. The study of simultaneous NP effects in neutrino production and detection is particularly relevant in setups with explicit electroweak symmetry, since neutrinos and charged leptons form SU(2) L gauge doublets. As a result, non-standard contributions to ν α N → ν β N come in general with non-standard effects in leptonic pion decay, π → ν α β . Let us consider for instance the SMEFT operators [O (along with their conjugates so that the Lagrangian is Hermitian), and let us abbreviate their associated Wilson coefficients as c 3 ≡ v 2 [C  In these fits only two couplings are allowed to be present at a time. Finally, we also include the constraint from the ratio Γ(π → eν)/Γ(π → µν) [37,58,59]. COHERENT: We show in figure 5 (right panel) the bounds that we obtained on the coefficients c 1 and c 3 using COHERENT data. We can also constrain charged-current NSIs using the measurements of leptonic pion decay widths. To make things simpler, we work with the ratio R π ≡ Γ(π → eν)/Γ(π → µν) that, in the specific cases discussed above, is modified as R π = R SM π /(1 + ( ud µτ ) 2 ). Figure 5 shows the interplay between this constraint and the one obtained from COHERENT data.

Comparison and combination with other precision observables
In this section we discuss the place of the COHERENT experiment in the larger landscape of electroweak precision observables. To this end we will employ the SMEFT framework [16,17], which will allow us to combine information from COHERENT and other experiments below the Z-pole, with that obtained by the high-energy colliders at or above the Zpole. 15 It will also allow us to combine the information from NC and CC processes, which are related by the SU(2) L gauge symmetry. We consider operators up to dimension six, using the standard SMEFT power counting where the corresponding Wilson coefficients are O(Λ −2 ) in the new physics scale Λ. Consequently, we expand observables to order 1/Λ 2 , ignoring higher order corrections. This implies that we only need the linearized version of the COHERENT results obtained in the WEFT formalism, cf. eq. (5.10).
The COHERENT experiment probes contact 4-fermion interactions between the lefthanded lepton doublets l L and quark doublets q L , and singlets u R , d R . The relevant JHEP05(2023)074 dimension-6 operators are [17] L SMEFT ⊃ C (1) Here, the capital letter Wilson coefficients are dimensionful, [C X ] = mass −2 . For the numerical analysis it is more convenient to work with dimensionless objects c X = v 2 C X , where v 246.22 GeV. The SM fields are 3-component vectors in the generation space, however the flavor index is suppressed here to reduce clutter. In this section we assume that the Wilson coefficients are flavor universal, more precisely, that they respect the U(3) 5 flavor symmetry acting on the three generations of q L , u R , d R , l L , e R . This is by far the most studied SMEFT setup, especially in the context of global fits, and often described simply as the EWPO fit. As we will show below, even in this restricted framework, COHERENT has a significant impact on the global fit. Later in appendix C we will relax this assumption, in which case the impact of COHERENT will be even more spectacular thanks to lifting degeneracies in the multi-dimensional parameter space of Wilson coefficients.
From the SMEFT point of view, the COHERENT experiment also probes the coupling strength of the Z boson to quarks and neutrinos. For these interactions we will use the Higgs basis parametrization (see [63] for a recent summary): where g L and g Y are the gauge couplings of the SU(2) L × U(1) Y local symmetry, . Above, the effects of the dimension-6 SMEFT operators on the Z couplings to quarks are parametrized by the four dimensionless vertex corrections δg Zu/d L/R . There is one more parameter δg Zν L describing the dimension-6 effects on the Z coupling to neutrinos. In the Higgs basis it is expressed by other leptonic vertex corrections: δg Zν L = δg W L + δg Ze L . The COHERENT results analyzed in this paper constrain the 4-fermion Wilson coefficients in eq. (6.1) and the vertex corrections in eq. (6.2). These Wilson coefficients are related to the NC WEFT Wilson coefficients in eq. (2.1) by [40] uu where there is no implicit sum over the repeated index α. Because of our assumption of U(3) 5 symmetry, the expression is the same for any value of the index α, that is to say, the quarks interact with the same strength with all flavors of the neutrino. Therefore, it  In the reminder of this section we will compare the strength of the COHERENT constraints on SMEFT coefficients to that of the other electroweak precision measurements. Ref. [40] compiled the input from experiments sensitive, much as COHERENT, to flavor-conserving Zf f vertex corrections and 4-fermion operators with two leptons and two quarks. That analysis included Z and W pole measurements, e + e − → ff data, (noncoherent) neutrino scattering on nucleon targets, atomic parity violation (APV), parityviolating electron scattering, and the decay of pions, neutrons, nuclei and tau leptons. Ref. [40] also included purely leptonic observables since, in addition to their dependence on SMEFT 4-lepton operators, they are sensitive to some of the vertex corrections in eq. (6.2). We use the likelihood for the SMEFT Wilson coefficients constructed in ref. [40], updated to include new theoretical and experimental developments (the full list of observables used and the updates can be found in appendix C). In the following we compare and combine these constraints with the ones obtained in this paper using the COHERENT data.
The first comment is that the impact of COHERENT is negligible if only a single Wilson coefficient appearing in eq. (6.3) is present at a time. In table 3 we show the uncertainty obtained in such a one-at-a-time fit. We can see that the sensitivity of COHERENT is inferior by 1-2 orders of magnitude compared to that achieved by a combination of other electroweak precision measurements. This is not surprising, given that the latter contain a number of observables that have been measured with a (sub)permille precision (namely LEP1, APV, or baryon decays), while COHERENT currently offers a percent level precision.
However, most new physics models generate several operators simultaneously, and thus a global analysis is required to assess the importance of COHERENT data. With this in mind, we turn to analyzing the situation when all flavor-universal dimension-6 SMEFT Wilson coefficients are allowed to be present with arbitrary magnitudes within the regime of validity. Now we are dealing with a multi-dimensional parameter space, where certain directions may not be constrained by the most precise observables, and where the input from COHERENT may be valuable. More precisely, in the flavor universal case the JHEP05(2023)074 observables taken into account in our analysis probe 18 independent Wilson coefficients. In addition to the six defined in eqs. (6.2) and (6.2), our analysis probes 11 more four-fermion operators as well as the vertex correction to the Z boson coupling to right-handed leptons. For their definition see appendix C, in particular eqs. (C.2), (C.3) and (C.7). We find the fully marginalized constraints We highlighted the constraints on c lu , c ld , and c eq , which improve significantly, by about 30-40%, after including the COHERENT input. The improvement is visualized in the left panel of figure 6. for the fit including COHERENT is (6.6) In the U(3) 5 -invariant SMEFT, one can obtain bounds on 10 new (combinations of) Wilson coefficients from diboson production at LEP2 and Higgs measurements [64]. The rest of U(3) 5 -invariant SMEFT coefficients, which are not probed by these 2 global fits, are made up of only Higgs doublets, only quarks, or only gluons, or violate CP.
Another way to illustrate the impact of COHERENT is to consider global constraints on the combinations u and d of the SMEFT Wilson coefficients, defined in eq. (6.3). Let us recall that COHERENT alone constrains one linear combination at a percent level: 0.67 u + 0.74 d = 0.003 (10), while the orthogonal combination is poorly constrained. This is shown in the right panel of figure 6 as the diagonal beige band (this is simply a zoomed JHEP05(2023)074 in version of figure 5). We can use the results of our flavor-diagonal global SMEFT fit of electroweak precision observables in eq. (6.5) (without COHERENT) to constrain the u,d coefficients, finding u = 0.003(10) and d = 0.019 (21). This constraint is also percent level, indicating that COHERENT has an important impact on the global fit. Indeed, the combination of the COHERENT results with other precision observables leads to u = −0.0037(54), d = 0.0054 (93), which represents a factor of two improvement. These results are represented in the right panel of figure 6.
All in all, the results of this section demonstrate that COHERENT has become an indispensable ingredient in the family of electroweak precision observables constraining the SMEFT Wilson coefficients.

Conclusions and discussion
In this paper, we have laid down a new theoretical framework based on effective field theories to describe non-standard effects in coherent neutrino scattering on nuclei. The framework is very versatile, allowing us to handle simultaneous new physics contributions to neutrino production and detection, non-linear effects of non-standard Wilson coefficients, different input schemes for the SM parameters, as well as an arbitrary flavor structure of neutrino-matter CC and NC interactions. It can also be readily applied to EFTs at different energy scales (nucleon-or quark-level, below or above the electroweak scale. . . ). This generalizes the NSI language that is followed by the greater part of the literature, and facilitates connection to specific BSM models with new particles heavier than the characteristic experimental scale.
An important element of our analysis is that we have included new physics effects coming from both the neutrino production and detection processes. There exists a nontrivial interplay between these two pieces, which cannot be reduced to a simple factorization of the neutrino flux and the detection cross section. We remark that correlated effects in production and detection are rather generic in new physics models. In particular, in the SMEFT framework several dimension-6 operators contribute to both, due to the SU(2) L gauge symmetry relating the CC and NC interactions. The consideration of these correlated effects opens the door to using the measurements at COHERENT to probe new physics models that relate CC and NC interactions.
We have introduced three generalized weak chargesQ f , which, in a certain sense, can be associated to the production and scattering of ν e , ν µ andν µ on the nuclear target, cf. eq. (3.17). 16 These can be extracted from the COHERENT data, given the recoil energy and time distribution of the nuclear recoil events. They contain full information about the contributions of new physics affecting the effective contact interactions between neutrinos and matter. The framework can be simplified by folding in further assumptions. If new physics affects only detection, or if one restricts to a linear order in the EFT Wilson coefficients, there remain two independent generalized weak charges (electronic

JHEP05(2023)074
and muonic), similarly as in the prior analyses within the NSI framework. In the SM limit there is a single nuclear weak charge that needs to be extracted from experiment.
We have applied this framework to obtain novel constraints on new physics based on the results of the COHERENT experiment. We have examined the full dataset available at the moment, which includes distributions of energy and time of recoils measured in cesium iodine and argon nuclear targets. The central result of our numerical analysis is given in eqs. (5.1) and (5.5), where we give the confidence intervals for the three generalized weak charges, including the correlations. These encode the state-of-the-art description of the CO-HERENT constraints on EFTs. Let us stress the approximate Gaussianity of the results when the squared charges are used. This means that our entire analysis of LAr and CsI data, which contains 664 bins, with thousands of associated backgrounds, and 13 nuisance parameters, can be expressed in terms of three central values, three diagonal errors and one correlation matrix (per target). Any BSM practitioner can then trivially apply these results to their particular NSI setup, EFT approach, or New Physics model. We encourage the community to provide the result of their analyses of COHERENT data in this convenient form.
We have recast these limits into bounds on EFT Wilson coefficients, both in the WEFT and in the SMEFT. We find that the COHERENT data so far provide percent level constraints on two particular combinations of Wilson coefficients in the EFT parameter space. We demonstrate that these constraints are a valuable ingredient in the grander scheme of electroweak precision observables. The impact of the coherent neutrino scattering information is most relevant when performing a global fit, both in the constrained flavor-blind (U(3) 5 -symmetric) setup, and in the completely generic scenario. From the point of view of the SMEFT fits, COHERENT is the most relevant neutrino-detection experiment, clearly superior in comparison to previous neutrino scattering experiments at higher energies.
All in all, in this work we have carried out a complete study of New Physics effects at the COHERENT experiment within an EFT approach. This includes for the first time simultaneous effects in neutrino production and detection, and its addition to the global SMEFT fit of electroweak precision data. Our work enables the study of the impact of COHERENT to new BSM setups and opens exciting possibilities for future developments. For instance, our work can be extended to other scenarios with light exotic particles, or to the study of additional experiments in the broad and flourishing CEνNS landscape.

A Details of the calculation of the event rate
In this section we elaborate on some details of our calculation of the CEνNS event rate. First, we show the expressions for the f f (T ) functions, which contain part of the kinematic dependence of the predicted prompt and delayed number of events in eq. (3.16): where, as was mentioned in the main text, E min ν (T ) denotes the minimum energy of the (anti)neutrino required to produce CEνNS with a recoil energy T , which is given by Next, we expand the information presented in eq. (3.19), which we reinstate here: Here, the fluxes are defined in the usual form: As for the cross sections, they are also defined in the usual form but using the generalized chargesQ f , that is

B Numerical analysis: further details
In this section we describe in detail the input used in our numerical analysis, which is chosen in every case following closely the corresponding COHERENT prescription.
Before discussing the details that are specific to each measurement, let us show the expression that we use for the form factor, F (q 2 ), since that is a common input to all cases. We use the Helm parametrization [47], which gives the following expression for the neutron and proton form factors: Here j 1 (x) is the order-1 spherical Bessel function of the first kind, s = 0.9 fm is the nuclear skin thickness [65] and R 0,p/n is a function of s and the proton/neutron root-mean-square (rms) radius R p/n given by The proton and neutron rms radii for the studied nuclei are taken from refs. [38,66,67] R p (Cs) = 4.821 (5)  Following the COHERENT prescription, the uncertainty associated to this description of the form factor is included in our analysis through a nuisance parameter, as described below.
For the CsI analysis, we will take the average of R p/n (Cs) and R p/n (I) and of the nuclei masses. Moreover, as discussed in the main text, in both the LAr and the CsI cases we approximate neutron and proton form factors to be equal by taking the average of R p and R n . This simplifies significantly the presentation of intermediate results and it is not expected to have any impact in the final results for the WEFT Wilson coefficients, taking into account current COHERENT uncertainties.
Finally, we do not consider the contribution from neutrino electron scattering, which acts as an additional background source. This contribution was separated from the signal in the LAr analysis, but not in the CsI analysis. In both cases, its effect on the SM event rate is negligible, while for heavy NP scenarios its impact is also irrelevant [32,34].

B.1 LAr measurement
The LAr dataset consists of a 3D distribution in recoil energy, time and the fraction of integrated amplitude within the first 90 ns after trigger (F 90 ) [35]. The latter plays no direct role in our analysis, so we integrate over F 90 and work with the resulting 2D recoil and time distribution. Our analysis covers the range 0 < T rec ee (keV) < 40 and 0 < t(µs) < 5, using 4x10 bins of equal width. These are the bins with a significant amount of CEνNS events.

JHEP05(2023)074
The expected number of events per bin is given by where pBRN, dBRN, SS correspond to prompt BRN, delayed BRN and SS backgrounds (NIN contribution is neglected in this analysis). Their predicted values, N bkg,a ij , are readily provided in the LAr measurement data release [52], with efficiencies already applied to them.
Thus, the nuisance parameters in this analysis are x = {α, β pBRN , β dBRN , β SS , z }, with 1 ≤ z ≤ 5, with uncertainties equal to {13%, 32%, 100%, 0.8%} and 100% for all five z parameters. We see that the generic h ij functions introduced in eq. (4.5) are in this case We have nuisance parameters associated to the systematic uncertainties of the overall normalization of backgrounds (β a ) and signal (α). The latter includes errors associated to detector efficiency, energy calibration, F 90 calibration, quenching factor, nuclear form factor and neutrino flux. In addition, this analysis includes systematic uncertainties affecting the shape of the distributions ( z ). In particular, we have two systematic errors affecting the signal distribution, coming from the energy dependence of the F 90 distribution and the trigger time mean, and three systematic errors affecting the prompt BRN distribution, with the energy distribution, the trigger time mean and the trigger time width as their sources. These bin-dependent systematics are included through the following functions: where N 1σ z,ij is the predicted number of events with a 1-σ shift due exclusively to the zth systematic error and N CV ij is the predicted central value. Note that these quantities include the total number of events and not only the (signal/pBRN) events affected by the z-th systematic error. We take the five 1-σ distributions (1 ≤ z ≤ 5) and the three background CV's (pBRN, dBRN,SS) from the COHERENT data release [52].
The number of expected CEνNS events is obtained using eq. (4.4), which we repeat here The timing information (i.e., the g j factors) is extracted from the neutrino flux characterization presented in the COHERENT data release [52,68]. The energy distributions, N i , are calculated using eq. (4.3), which involves an efficiency function, energy resolution and quenching factor that we describe below. The quenching factor is parametrized through a polynomial expression, given by where a QF = 0.246 keV −1 and b QF = 0.00078 keV −2 .

JHEP05(2023)074
The detector resolution function is where the T ee -dependent width is given by σ ee = 0.58 keV T ee /keV. For the lower limit of the T integration in eq. (4.3) we do not use zero but T min = 79 eV (average energy to produce a scintillation photon in Ar [69]), but this has a negligible impact in our results. The efficiency function, (T rec ee ), used in the calculation of the CEνNS events is given by COHERENT as a T rec ee bin dependent quantity [52].

B.2 CsI measurement
Our CsI analysis uses the 2D distribution in recoil energy and time covering the ranges 8 < PE < 60 and 0 < t(µs) < 6 with 1 PE and 0.5 µs of width respectively, which yields a total of 52x12 bins. This is the same bin set as the one used in the original COHERENT analysis [36].
The formula for the expected number of events per bin is given by The α parameter encodes the systematic uncertainties associated to the signal (due to the QF, neutrino flux and form factor), whereas the β ones encode the uncertainty associated to the normalization of the backgrounds.
As in the LAr case, the expected number of CEνNS events, N signal ij , is calculated using eq. (4.4) (repeated above), and the prompt and delayed energy distributions, N a i , are obtained using eq. (4.3). The latter involves an efficiency function, energy resolution and QF that we describe below. Finally, the time information, i.e., the g j factors in eq. (4.4), is described afterwards.
The COHERENT prescription for the QF in this analysis is the following [36] The measured events are organized in PE bins, so we apply the light yield, given in this case by LY=13.35 PE/keV. The energy resolution function is given by where a and b encode the T ee dependence: a = 0.0749 keV/T ee and b = 9.56 keV −1 × T ee .

JHEP05(2023)074
The energy-dependent efficiency applied in this measurement is where a = 1.320 ± 0.023, b = (0.28598 ± 0.00061) PE −1 , c = (10.9 ± 1.0) PE, d = −0.333 ± 0.023. We have checked that these uncertainties have a negligible effect in our fit and thus we have neglected them. As for the timing information (g j factors), we extract them from the information about the flux for every neutrino flavor, which is binned in time and recoil energy and provided in the data release [52]. Integrating over the recoil energy we obtain the prompt and delayed distributions. Then we take into account the timing efficiency of the detector, which is available for this measurement. Namely where a = 0.52 µs and b = 0.0494 µs −1 . Once again, the impact of the uncertainties of these parameters in our fits is negligible. We apply this timing efficiency to the projected time distributions, and then we normalize each of them to recast them as probability distribution functions. The ν µ flux gives the prompt distribution, and the ν e (orν µ ) flux can be used to obtain the delayed distribution.
Concerning the backgrounds, the BRN and NIN distributions are provided in the COHERENT data release [52]. They can be normalized to produce 1D distributions in the recoil energy and timing directions, denoted by g a (P E) and f a (t) respectively (a =BRN, NIN). The full 2D distributions are obtained just by taking N a tot g a (PE) f a (t) t (t), where N BRN tot = 18.4 and N NIN tot = 5.6 are the total number of predicted BRN and NIN events. The SS background can be estimated from the anti-coincidence data (AC), which is also given in a 2D distribution in recoil energy and time. The projection onto the PE axis provides directly the recoil-energy distribution, while for the description of the time evolution of this background, the collaboration advises the use of an exponential model. That exponential is fitted to the projection of the AC data on the time axis and then normalized, yielding an expression f SS (t) ∝ e −a SS t , with a SS = −0.0494 µs −1 . This procedure (instead of working directly with the 2D AC distribution) avoids possible biases in the fit due to limited statistics in the sample [3]. Since this distribution is inferred directly from the data, it is not necessary to apply efficiencies here.

C Details of the SMEFT analysis
In this section we provide additional details regarding our SMEFT fit combining the CO-HERENT results with other low-and high-energy precision measurements. We list the observables used in the fit, and discuss the updates with respect to the analogous fit in ref. [40]. We also present the global fit to the Wilson coefficients without assuming the U(3) 5 flavor symmetry.

C.1 Operators
If the BSM particles are not only heavier than the characteristic scale of COHERENT but also heavier than the electroweak scale, then above that scale their dynamics can be described by the effective theory called SMEFT [16,17]. The latter contains all the SM fields (including the heaviest ones, such as the Higgs and the top quark) and the SU(3) × SU(2) × U(1) gauge symmetry is (linearly) realized. The leading effects of leptonnumber-conserving BSM particles are encoded in the dimension-6 operators in the SMEFT Lagrangian. We parametrize the dimension-6 operators using the so-called Higgs basis [70] (see [63] for a recent review) as it is very convenient for the sake of electroweak precision fits. The operators involved in our fit are a subset of 4-lepton and 2-lepton-2-quark operators, as well as vertex corrections to electroweak gauge bosons interactions with fermions: The operators in the three categories entering our fit are defined as  17 The four-fermion pieces, L 2l2q and L 4l , are manifestly invariant under the SM gauge group, and are actually the same as the analogous terms in the Warsaw basis. Concerning L vertex , although it is manifestly invariant only under SU(3) C × U(1) em , it is obtained from an SU(3) C × SU(2) L × U(1) Y invariant Lagrangian after electroweak symmetry breaking. In this case the imprints of the full SM gauge symmetry are the correlations between the W and Z couplings to fermions. The map between the vertex corrections and the Wilson coefficients of the manifestly SM gauge-invariant Warsaw basis is In the U(3) 5 symmetric limit, assumed in the main body of this paper, our Wilson coefficients reduce to

C.2 Experimental input
Here we list the observables that are included in our update of the global SMEFT fit carried out in ref. [40], with special emphasis on the changes with respect to that work: • e + e − collisions at energies above [71,72], at [37,73,74], and below [75][76][77] the Z-pole.
Concerning the Z-pole results we follow ref. [78], where we took into account recent theoretical calculations [79,80] that lead to minor modifications of the Z width, the hadronic cross section, and the forward-backward asymmetry of b quarks.
• For W boson data we follow ref. [78] as well: we include the mass and total width [37], leptonic branching ratios from LEP [72], Tevatron [81] and LHC [82][83][84], and the ratio R W c ≡ Γ(W → cs)/Γ(W → ud, cs) [37], which has a limited precision, but helps to break a flat direction. These data include three significant updates with respect to the fit in ref. [40]: (i) the recent LHC results concerning the leptonic branching ratios [82][83][84], which play an important role in removing certain LEP-2 tensions and improving the constraints on the W ν vertices; (ii) the W boson mass, which is updated using the current PDG combination m W = 80.377 (12) GeV [37,72,[85][86][87][88]; 18 (iii) we no longer use ref. [90] to extract the W t L b L vertex (and hence [δg Zu L ] 33 ), as this measurement is also sensitive to other operators in the top sector.
The LHC bounds, which were obtained in ref. [78], represent a novel addition with respect to the fit in ref. [40] and they tighten the LEP constraints on the Z boson couplings to first generation quarks.
• Muon-neutrino scattering on nuclei in the CHARM [94], CDHS [95], and CCFR [96] experiments. As in ref. [40], we use the PDG combination of these data, which includes as well additional experimental input (with very limited precision) from elastic neutrino-proton scattering and neutrino-induced coherent neutral-pion production from nuclei.
• Parity violation in atoms and in scattering: (i) measurements of atomic parity violation in cesium [97] and thallium [98,99] atoms; (ii) the weak charge of the proton measured in scattering of low-energy polarized electrons by QWEAK [100]; and (iii) deep-inelastic scattering of polarized electrons on deuterium by the PVDIS experiment [101]. As in ref. [40], we use the PDG combination of these data, supplemented by the SAMPLE measurement of the scattering of polarized electrons on deuterons in the quasi-elastic kinematic regime [102]. The only difference with ref. [40] is that we use the 2022 PDG combination (table 10.9 of ref. [37]), which includes the updated measurement of the proton weak charge by the QWEAK experiment [100].
• Parity-violating scattering of electrons at low energies in the SLAC E158 experiment [112]. As in ref. [40], we use the PDG combination for the low-energy parityviolating electron self-coupling [37].
• Leptonic decays of taus and muons [37]. Unlike [40], we no longer use the ratio of effective Fermi constants G τ µ /G F from the τ -Lepton Decay Parameters review in ref. [115], since its large correlation with the poorly known [c le ] µτ τ µ coefficient was not provided. Instead, we directly use the measured τ → eνν branching fraction and the associated Michel parameter [37]. For completeness we also include the measurement of the Michel parameter in muon decay [37]. These modifications allow us to target additional 4-fermion semileptonic operators, not constrained by the fit in ref. [40], namely the [c le ] µτ τ µ and [c le ] eµµe coefficients.
• We include a set of hadronic tau decay observables described in ref. [105]. These constraints, which were obtained subsequently to the global fit of ref. [40], give access to contact interactions between first generation quarks and third generation leptons. They also improve our knowledge of the vertex corrections.
We refer to refs. [40,78,105] and to the original experimental papers for the central values and errors used in our fit. We calculate these observables at tree level in the SMEFT in terms of the Wilson coefficients in eqs. (C.2), (C.3) and (C.4), and expand them to linear order. In other words, we keep track of the O(1/Λ 2 ) effects in the SMEFT power counting, while ignoring O(1/Λ 4 ) effects due to dimension-8 operators and squares of dimension-6 operators. We also ignore all loop effects (except for QCD running of certain WEFT operators below the electroweak scale). 19 There is a nonzero correlation between up-down effective couplings (the focus of this work) and upstrange couplings (which we marginalized over). This must be taken into account when going to specific scenarios, such as the one-a-time results in table 3 or the flavor-blind fit in eq. (6.5). The full likelihood is available in ref. [105].

JHEP05(2023)074
In addition to the inputs above, we consider the COHERENT constraints, which represent the main topic of this paper. Below we incorporate them for first time in this kind of global SMEFT fit. For that, the constraints on the WEFT Wilson coefficients listed in eq. (5.10) are translated into constraints on the SMEFT Wilson coefficients using the map uu αα = δg Zu L + δg Zu R + 1 − lq + c ld ] αα11 . (C.8)

C.3 Fit results
As discussed at length in ref. [40],  Using these variables, the global likelihood depends on the Wilson coefficients on the righthand sides of eqs. (C.9) only via theĉ and d P (2 GeV) combinations. Let us stress that the Wilson coefficients in the r.h.s. of d P (2 GeV) in eq. (C.9) are defined at µ = m Z .

JHEP05(2023)074
All in all, we simultaneously fit 65 independent (combinations of) SMEFT Wilson coefficients, including all correlations. We find the following 1σ intervals     18(20) We restrain ourselves from copying here the 65 × 65 correlation matrix. This, as well as the full Gaussian likelihood function in the Higgs or Warsaw basis, is available in the numerical form on request. We highlighted in red the Wilson coefficients whose bounds have improved significantly thanks to including the COHERENT results. The improvement is illustrated in figure 7. The most spectacular effect is observed for the combination [ĉ eq ] ee11 = [c eq ] ee11 + [c (1) lq ] ee11 . In the fit of ref. [40] it was constrained only by the CHARM measurement of electron neutrino scattering on nuclei [93], which is however very imprecise. The COHERENT results analyzed in this paper reduce the error bars on [ĉ eq ] ee11 by a factor of 5.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.