Reactor neutrino oscillations as constraints on Effective Field Theory

We study constraints on the Standard Model Effective Field Theory (SMEFT) from neutrino oscillations in short-baseline reactor experiments. We calculate the survival probability of reactors antineutrinos at the leading order in the SMEFT expansion, that is including linear effects of dimension-6 operators. It is shown that, at this order, reactor experiments alone cannot probe charged-current contact interactions between leptons and quarks that are of the (pseudo)vector (V$\pm$A) or pseudo-scalar type. We also note that flavor-diagonal (pseudo)vector coefficients do not have observable effects in oscillation experiments. In this we reach novel or different conclusions than prior analyses of non-standard neutrino interactions. On the other hand, reactor experiments offer a unique opportunity to probe tensor and scalar SMEFT operators that are off-diagonal in the lepton-flavor space. We derive constraints on the corresponding SMEFT parameters using the most recent data from the Daya Bay and RENO experiments.


Introduction
Since the discovery of neutrino oscillations, the field of precision neutrino physics has experienced a formidable rate of progress. Assuming the standard 3-flavor picture, the mass squared differences between the neutrino eigenstates and all three angles in the mixing matrix have been determined with a good precision, see ref. [1] for a recent update. The standard parameters are now overconstrained by multiple independent measurements, with overall a good consistency. In a way, the situation is similar to that in electroweak precision physics in the 1990s when, given the wealth of precise and theoretically clean information from LEP-1, the initial focus on measuring the parameters of the Standard Model (SM) could be extended to constraining hypothetical phenomena (technicolor, supersymmetry, JHEP05(2019)173 etc.). By the same token, neutrino experiments now have a potential to systematically explore new physics beyond the neutrino masses and mixings.
One such area of exploration are the so-called non-standard interactions (NSI). Oscillation experiments are sensitive not only to neutrino masses and mixing, but also to how neutrinos interact with matter. The SM makes precise predictions about these interactions, which however can be perturbed by physics beyond the SM (BSM). In particular, new effective 4-fermion interactions between leptons and quarks may give observable effects in neutrino production, propagation, and detection, and thus they can be constrained by experiment. These studies have a long history in the literature, see e.g. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] and [20] for a recent review.
Previous works were contended with an ad-hoc parametrization of NSI, using the effective couplings s,d,m to describe the non-standard effects in production, detection, and propagation of neutrino states (see appendix B for more details). There has been little emphasis on connecting these couplings to parameters of concrete BSM models, or to Wilson coefficients of a well-defined and systematic framework of effective field theory (EFT). Consequently, full attention has not been paid to such issues as power counting of NSI effects, extraction of the mixing angles in the presence of general new physics, or comparison between the sensitivity of oscillation and other experiments. We argue here that there are distinct advantages in embedding NSI in a solid EFT. First, consistent EFTs come with an expansion parameter, and the Lagrangian, amplitudes, and observables can be systematically constructed order by order in that expansion. This allows one to compare different NSI effects in neutrino oscillations, and unambiguously identify the leading order contributions. Moreover, EFTs may predict correlations between the magnitude of effects in neutrino oscillation and in other precision experiments, such as nuclear beta transitions, meson decays, Drell-Yan production at the LHC, etc. In this picture, oscillation experiments become an ingredient in the broad program of precision measurements. Moreover, sensitivities of the different precision probes can be meaningfully compared.
In this paper we propose a systematic EFT approach to neutrino oscillations. We focus on short-baseline reactor neutrino experiments, however the formalism can be readily applied to experiments with longer baselines and for different neutrino production and detection processes. Our point of departure is the so-called SMEFT, where higher-dimensional interactions invariant under the local SU(3) C × SU(2) L × U(1) Y symmetry are added to the SM. They are organized in an expansion in 1/Λ, where Λ can be interpreted as the BSM scale suppressing the higher-dimensional operators. Given the SMEFT Lagrangian, we derive the effective charged-current interactions between neutrinos, charged leptons, and nucleons in the low-energy EFT relevant for reactor experiments. This allows us to calculate the survival probability of an electron antineutrino in short-baseline experiments in the presence of general dimension-6 (order Λ −2 ) SMEFT interactions. We point out that these experiments offer a unique opportunity to probe, at the linear level, certain SMEFT operators that are off-diagonal in the lepton-flavor space. We identify the linear combinations of Wilson coefficients that are constrained by reactor experiments at the leading order in the SMEFT expansion. We then proceed to obtain numerical constraints on these combinations using the most recent data from the Daya Bay [21] and RENO [22] experiments.

JHEP05(2019)173
Our systematic approach puts into perspective some of the conclusions reached in the prior NSI literature. We will argue that, at the leading order in the SMEFT expansion, NSI interactions diagonal in the lepton-flavor space are currently not probed by oscillation experiments. More precisely, any modifications of diagonal V±A interactions are fully absorbed in the phenomenological extraction of SM parameters (the CKM element V ud and the neutron axial charge g A ), while for the scalar and tensor ones stringent modelindependent constraints from nuclear beta decays exclude observable signals in reactor experiments, given the current precision of the latter. As for NSI off-diagonal in the lepton-flavor indices, those involving only left-handed leptons and quarks (V-A type) cannot be constrained by the reactor experiments alone, as they merely renormalize the a-priori unknown mixing angle θ 13 . Off-diagonal NSI with right-handed quarks (V+A) are in principle observable in the reactor oscillation experiments, however they do not arise at O(Λ −2 ) from dimension-6 SMEFT operators. On the other hand, reactor experiments show an interesting sensitivity to off-diagonal tensor and scalar NSI, which were actually neglected in most prior studies. This paper has the following structure. In section 2 we review the formalism of the SMEFT and the resulting EFT below the weak scale. In section 3 we derive the dependence on the SMEFT parameters of the anti-neutrino survival probability in reactor experiments. The constraints on these parameters from the Daya Bay and RENO observations are presented in section 4, and compared in section 5 to the constraints from other precision experiments. We summarize our findings in section 6 and comment on the significance of our results on the NSI program in oscillation experiments.

EFT ladder
As mentioned above, the oscillation pattern of neutrinos depends not only on their mass differences, but also on their interactions with other particles. We are interested in the situation where these interactions deviate from the SM predictions. Our goal is to derive new constraints on fundamental theories with heavy BSM particles, without referring to a specific model. For that reason we will use the language of EFT. In this section we review the crucial elements of EFTs relevant for our analysis.

SMEFT
If new particles beyond the SM are much heavier than the Z boson and the electroweak symmetry breaking is linearly realized, then the relevant effective theory above the weak scale is the so-called SMEFT [23,24]. It has the same local symmetry and particle content as the SM, which in particular entails the absence of right-handed neutrinos. But the SMEFT differs from the SM by the presence of higher-dimensional (non-renormalizable) interactions in the Lagrangian, which provide an effective description of physical effects of heavy BSM particles. They are organized in a systematic expansion in operator dimensions, with each consecutive terms suppressed by a higher power of the new physics scale Λ. Dimension-5 operators are essential as they give rise to Majorana masses of the SM neutrinos. The formulas presented in this work assume the normal ordering of neutrino masses, JHEP05(2019)173 but the changes in the case of inverse ordering are trivial, as we will explicitly discuss. Due to the smallness of the neutrino masses, dimension-5 operators have negligible effects on production and detection amplitudes of relativistic neutrinos. On the other hand, these can be significantly affected by dimension-6 operators suppressed by Λ −2 . In particular, some dimension-6 operators lead to deviations of the couplings of the SM quarks and leptons to the W boson from the SM prediction; others introduce new contact interactions between quarks and leptons. In our study we will ignore the effects of operators with dimensions higher than six, which are suppressed by more than two powers of Λ. Consequently, we will only trace new physics corrections linear (order Λ −2 ) in Wilson coefficients of dimension-6 operators, and ignore the quadratic effects that are O(Λ −4 ).

WEFT
The SMEFT is a convenient tool when it comes to studying high-energy physics above the weak scale. However, neutrino oscillation experiments are performed at energies well below the weak scale. At the scale µ m W , the W and Z bosons, as well as the Higgs boson and the top quark, can be integrated out from the SMEFT, leading to another effective theory that we refer to as the weak EFT (WEFT). 1 It has a smaller particle content and different interactions than the SMEFT. Below we focus on the charged-current 4-fermion interactions between the up and down quarks and the 3 generations of charged leptons and neutrinos. At the leading order in the WEFT we can parametrize them as where v is the VEV of the Higgs doublet, V ud is a CKM matrix element, α = e, µ, τ is a charged lepton field, σ µν = i[γ µ , γ ν ]/2, and P L,R are the usual chirality projectors (1 ∓ γ 5 )/2. 2 Above, the fields u, d, α are written in the basis where their mass terms are diagonal. The flavor neutrino states ν α are connected to the mass eigenstates by ν α = U αJ ν J , where α = e, µ, τ , J = 1, 2, 3, and U is the unitary PMNS matrix parametrized by three mixing angles (θ 12 , θ 13 , θ 23 ) and one CP-violating phase δ CP : 1 In this work we consider the WEFT as the low-energy theory of the SMEFT. However, the WEFT is a consistent EFT in its own right, which can be valid even if it is not completed by the SMEFT at higher energies. See appendix A for the discussion of such a set-up. 2 The hat overˆ T indicates that the normalization differs by a factor of 4 from T used e.g. in [25].
The present normalization is more natural in the sense that typical new physics models generating tensor interactions will give comparable contribution toˆ T and S,P , see e.g. [26].

JHEP05(2019)173
and s ij ≡ sin θ ij , c ij ≡ cos θ ij . As mentioned earlier, the neutrinos have only the left-handed components, while right-handed neutrinos are absent in this effective theory. The leading NSI corrections to the standard neutrino interactions are summarized by the parameters X in eq. (2.1), which are 3 × 3 matrices in the lepton flavor space. In the neutrino literature these are customarily referred to as the NSI. Apart from the SM-like V-A interactions (1 + L ), right-handed ( R ), scalar ( S ), pseudoscalar ( P ), and tensor (ˆ T ) interactions between leptons and quarks are allowed at the same order of the WEFT expansion. The matching between X and the Wilson coefficients of dimension-6 operators in the Warsaw basis [24] of the SMEFT at the renormalization scale µ ∼ m W is given by [27][28][29][30] [ where SMEFT operators are defined in the flavor basis where the up-quark Yukawa matrices are diagonal. There are three important conclusions from this matching exercise. Firstly, all the X parameters in eq. (2.1) arise at O(Λ −2 ) in the SMEFT, thus a priori they are equally important. Secondly, the right-handed interactions are proportional to the unit matrix in the lepton flavor space, up to corrections from dimension-8 and higher SMEFT operators [27]. Indeed, at the dimension-6 level R can originate only from the operator O Hud = iH T D µ H(ū R γ µ d R ) and its conjugate, which induce the W boson coupling to right-handed quarks. Integrating out the W exchange between the quarks and leptons generates R in eq. (2.3). Since the SM W couplings to leptons are diagonal and flavor universal, so is R at O(Λ −2 ). Off-diagonal and flavor non-universal contributions to R can appear only at O(Λ −4 ), either from the W exchange (if the W couples to right-handed quarks and non-universally to leptons at order Λ −2 ), or from dimension-8 contact operators such as e.g (L α Hγ µ L β H)(ū R γ µ d R ), where L α = (ν L , L ) α are the lepton doublets, and H is the Higgs doublet. On the other hand, L,S,P,T do contain off-diagonal and nonuniversal terms already at the dimension-6 level, in general. Finally, L is approximately a Hermitian matrix in the lepton flavor space, up to corrections suppressed by off-diagonal CKM matrix elements. This directly follows from the hermiticity properties of the SMEFT Wilson coefficients: [c Hq ] kj , and [c

Lee-Yang
At the energy scale characteristic for reactor neutrino experiments the relevant degrees of freedom are not quarks, but rather their bound states such as nucleons and nuclei. Therefore, it is advantageous to descend one more step in the EFT ladder, into an effective theory JHEP05(2019)173 of protons and neutrons interacting with charged leptons and neutrinos. Matching this EFT to the WEFT Lagrangian in eq. (2.1) we obtain the so-called Lee-Yang Lagrangian [31]: where p, n are relativistic proton and neutron fields. The couplings g V,A,S,P,T are vector, axial, scalar, pseudoscalar, and tensor charges of the nucleon, which can be calculated on the lattice or from symmetry considerations. For the vector coupling, one can prove that g V = 1 up to quadratic corrections in isospin-symmetry breaking [32]. For the remaining charges we use the numerical values collected in table 1 of [33] (which are taken from refs. [34,35]), except for g A , which is taken from the fit in eq. (84) of [33]: g A = 1.2728 ± 0.0017, g S = 1.02 ± 0.11, g P = 349 ± 9, g T = 0.987 ± 0.055, (2.5) at µ = 2 GeV in the M S scheme, a choice of scale and scheme that will apply as well to the X bounds obtained in this work. For our purpose, the charges are known with sufficient precision, and we will use their central values ignoring the errors. The processes relevant for neutrino production and detection in reactor experiments are (inverse) beta decays. While in these reactions neutrinos and electrons are typically relativistic, the exchanged momenta are much smaller than the masses of nucleons. Therefore, one can describe the latter using non-relativistic fields ψ in an effective theory expanded in powers of the spatial derivatives ∇ k ψ. At the leading (zero-derivative) order, the Lee-Yang Lagrangian in eq. (2.4) reduces to where ψ p and ψ n are non-relativistic fields annihilating protons and neutrons, respectively, and Σ k = 0 σ k σ k 0 with σ k being the Pauli matrices. Note that P does not appear in eq. (2.6), hence the pseudoscalar interactions do not affect beta transitions at the leading order. Moreover there are only two independent hadronic structures at this order: ψ p ψ n andψ p σ k ψ n , which mediate the Fermi and Gamow-Teller nuclear transitions, respectively. Continuing the non-relativistic expansion, at the next order one would obtain the interactions with one derivative acting on ψ, which lead to the so-called first-forbidden beta transitions.
It is worth noting that the same effective interactions parametrized by X that can be probed in neutrino oscillation experiments also affect the phenomenological extraction of JHEP05(2019)173 V ud and g A from nuclear and neutron decays, respectively [33]. Since the latter quantities are needed to calculate the predicted number of produced and detected neutrino events in oscillation experiments (see section 4.1 for more details), these effects have to be taken into account consistently at the chosen order in the EFT expansion. Such a consistent analysis has never been done in previous NSI literature, to the best of our knowledge. In particular, it can be shown that neutrino oscillation data does not depend on the flavor-diagonal vector EFT couplings [ L ] ee and Re [ R ] ee at any order. This is so because their direct effect is completely cancelled by the indirect effect entering through the phenomenological determination of V ud and g A . This can be shown at the Lagrangian level, since these nonstandard contributions only appear in the following two combinations [33] Precision experiments that provide the numerical values of V ud and g A in fact measure the above combinations of the SM and NSI parameters, when interpreted in the EFT context. 3 [20]). Let us note that the [ L,R ] ee coefficients can be probed in precision beta decay measurements, through a (firstrow) CKM unitarity test and through the comparison of lattice and experimental values of g A . The resulting bounds are below the permil and percent level, respectively [33]. In the following we will use the Lagrangian in eq. (2.6) to calculate amplitudes of beta decay processes relevant for reactor neutrino oscillations. We will treat X as small parameters of order Λ −2 , as derived from the matching to the SMEFT, and we will ignore any contributions to observables that are O(Λ −4 ) or smaller.

Oscillations in EFT
In this section we review the theory of neutrino oscillations in the presence of NSI. We focus on providing a systematic EFT description of new physics effects in neutrino production and detection in short-baseline reactor experiments. We neglect matter effects in neutrino propagation, which would be relevant for long-baseline experiments.
Consider an antineutrino produced with energy E ν in the process X P → − αν Y P and detected in the processνX D → + α Y D , where α is a charged lepton: electron, muon, or tau. Given a neutrino produced in association with − α , the survival probability is defined as the probability of it being detected at the distance L from the source in association with + α of the same lepton flavor. Quite generally, the survival probability is given by the formula [36]:

JHEP05(2019)173
where the indices J, K, . . . label neutrino mass eigenstates ν J , ∆m 2 JK ≡ m 2 ν J − m 2 ν K , and A P αJ and A D Jα denote the amplitudes for the production and detection of ν J : In neutrino oscillation experiments, polarization of particles involved in production and detection is not measured, therefore summation over spins (and any other internal indices) is implicit in each bracket in eq. (3.1). Likewise, there is an integration over all kinematic variables (except the neutrino energy), as indicated by in eq. (3.1). We now derive a general expression for the coefficients C α JK as a function of the Wilson coefficients X in the WEFT Lagrangian eq. (2.1). The amplitudes in eq. (3.2) can be decomposed as (3.3) Here M P X and M D X are independent of the mass index of the emitted/absorbed antineutrino, up to totally negligible corrections due to the neutrino masses. Then, keeping only the linear effects in X , we can approximate (3.5) The first line in eq. (3.4) encapsulates the standard oscillations in the absence of BSM effects other than the neutrino masses. The second and third lines in eq. (3.4) describe corrections to the survival probability due to NSI affecting, respectively, the neutrino production and detection processes. The coefficients p X and d X depend on the processes in which neutrinos are produced and detected, and in general they may be functions of the neutrino energy. Note that the diagonal elements X do not enter eq. (3.4); in fact they cancel out between the numerator and denominator of eq. (3.1). Therefore, only the off-diagonal (in the charged-lepton flavor basis) Wilson coefficients of the effective Lagrangian eq. (2.1) affect the survival probability at the leading order. Recall that if the WEFT Lagrangian is derived from the underlying SMEFT (that is, if new physics is heavier than m W and respects the full SM local symmetry) then R is a diagonal matrix, leading to the conclusion that the charged currents involving right-handed quarks do not affect neutrino oscillations at O(Λ −2 ). Specializing to reactor experiments such as Daya Bay and RENO, neutrinos are detected via inverse beta decay on water (practically, proton) targets, with a positron and a JHEP05(2019)173 neutron in the final state:νp → ne + . Calculating the amplitude for this process starting from the non-relativistic effective Lagrangian in eq. (2.6) we find the following detection coefficients MeV and m e ≈ 0.511 MeV is the positron mass. The same result is obtained starting from the relativistic eq. (2.4) in the limit where the proton recoil is neglected. That calculation reveals that the contribution proportional to P is suppressed by the small factor g P m e /m p ∼ 0.1 in spite of the large value of the pseudoscalar charge g P . In the following we neglect these subleading pseudoscalar contributions. Note that d S and d T depend on the neutrino energy, which will be an important handle for constraining the scalar and tensor Wilson coefficients in reactor experiments. The factor in the amplitude proportional to me Eν −∆ ≈ me Ee goes under the name of the Fierz interference term [37], and is due to the lepton-chirality flip in the corresponding Lagrangian terms in eq. (2.6).
While non-standard effects on the detection side are calculable to a good accuracy, the production side is far more involved. There are hundreds of different beta decay processes contributing to the antineutrino flux in the reactor [38,39], and the NSI effects on their amplitudes may be subject to relatively large uncertainties. To tackle that problem, we have to resort to certain crude approximations. First, we assume that all beta decays contributing to the reactor antineutrino flux above the detection threshold E ν = 1.8 MeV are of the Gamow-Teller type. With that assumption, the production coefficients are given by .
As before, the pseudoscalar interactions can be neglected at the leading order. In addition the scalar ones do not contribute to Gamow-Teller transitions. The form factor in the tensor coefficient is given by . (3.8) The sum in the first line goes over all β decays resulting from nuclear fission processes in the reactor, with appropriate weight factors w i determined by the fission yield. Furthermore, ∆ i are the mass differences of the initial and final state nuclei participating in the beta processes. As shown in the second line of eq. (3.8), rather than using a detailed reactor model with all distinct processes explicitly included, to calculate f T (E ν ) we replace the sums by integrals over endpoint energies. We use a gaussian distribution for W (∆) [40] peaked at 1.7 MeV and with σ = 2.5 MeV, which approximates well the phenomenological distribution (see e.g. refs. [41,42]).

JHEP05(2019)173
In reality, only about 70% of beta transitions in reactors are of the Gamow-Teller type [43]. Most of the remaining ones are the first-forbidden transitions, whose neutrino spectrum has considerable uncertainties even in the SM limit, and whose dependence on non-standard interactions is poorly known (see ref. [44] for recent work in this direction). These are expected to give a non-negligible contribution to the reactor antineutrino flux, especially for E ν far above the detection threshold [45]. In particular, the first-forbidden decays may reintroduce some sensitivity to the pseudoscalar interactions. In this paper we ignore this complication, however we will check the robustness of our results by testing how much they rely on the events at the high end of the reactor antineutrino spectrum.
In the SMEFT approach one assumes no new degrees of freedom beyond those of the SM, therefore the sum in eq. (3.1) goes over the 3 neutrino states, and the oscillation probability in general depends on the two independent mass squared differences ∆m 2 21 and ∆m 2 31 . However, in short-baseline neutrino experiments one can typically neglect the effects proportional to ∆m 2 21 L/E; in particular, this is a good approximation in Daya Bay and RENO. In such a case eq. (3.1) simplifies to For the reactor experiments the relevant observable is the electron antineutrino survival probability. Taking α = e, and plugging in the expression of C e JK in eq. (3.4) we obtain where we defined the following combinations of the PMNS and NSI parameters: Using the detection and production coefficients in eqs. (3.6) and (3.7), the survival probability takes the form The oscillation formula in eq. (3.12) is valid away from ∆m 2 31 L/E ν ≈ 0, when the Wilson coefficients X obey the SMEFT scaling: 2 θ 13 . In its derivation we assumed that the off-diagonal elements of R vanish, which is true up to O(Λ −4 ) corrections when the SMEFT is a valid effective theory in some energy regime above m W (see eq. (A.1) for a more general formula). The expression for P νe→νe would be analogous with the reversed sign of the second line of eq. (3.12). Our formula agrees with the survival probability written down in ref. [3] after expressing their effective couplings s,d by the Wilson coefficients of the WEFT Lagrangian, see appendix B.
There are several important conclusions one can draw from eq. (3.12): • As mentioned before, the neutrino survival probability at the leading order depends only on off-diagonal Wilson coefficients X . We remark that the total number of produced and detected events (rather than the survival probability) is in principle sensitive to the diagonal scalar and tensor [ S,T ] ee , which we discuss in more detail in section 4.1. However, this caveat has no practical consequences due to the very stringent model-independent constraints on these coefficients from nuclear and meson decays [33]. As discussed below eq. (2.7), the effects of [ L ] ee and [ R ] ee are completely absorbed into the phenomenological values of the CKM element V ud and the axial charge of the nucleon g A , and are unobservable in neutrino oscillation experiments.
• The sensitivity of reactor experiments to pseudoscalar NSI ( P = 0) vanishes in the zero-recoil limit of beta decays, and when first-forbidden transitions in the reactor are neglected.
• At the leading order, reactor experiments alone are not sensitive to off-diagonal NSI of the V-A type ([ L ] eα = 0). The reason is that, as evident in eq. (3.13), their effects can be fully absorbed into a redefinition of the PMNS mixing angle θ 13 into the effective mixing angleθ 13 . 4 That redefinition can in fact be performed including also quadratic corrections in L [5]. Since θ 13 is an unknown parameter, which in the standard context was actually measured by Daya Bay, RENO, and Double Chooz, these experiments cannot separate the effect of the PMNS mixing parametrized by θ 13 from the new physics corrections contained in [ L ] eα . To that end, it is necessary JHEP05(2019)173 to measure another observable that is sensitive to a different combination of θ 13 and L than the one defined byθ 13 . This conclusion continues to hold when subleading terms in ∆m 2 21 are taken into account in the survival probability.
• On the other hand, reactor experiments are sensitive to scalar and tensor chargedcurrent interactions between leptons and quarks. The survival probability depends on the real and imaginary parts of the [S] and [T ] combinations defined in eq. (3.11). Two handles allow us to explore that dependence in practice. One, in the presence of CP violation (due to δ CP in the PMNS matrix or imaginary components of S,T ), the survival probability acquires a different oscillatory dependence on L/E ν than in the standard case. Secondly, the knowledge of the dependence of the survival probability on neutrino energy E ν allows one to disentangle CP-conserving effects of scalar and tensor interactions from each other, and from the (energy-independent) effective mixing angleθ 13 .
• The survival probability in eq. (3.12) manifestly satisfies 0 < P (ν e →ν e ) ≤ 1 in its regime of validity specified below eq. (3.13). Naively, for ∆m 2 31 L/E ν 1 one could obtain P (ν e →ν e ) > 1 or P (ν e → ν e ) > 1 (depending on the sign and magnitude of β X ) due to the contribution in the second line in eq. (3.12). In this regime, however, one can show that the O( 2 X ) contributions cannot be neglected; including the full non-linear X dependence in eq. (3.4) one recovers P (ν e →ν e ) ≤ 1 independently of the magnitude of β X . Note that the Daya Bay and RENO experiments are designed such that ∆m 2 31 L/E ν ∼ 1 for typical E ν , therefore this caveat has no practical consequences for our analysis. Note also that there are no non-oscillatory terms in eq. (3.12), therefore the so-called zero-distance effects [3,48] sometimes discussed in the NSI literature are absent in our approach. This is reassuring, as zero-distance effects at the linear level in X would also lead to P (ν e →ν e ) > 1 for some parameter choices.
• The last term in the survival probability in eq. (3.12) is proportional to sin(∆m 2 31 L/(2E ν )), which clearly depends on the choice of mass ordering. Throughout this analysis we assume ∆m 2 31 > 0. Choosing the inverted mass ordering would result in the opposite signs for the best fit values of Im [S] and Im [T ] compared to that determined in the next section. the scalar and tensor interactions compared to the SM one. Note that, as can be seen from the right panel of figure 1, reactor neutrino oscillations are more sensitive to tensor interactions. This is because they interfere with the SM axial interactions, which typically give larger contributions than the SM vector ones in reactor transitions and inverse beta decay. In the scalar case (the left panel) the far to near probability is similarly sensitive to the change in the real and imaginary parts. This is because the contributions of the real and imaginary parts consist of single terms (coming only from the detection side) which are of the same order and have similar effects on the probability. This is not the case for the tensor case (right panel), for which the real and imaginary parts contributions appear as the sum of two terms: those two terms have opposite sign for the real part but the same sign for the imaginary part. For this reason the survival probability is more sensitive to the imaginary part of [T ]. These comments are illustrative, as they they are valid for a fixed value ofθ 13 , and they can change in a complete analysis where the latter is also a floating parameter, as we will see in section 4.3.

Observables and NSI sensitivity
Typical reactor experiments detect antineutrinos via the inverse beta decay (IBD) process νp → e + n. They can measure not only the number of events but also the antineutrino energy. In the EFT framework the number of detected IBD events with an antineutrinoenergy E ν at a distance d is given by where P is the survival probability given in eq. (3.12), and dN no−osc ν is the differential number of IBD events that would be detected in the absence of oscillations. 5 5 We note that zero-distance effects are included in P and not in dN no−osc ν , since they are due to a mismatch between the source and detector neutrino flavor eigenstates. We recall that these effects do not appear at O(Λ −2 ).

JHEP05(2019)173
Note that the latter depends on the distance in a purely geometric form (∼ 1/d 2 ). It also depends on the nonstandard EFT coefficients through the nuclear decays widths and the IBD detection cross section. However, the situation is the opposite as in the survival probability, since N no−osc ν does depend linearly (order Λ −2 ) on flavor-diagonal coefficients [ X ] ee , but not on flavor-nondiagonal ones. Moreover, since there is no dependence at any order in the vector EFT couplings [ L,R ] ee (see discussion at the end of section 2.3), the only linear BSM corrections in N no−osc ν comes from flavor-diagonal scalar and tensor coefficients, [ S ] ee and [ˆ T ] ee . We have taken into account once again that the pseudoscalar contribution vanishes in the non-relativistic limit and we neglect the contribution of forbidden decays.
To reduce systematic errors, reactor experiments use near and far detectors at different distances from the reactor sources. The ratio of the number of IBD events in the energy bin around E ν =Ē i ν in two such detectors is given by where the last approximation is valid for small enough energy bins ∆. In that case, the ρ contributions to the numerator and denominator cancel except for the geometric ddependence (flux), and only the EFT corrections entering via the survival probability Pν e→νe have an effect on the ratio T j/k i . Note however that for large energy bins, EFT corrections to ρ with an energy dependence different to the SM one (as is the case for scalar and tensor interactions) do not cancel. Obviously, this is also the case for the ratio of the inclusive number of events in two detectors, since one has to integrate over all energies. Finally, we note that there is an additional linear effect of flavor-diagonal interactions in the total number of events, which (i) cancels in the far/near ratios and (ii) suffers the large uncertainty of the total reactor flux.
This introduces a dependence of the inclusive detector rates (and their far/near ratios) on the diagonal coefficients [ S ] ee and [ˆ T ] ee . However, given the relatively large (at least percent level) uncertainties in reactor nuclear processes, and the strong (per-mille level or better) constraints from "cleaner" beta decays and other precision experiments on these diagonal coefficients [33], we will simply ignore this effect in our analysis. In fact, for this same reason, such diagonal EFT coefficients cannot explain the observed deficit of detected reactor antineutrino fluxes relative to the SM predictions [38,39], which is often referred to as the reactor antineutrino anomaly [49].

Setup and analysis
In our numerical analysis we use the results from the Daya Bay [21] and RENO [22] experiments with 1958 days and 2200 days of data taking, respectively. The Daya Bay experiment has 4 near detectors located at Experimental Halls 1 and 2 (EH1 and EH2) and 4 far detectors located at Experimental Hall 3 (EH3). The weighted distances from the reactor cores are respectively 516 m, 555 m and 1571 m. The RENO experiment has one near and one far detector located at 367 m and 1440 m, respectively.

JHEP05(2019)173
First we define the following χ 2 function that only uses spectral information is the ratio of far to near IBD events in the i-th bin of energy. For Daya Bay, since there are two sets of near detectors at different distances, one defines where ω EH1 = 0.05545 and ω EH2 = 0.2057 are the weights that sample the different fluxes of the different reactors in equal proportions to the two near experimental halls [50]. The statistical uncertainty δR F/N i is given by where N i,bkg is the background expected in each energy bin. The systematic uncertainties of the background, reactor flux, and efficiency are taken into account by the pull parameters b d , f r , and , respectively, which we take from the original Daya Bay and RENO publication [22,50]. The d and r indices refer to the different detectors and reactors. Finally, we construct χ 2 spectral separately for RENO and Daya Bay, and combine the two in our analysis.
In addition to the spectral information in χ 2 spectral , we also take into account the ratio of the total IBD rate measured in the near and far detectors of Daya Bay and RENO, following closely the method described in ref. [51]. For our analysis we simply sum the two likelihoods: χ 2 = χ 2 spectral + χ 2 rate .

Results
We are ready to extract constraints on the mixing angle θ 13 and NSI parameters appearing in eq. (3.12) from a combination of Daya Bay and RENO data. 6 In our analysis we do not treat ∆m 2 31 as a free parameter, but rather use the best-fit value ∆m 2 31 = 2.52×10 −3 eV 2 [1] assuming normal ordering. This is justified because that result is dominated by other oscillation experiments than the reactor ones, and the new physics effect on the best-fit value and error is expected to be negligible. We have checked the stability of our results by letting ∆m 2 31 vary within its 1σ uncertainty. Consider first the case when only left-handed NSI are present: L = 0, R,S,P,T = 0. This corresponds to vanishing α and β parameters in eq. (3.12) and the only free parameter that remains isθ 13 . As explained previously, the reactor experiments alone are not sensitive to new physics parametrized by [ L ] eµ and [ L ] eτ , as these parameters can be absorbed into the unknown mixing angle θ 13 , and only theθ 13 combination defined in eq. (3.13) is probed. After marginalizing the χ 2 with respect to all the pull parameters we find the following result (all the uncertainties in this section are 68% CL): sin 2 (2θ 13 ) = 0.0841 ± 0.0027. In this simple case, leaving ∆m 2 31 as a free parameters would have a negligible impact on the confidence interval.
Next, we allow the scalar NSI to be non-zero: S = 0, R,P,T = 0. This implies α P = β P = 0 in eq. (3.12), but now α D and β D can be non-zero, that is to say, NSI effects can appear at the detection side. Of course, now we cannot use the value forθ 13 in eq. (4.5), as it was obtained under assumption that S = 0. Instead, we need to derive simultaneous constraints onθ 13 and the combination of NSI parameters [S] defined in eq. (3.11). We

JHEP05(2019)173
consider three different cases: 1) only Re [S] is non-zero, 2) only Im [S] is non-zero, and 3) both are non-zero and independent. We present our results in figure 2. The orange, blue, and purple contours are the 1-, 2-, and 3-σ allowed regions, respectively. For the real part we see some degeneracy betweenθ 13 and Re [S]. This is expected: while the two lead to a different energy dependence of the survival probability in eq. (3.12), they carry the same oscillatory dependence on L/E ν . Setting Im [S] = 0 (upper left panel in figure 2), we find the following constraint Re [S] = 0.54 ± 0.39 .  We emphasize that the sign of the best fit value for Im [S] depends on choosing the mass ordering, and would be flipped for the inverted ordering. Finally, we allow tensor NSI to be non-zero:ˆ T = 0, R,S,P = 0. The effects of tensor NSI appear on both the production and detection sides. We consider again three distinct cases, one where only Re [T ] is non-zero, another where only Im [T ] is non-zero, and where both are free parameters. Before presenting our results we recall that for calculating the production coefficients in eq. (3.7) we have assumed that the beta transitions in the reactors are of the Gamow-Teller type, while in fact almost 30% of the decays are first-forbidden transitions [43]. These transitions are expected to be more important at the high end of the neutrino spectrum [45]. Therefore, to test the robustness of our conclusions, we compare the results obtained using the whole neutrino spectrum with the ones where the neutrino energies are restricted to E ν < 5 MeV. We show the results in figure 3, using the same color coding as in figure 2. For the 3σ regions, we show the results using the entire neutrino energy spectrum (solid contours) and with the 5 MeV cut (dashed contours). The cut has limited impact on the preferred parameter regions, which suggests that the presence of forbidden transitions in the reactors should not affect our EFT constraints significantly.    23 [ P ] eτ ) of the pseudo-scalar couplings in the effective Lagrangian. However, since the effects of P are velocity-suppressed in our approximation, our analysis would be sensitive only to [P ] 1, outside the validity range of the effective theory. As we discuss below in section 5.3, much stronger constraints on P can be derived from pion decays.
It is worth mentioning that our constraints on Im [S] and Im [T ] are dominated by χ 2 spectral , with χ 2 rate having a small impact on the confidence intervals. On the other hand, using only the spectral information we find that the degeneracy between Re [S] andθ 13 is worsened, which translates in weaker marginalized bounds. Last, we also note that the O(0.1 − 0.4) bounds on Re [S] and Re [T ] obtained above do not do justice to the sensitivity and potential of these measurements. One should keep in mind the large correlation with θ 13 shown in figure 2. This translates to much more precise measurements in (i) less general scenarios, like the SM case in eq. (4.5) shows; or (ii) after combination with other measurements that are sensitive to the same coupling with different correlation.

Non-oscillation constraints on EFT parameters
In the previous section we derived a-few-percent-level constraints on linear combinations of Wilson coefficients [ S ] eα and [ˆ T ] eα , α = µ, τ , from neutrino oscillations in reactor experiments. To see these results in a wider context, in this section we discuss precision observables that do not involve neutrino oscillations but are sensitive to the same parameters. There is one important difference between these two classes. While the oscillations are sensitive to linear effects [ X ] eα , the observables discussed below are sensitive to absolute values squared of these parameters (or their combinations). One consequence is that they cannot distinguish between real and imaginary parts. Furthermore, the dependence on [ X ] eα enters at O(Λ −4 ) in the SMEFT expansion. It is in principle possible that their effects cancel against linear effects in [ X ] ee , coming from dimension-6 or dimension-8 SMEFT operators. For illustration we neglect such terms in the discussion below and set bounds when only one [ X ] eα term is present in the Lagrangian at a time. Because of these JHEP05(2019)173 assumptions, the bounds obtained in this section are less robust than the ones from oscillations. Although they may be valid for the SMEFT derived from particular UV models, it is important to keep in mind that they can be relaxed significantly if several interactions are present at the same time. A thorough analysis of such scenarios is however beyond the scope of this work, which is focused on neutrino physics.

Neutron and nuclear beta decay
Instead of the plethora of β-decay transitions happening inside nuclear reactors, one can search for nonstandard effects in specific decays that happen to be very clean both experimentally and theoretically [33]. One expects strong bounds on nonstandard interactions involving wrong-flavor neutrinos from such studies, which has been used sometimes in the past as an argument to neglect e.g. scalar and tensor interactions. However, to best of our knowledge, all available beta-decay analysis have focused on interactions involving electron neutrinos. We amend this situation here, deriving the bounds on scalar and tensor operators with a wrong-flavor neutrino. Our experimental input are the so-called F t values of 0 + → 0 + transitions [52] and neutron data (lifetime and correlation coefficients). We use the same statistical approach and dataset as in the recent review in ref. [33] (table 4, 5 and 7), but also including the recent PERKEO-III measurement of the beta asymmetry in neutron decay, A n = −0.11985(21) [53]. We note that the errors of the average lifetime and beta asymmetries are re-scaledà la PDG to take into account tensions among various measurements. Assuming only one interaction is present at a time we find the following 90% CL bounds: We note that the observables depend quadratically on these WEFT coefficients and thus the error distribution is highly non-gaussian.

CKM unitarity
In the beta decay fit discussed above one extracts simultaneously the non-standard scalar or tensor coupling and the SM parameters (|V ud | and the axial charge g A ). It is important to note that significant correlations between the scalar coupling and |V ud | appear. Thus, adding to this analysis the very precise |V ud | value obtained from CKM unitarity: |V unit. ud | ≡ 1 − |V us | 2 − |V ub | 2 1/2 , has a drastic impact on the bound on scalar interactions [33]. Namely: where we used V us = 0.2243(5) and V ub = 0.00394(36) [54]. Roughly speaking, the bound above comes from the comparison of |V ud | 2 1 + g 2 S |[ S ] eα | 2 , extracted from 0 + → 0 + transitions, and |V unit. ud | 2 . The former happens to be currently a bit smaller than the latter, which translates in a bound on |[ S ] eα | more stringent than naively expected (because this interaction can only contribute positively). 7

JHEP05(2019)173
CKM unitarity constrains also the offdiagonal vector coefficients, [ L,R ] eα . In a oneoperator analysis we find |[ L,R ] eα | ≤ 1.9 × 10 −2 (90% CL) . (5. 3) The stringent bounds in eq. (5.2) and eq. (5.3) assume (i) the absence of other nonstandard β-decay couplings, as in the previous section; (ii) the absence of new physics effects in the extraction of V us , V ub , and the Fermi constant G F ; and (iii) the 3-family setup, which is not an extra assumption in the SMEFT.

Leptonic pion decays
The π → eν e channel is extremely sensitive to pseudo-scalar couplings because the latter do not suffer the strong chiral suppression of the SM contribution (the SM width vanishes for zero electron mass). For a pseudo-scalar interaction with a wrong-flavor neutrino, a bound on |[ P ] eα (µ = 2 GeV)| can be derived from the clean ratio R π of the π → eν e and π → µν µ widths [58][59][60]. Using the experimental and SM values R π = 1.2327(23) × 10 −4 [54] and R SM π = 1.2352(1) × 10 −4 [60], the 90% CL constraint is given by BSM models generating scalar/tensor interactions often generate pseudoscalar interactions of similar magnitude. More importantly, even if this is not the case at tree-level, the pseudoscalar interactions is generated radiatively [61]. For instance, the connection between P (2 GeV) and the coefficients at the EW scale and 1 TeV are given by [30]  The larger mixing found in the 1-TeV case is due not only to the trivial larger running but also because the mixing happen to be larger in the SMEFT than in the WEFT [30]. In full generality P (2 GeV) represents a different direction in the parameter space with respect to S,T (2 GeV). However, the mixing relations given above imply strong constraints on the tensor coupling in simplified scenarios where (pseudo)scalar and tensor couplings are not independent degrees of freedom, since severe cancellations among them are not possible anymore. For example, if we assume that the pseudo-scalar operator is not generated at tree-level at the high-scale, we obtain the following 90% CL bounds where the ranges correspond to the one-operator and global analysis discussed in eq. (5.4). Let us note that the derivation of these bounds only takes into account log-enhanced oneloop corrections and it can be altered by finite pieces, especially if the running is not carried out to very high-energy scales.

LHC (pp → e + MET + X)
One can look for the same nonstandard charged-current interactions (or more precisely, for their SMEFT counterparts) in the Drell-Yan process pp → e + MET + X [28,62]. This connection requires additional assumptions such as the validity of the SMEFT at such high-energies, approximating V ij = δ ij in the SMEFT-WEFT mapping, and, especially, neglecting the contributions from dim-8 operators (since LHC bounds are dominated or very sensitive to dim-6 squared contributions). Chirality-flipping interactions do not interfere with the SM and then the usual bounds on operators involving electron neutrinos actually to the incoherent sum over all three flavor neutrinos. Thus we can reinterpret the results from figure 8 of ref. [63], which used the 13-TeV ATLAS search with 36 fb −1 [64]: at 90% CL and at µ = 2 GeV.

Charged-lepton-flavor violation
In the SMEFT, dimension-six operators that give rise to charged-current interactions between quarks and leptons also yield neutral-current interactions between quarks and pairs of charged leptons = e, µ, τ . In consequence, neutrino interactions parametrized by offdiagonal [ X ] eα appear in the Lagrangian together with 4-fermion charged-lepton-flavor violating (CLFV) interactions. The latter mediate at tree-or loop-level such processes as → γ, → 3 , or N → N , which for = have not been observed so far and are stringently constrained by experiment. The contribution of dimension-6 operators to CLFV observables arises at O(Λ −4 ), which is the leading order in this case because the SM contributions are absent. The resulting constraints on lepton-flavor off-diagonal SMEFT operators are typically very severe [65][66][67][68]. Using the analytical formula for the µ → e conversion rate in ref. [69] and the experimental bound Br(µ → e) Au ≤ 7 × 10 −13 [70] WEFT were not UV-completed by the SMEFT (because new physics is not much heavier than m Z , or because electroweak symmetry is realized non-linealy in the UV theory), as then off-diagonal X are not correlated in general with CLFV interactions.

Conclusions
We have proposed a systematic approach to neutrino oscillations in the situation when neutrino interactions with matter are modified by heavy new physics. To this end we employed the model-independent framework of the SMEFT, with the Lagrangian organized into an expansion in powers of 1/Λ, where Λ is the mass scale of new particles affecting the neutrino interactions. The SMEFT framework enables consistent power-counting of new physics effects and identifying the leading order corrections to the neutrino oscillations probability. In this paper we applied it to oscillations in short-baseline reactor experiments, however the formalism can be readily extended to other types of neutrino experiments.
We calculated the survival probability of reactor antineutrinos at O(Λ −2 ) in the SMEFT expansion, that is including linear effects of dimension-6 operators. The main result of this paper is given in eq. (3.12), from which the dependence of the survival probability on the PMNS parameters and dimension-6 Wilson coefficients can be read off. We have taken into account all SMEFT operators that contribute at the leading order. In addition to operators affecting the SM-like (V-A) charged-current interactions between quarks and leptons, those mediating scalar and tensor contact interactions contribute at the same order Λ −2 . The latter lead to a different energy dependence of the neutrino production and detection amplitudes, which is reflected in eq. (3.12). We also paid due attention to the interplay between the effects of dimension-6 operators and of the PMNS mixing. It is a familiar fact in electroweak precision measurements and flavor physics that some dimension-6 corrections to physical observables can be absorbed into SM parameters, such as the Higgs vacuum expectation value or the Wolfenstein parameters. In the present case one observes an analogous effect. It turns out that corrections to the survival probability due to lepton-flavor off-diagonal V-A interactions parametrized by [ L ] eα can be absorbed into a redefinition of the mixing angle θ 13 [5]. Since θ 13 is not known a-priori (other than from the very reactor experiments we consider here), the existing data only constrain a linear combinationθ 13 of the original mixing angle θ 13 and dimension-6 Wilson coefficients, cf. eq. (3.13), but not the two separately. We conclude that, at the leading order in the SMEFT expansion, reactor neutrino experiments alone only constrain the scalar and tensor interactions, but not the V-A ones. This explains the origin of the degeneracy between θ 13 and off-diagonal V-A NSI found in the previous literature.
We pointed out that neutrino oscillations can probe, at the linear (order Λ −2 ) level, the dimension-6 tensor and scalar SMEFT operators that are off-diagonal in the leptonflavor space. To our understanding, it is the unique class of observables where this is the case. Consequently, the oscillation constraints on these operators are robust as long as the expansion of the SMEFT Lagrangian in powers of 1/Λ is quickly convergent. The same operators can be probed by meson and nuclear decays or by production processes at the LHC, however in those cases they enter quadratically (at order Λ −4 ). Such constraints are JHEP05(2019)173 then subject to model-dependent assumptions about other dimension-6 and dimension-8 contributions to the same observables, and are thus less robust.
We identified the linear combinations of Wilson coefficients of scalar and tensor SMEFT operators that can be constrained by reactor oscillation experiments, cf. eq. (3.11). Using the most recent data from Daya Bay and RENO, we derived numerical constraints on these combinations. Under various more or less constraining assumptions, they are presented in section 4 and illustrated in figures 2 and 3. As of today, constraints at a few percent level can be extracted from the publicly available reactor experiment data. This is competitive with the constraints extracted from nuclear decays (which are less robust, as discussed above). At face value, the LHC constraints on the same operators are at least an order of magnitude more stringent. However, they are more model-dependent, and rely on the assumption that the SMEFT is a valid framework at TeV energy scales. We also derived stringent constraints on scalar and, especially, tensor Wilson coefficients arising due to renormalization group mixing with the pseudoscalar ones. The latter are strongly constrained by pion decays thanks to chiral enhancement. Again, those constraints are less robust, in particular they depend on the starting point of the running, and assume the absence of cancellations between different contributions. Finally, CLFV processes typically place severe constraints on SMEFT operators that are off-diagonal in lepton-flavor indices. In particular, µ → e conversion on nuclei strongly constrains operators contributing to the NSI parameter [ S ] eµ , while τ → eπ + π − decays constrain the operators contributing to [ S ] eτ . To our knowledge, however, CLFV constraints on the operators contributing at tree level to the tensor parameters [ T ] eµ , and [ T ] eτ (that also affect reactor neutrino oscillations at the leading order) are not given in the literature.
We note that the main goal of reactor experiments so far has been a precise determination of the mixing angle θ 13 in the standard context, and the analyses were certainly not optimized for constraining SMEFT operators. We believe that with more data, more detailed spectral information, and targeted analyses, the constraints obtained in this paper can be significantly improved. Furthermore, in our analysis the constraining power of reactor experiments is weakened by a partial degeneracy between NSI and the effective mixing angleθ 13 . Designing new observables sensitive to a different linear combination ofθ 13 and [ X ] eα may be another path to increasing senisitivity to new physics.

A WEFT without SMEFT
The WEFT is an effective theory below the scale µ m W describing interactions of the SM particles with the exception of the W , Z and Higgs bosons and the top quark, which have been integrated out. In the main body of this paper we treated the WEFT as the low-energy theory of the SMEFT. However, the WEFT is a consistent EFT in its own right, which can be valid even if it is not completed by the SMEFT at higher energies. This caveat is relevant if the masses of BSM particles are between a few and a 100 GeV, or if the electroweak symmetry in the BSM theory is realized non-linearly. In this situation the leading order Lagrangian relevant for neutrino oscillations is still the one of eq. (2.1), however the parameters X are no longer related by matching to the SMEFT parameters at higher energies. The most important practical consequence is that R may contain non-diagonal elements at the leading order in the WEFT. Then the reactor antineutrino survival probability in eq. (3.12) is generalized to where the definitions for γ R and the (new) effective mixing angle arê In this setting the survival probability acquires dependence on one linear combination of the R parameters, which is in principle distinguishable from other parameters in eq. (A.1) due to the distinct L/E ν and E ν dependence. We show in figure 4 the results of the combined analysis of the Daya Bay and RENO experiments allowing only the righthanded NSI to be non-zero: R = 0, S,P,T = 0. Under this assumption we find

B Traditional NSI formalism
In the neutrino literature, NSI are typically parametrized by the effective couplings s and d , which correspond to non-standard effects in neutrino production and detection, respectively (see e.g. ref. [3]). In this approach, neutrinos produced at the source and detected at the detectors are not pure flavor states:  where the first and second indices in s αβ correspond to the flavors of the charged lepton and neutrino, respectively, which is reversed for d αβ . The oscillation probability is given by where in the absence of NSI in propagation the Hamiltonian is given by Up to first order in 's, rewriting s eα = | s eα |e iφ s eα and d βe = | d βe |e iφ d βe , and taking the limits ∆m 2 21 L/E ν 1, cos θ 13 ≈ 1, the survival probability of electron antineutrinos becomes [3,11]

JHEP05(2019)173
Comparing the probability in eq. (B.4) with the WEFT formula in eq. (A.1), the two agree given the matching between the effective couplings for α = µ, τ . On the other hand, matching eq. (B.4) to eq. (3.12) we find Re s ee + d ee = 0. In our formalism, the diagonal EFT coefficients do not enter the survival probability, but rather the total rate of produced and detected neutrinos (see section 4.1). This way, Pνs e →ν d e can indeed be interpreted as a probability, since non-oscillatory terms linearly proportional to s,d ee would lead to Pνs e →ν d e > 1 for some choices of parameters. It is important to note that here we compare the "traditional" NSI formalism only to short baseline experiments, while the matter effects would still need to be included for long base-line experiments.

C Separate RENO and Daya Bay analyses
For completeness, in this appendix we show separate RENO and Daya Bay constraints on the mixing angle and the NSI parameters in eq. (3.12). This allows us to compare the constraining power of the two experiments and their relative weight in the combined fit. We follow the same presentation as in section 4. We consider first the case when only left-handed NSI are present: L = 0, R,S,P,T = 0, in which case the only free parameter is θ 13 . It is constrained as