Probing Non-Standard Interactions at Daya Bay

In this article we consider the presence of neutrino non-standard interactions (NSI) in the production and detection processes of reactor antineutrinos at the Daya Bay experiment. We report for the first time, the new constraints on the flavor non-universal and flavor universal charged-current NSI parameters, estimated using the currently released 621 days of Daya Bay data. New limits are placed assuming that the new physics effects are just inverse of each other in the production and detection processes. With this special choice of the NSI parameters, we observe a shift in the oscillation amplitude without distorting the $L/E$ pattern of the oscillation probability. This shift in the depth of the oscillation dip can be caused by the NSI parameters as well as by $\theta_{13}$, making it quite difficult to disentangle the NSI effects from the standard oscillations. We explore the correlations between the NSI parameters and $\theta_{13}$ that may lead to significant deviations in the reported value of the reactor mixing angle with the help of iso-probability surface plots. Finally, we present the limits on electron, muon/tau, and flavor universal (FU) NSI couplings with and without considering the uncertainty in the normalization of the total event rates. Assuming a perfect knowledge of the event rates normalization, we find strong upper bounds $\sim$ 0.1\% for the electron and FU cases improving the present limits by one order of magnitude. However, for a conservative error of 5\% in the total normalization, these constraints are relaxed by almost one order of magnitude.


Introduction
The recent discovery of the smallest neutrino mixing angle θ 13 by the modern reactor antineutrino experiments Daya Bay [1][2][3] and RENO [4] has firmly established the three-flavor neutrino paradigm [5][6][7] and signifies an important development towards our understanding of the structure of the neutrino mass matrix, whose precise reconstruction would shed light on the underlying new physics that gives rise to neutrino mass and mixing [8][9][10]. Another reactor electron antineutrino disappearance experiment: Double Chooz [11,12], and the two accelerator electron (anti-)neutrino appearance experiments: MINOS [13] (completed) and T2K [14,15] (presently running) have also confirmed the non-zero and moderately large value of θ 13 in the standard three-flavor oscillation scenario. It is quite remarkable to see that with 621 days of data taking and using the merit of identical multi-detector setup, the Daya Bay experiment reveals a non-zero value of θ 13 at more than 16σ and suggests a best-fit value of sin 2 2θ 13 = 0.084 ± 0.005 [3]. This data provides a relative 1σ precision of 6% on sin 2 2θ 13 which is already better than the precision achieved on sin 2 θ 23 . Undoubtedly, this high precision measurement of the 1-3 mixing angle has speeded up the search for the neutrino mass ordering and the possible presence of a CP-violating phase in current and future neutrino oscillation experiments [16][17][18][19][20].
To explain the presence of small neutrino masses and relatively large neutrino mixings as indicated by neutrino oscillation data, various neutrino mass models have been proposed. These neutrino mass models come in various categories such as the cases where neutrinos acquire mass via the popular seesaw mechanism [21][22][23][24][25][26][27][28][29]. We find also models where neutrinos get mass radiatively due to the presence of extra Higgs bosons [30,31] or low energy supersymmetric hybrid models with spontaneous or bilinear breaking of R-parity [32,33]. The structure of the standard electroweak neutral and charged currents gets affected by the presence of these mechanisms responsible for the neutrino mass generation [25]. In most of the cases, in the low energy regime, these effects are known as non-standard interactions (NSI). Various extensions of the Standard Model (SM), such as left-right symmetric models and supersymmetric models with R-parity violation, predict NSI of neutrinos with other fermions [34][35][36][37][38][39][40][41][42]. The NSI in these models are usually generated via the exchange of new massive particles at low energies.
Neutrino NSI may be of charged-current (CC) or neutral-current (NC) type, and they can be classified in two main categories: flavor-changing NSI, when the flavor of the leptonic current involved in the process is changed, or flavor-conserving non-universal NSI, when the lepton flavor is not changed in the process but the strength of the interaction depends on it, violating the weak universality. In the low energy regime, these new interactions may be parameterized in the form of effective four-fermion Lagrangians: where G F is the Fermi constant, α and β are neutrino or lepton flavor indices, f and f label light SM fermions, and the dimensionless coefficients ε parametrize the strength of the interaction. Here we denote the CC-NSI couplings as ε s d,f f αβ since they affect in general the source (s) and detector (d) interactions at neutrino experiments, while ε m,f f αβ refers to the NC-NSI couplings generally affecting the neutrino propagation in matter (m). Note that in Eqs. (1.1) and (1.2) we have assumed for the new interactions the same Lorentz structure as for the SM weak interactions, (V±A). Even though more general expressions are possible within a generic structure of operators, as shown in Ref. [43], these are the dominant contributions for reactor experiments, where we will focus our attention in this work.
NSI effects may appear at three different stages in a given neutrino experiment, namely neutrino production, neutrino propagation from the source to the detector, and neutrino detection. In a short-baseline reactor experiment, the effects on the neutrino propagation are negligible since it happens mainly in vacuum. Therefore, the new generation of short-baseline reactor experiments, such as Daya Bay, offers an excellent scenario to probe the presence of NSI at the neutrino production and detection, free of any degeneracy with NSI propagation effects. Actually, some work has already been done in the context of NSI at short-baseline reactor experiments. In particular, forecasts for the sensitivity of Daya Bay to NSI have been published, for instance, in Ref. [44]. More recently, another article presented constraints on neutrino NSI using the previous Daya Bay data set [45]. The results derived there are in general agreement with some of the cases we discuss in Sec. 4. Nevertheless, our paper also studies the phenomenology of some other interesting cases not considered in Ref. [45], and provides a detailed description of the effect of NSI in the neutrino survival probability. Finally, we also discuss the fragility of the bounds on the NSI couplings derived using Daya Bay data against the presence of an uncertainty on the total event rate normalization in the statistical analysis of reactor data.
Given the production and detection neutrino processes involved in short-baseline reactor neutrino experiments (β-decay and inverse β-decay), the NSI parameters relevant for these experiments are ε ud eα , i.e., the CC-NSI couplings between up and down quarks, positrons and antineutrinos of flavor α. In the literature we can find the following 90% C.L. bounds on these parameters [46]: coming from unitarity constraints on the CKM matrix as well as from the non-observation of neutrino oscillations in the NOMAD experiment. Here the couplings with different chirality are related by: Other constraints on neutrino NSI couplings using solar and reactor neutrinos have been given in [47][48][49][50][51][52]. On the other hand, NSI have also been studied in the context of laboratory experiments with accelerators in Refs. [50,[53][54][55], while bounds from atmospheric neutrino data have been presented in Refs. [56,57]. Recently a forecast for the sensitivity to NSI of the future PINGU detector has been presented in [58]. This paper is organized as follows. In Sec. 2, we describe the procedure for implementing the NSI effect in the modern reactor experiments. There, we derive the effective antineutrino survival probability expressions which we use later to analyze the Daya Bay data. We also give plots to discuss in detail the impact of the NSI parameters on the effective probability considering a special case where ε s eγ = ε d * γe . We also show the possible correlations between the NSI parameters and θ 13 with the help of iso-probability surface plots. Sec. 3 describes the numerical methods adopted to analyze the reactor data. Apart from this, a brief description of the Daya Bay experiment and the important features of its present data set which are relevant for the fit are also given in this section. Sec. 4 presents the constraints on the NSI parameters imposed by the current Daya Bay data assuming a perfect knowledge of the event rates normalization. Next we derive the bounds on the NSI parameters taking into account the uncertainty in the normalization of event rates with a prior of 5% in Sec. 5. In Sec. 6, we compare the constraints on the NSI parameters obtained 1 The NSI parameters probed in our analysis get contributions from the (V±A) operators in Eq. (1.1) and therefore they can be generally expressed as ε ud,V±A αβ or equivalently ε ud,L±R αβ . However, for simplicity, we have dropped the chirality indices all over the paper.
with the current 621 days of Daya Bay data with the limits derived using the previously released 217 days of Daya Bay data. Finally, we summarize and draw our conclusions in Sec. 7. In Appendix A, we give the effective antineutrino survival probability expressions for the physical situations where ε s eγ = ε d * γe .

Implementing NSI in Modern Reactor Experiments
According to the usual procedure followed in the non-standard analyses of reactor data [44,59], we start by re-defining the neutrino flavour states in the presence of NSI in the source and detection processes. For the initial (at source) and final (at detector) neutrino flavor states, we have [60][61][62][63]: while the redefinition of the antineutrino flavor states is given by: The normalization factors required to obtain an orthonormal basis can be expressed as: 3) and the neutrino mixing between flavor and mass eigenstates is given by the usual expressions: The correct normalization of the neutrino states in presence of NSI is a very important point, required to obtain a total neutrino transition probability normalized to 1. However, one has to consider that when dealing with a non-orthonormal neutrino basis, the normalization of neutrino states will affect not only the neutrino survival probability but also the calculation of the produced neutrino fluxes and detection cross sections. In this case, as shown in Ref. [64], all the normalization terms coming from N s α and N d β will cancel while convoluting the neutrino oscillation probabilities, cross sections, and neutrino fluxes to estimate the number of events in a given experiment such as Daya Bay. This is due to the fact that the SM cross sections and neutrino fluxes used in our simulation have been theoretically derived assuming an orthonormal neutrino basis and therefore they need to be corrected. Then, from here we can consider the following effective redefinition of neutrino and antineutrino states: where we have dropped the normalization factors that will cancel in the Monte Carlo simulation of Daya Bay data. From these effective neutrino states we will calculate an effective neutrino oscillation probability, that will be used all along our analysis. Note that when we will discuss the features of the probability prior to the simulation of a particular experiment, we will always refer to the effective probability, that might be greater than one.

Effective antineutrino survival probability in reactor experiments
The effective antineutrino transition probability from flavor α to β after traversing a distance L from source to detector is defined as: In terms of the neutrino mass differences and mixing angles, this transition probability in vacuum may be written as: In the case of standard oscillations, Y j αβ is defined as: In presence of NSI, however, according to the definition of neutrino states in Eq. (2.2), this expression is modified as follows [59]: To obtain theν e survival probability in a reactor experiment, where an electron antineutrino is produced at the source and a positron is detected inside the detector, one has to replace α and β by e in Eq. (2.8). For the NSI parameters, we adopt the following parametrization by splitting the new couplings into its absolute value and its phase: The linear coefficients of the terms of order |ε| in Eqs. (2.14, 2.15, 2.16) are the same as given in Ref. [43], with a good agreement between the calculated probabilities here and there. In Eq. (2.15), note the presence of a term linear in the sine of ∆m 2 31 L/2E which therefore depends on the choice of neutrino mass ordering. This term does not appear in the standardν e →ν e oscillation expression and it can affect the L/E dependence of the probability in the presence of neutrino NSI.

Special case of NSI
In this work 2 we assume that, likewise the mechanisms responsible for production (via β-decay) and detection (via inverse β-decay) of reactor antineutrinos are just inverse of each other, this is also true for the associated NSI [43,44]. This assumption allows us to write ε s γ = ε d * γ ≡ |ε γ |e iφγ where we drop the universal e index for simplicity. With these assumptions, Eq. (2.12) takes the form (keeping the terms up-to the second order in small 2 In the appendix, we have given the effective probability expressions for the physical situations where ε s eγ = ε d * γe . In such cases, the spectral analysis of the reactor data becomes inevitable since the NSI parameters not only cause a shift in θ13 i.e. the change of the depth of the first oscillation maximum but also modify the shape of the probability due to the shift in its energy. A detailed analysis of the Daya Bay data under such scenarios will be performed in [65]. and, second, there is a shift in the effective 1-3 mixing angle due to oscillatory NSI terms which can be written as (2.19) These two features, which are brought about by the NSI parameters, are responsible for a global shift of the oscillation amplitude without affecting the shape of the oscillation probability as can be clearly seen from Fig. 1 and Fig. 2 that we will discuss in the next section. Eq. (2.19) suggests that it will be quite challenging to discriminate the effect of true θ 13 and NSI parameters in the modern reactor experiments. It is also interesting to note that there are some CP conserving terms in Eq. (2.19) which come into the picture due to the presence of NSI parameters. One of the most important consequences of the new non-oscillatory NSI terms (see Eq. (2.18)) is that they can cause a flavor transition at the source (L = 0) even before neutrinos start to oscillate. In the literature, this feature is known as "zero-distance" effect [43,66]. In modern reactor experiments, this effect can be probed using the near detectors which are placed quite close to the source. For definiteness, in this work we have restricted our analysis to the following choices of the NSI parameters: • Lepton number conserving non-universal NSI parameters which depend on the flavor characterizing the violation of weak universality. Under this category, we study the following two cases: 1. Considering only the NSI parameters |ε e | and φ e which are associated withν e .
At first order in |ε e | and neglecting the effect of the solar mass splitting, the new non-oscillatory NSI terms appearing at the survival probability produce a total shift in the effective θ 13 mixing angle given by: This expression will be very useful to discuss the behavior of the effective probability as well as the correlations between θ 13 and the NSI parameters in the next subsections.
In the presence of these flavor violating NSI parameters, Eq. (2.17) takes the form: Note that, for the NSI parameters |ε τ | and φ τ which are associated withν τ , the effective survival probability will be exactly the same as Eq. (2.22) with the replacements |ε µ | → |ε τ | and φ µ → φ τ provided that the 2-3 mixing angle is maximal ,i.e., sin 2 θ 23 = 0.5. As before, the new non-standard oscillatory terms in the neutrino transition probability may be interpreted as a global redefinition of the effective θ 13 mixing angle: • Lepton number conserving universal NSI parameters which do not depend on flavor.
In this case, we have |ε e | = |ε µ | = |ε τ | = |ε| and φ e = φ µ = φ τ = φ and the probability in Eq. (2.17) takes the form: In this case, the effective mixing angle in the presence of oscillatory and non-oscillatory NSI terms will be given by:

Impact of the NSI parameters on the effective probability
Now we will study the possible impact of the NSI parameters at the effective probability level. Table 1 depicts the benchmark values of the various oscillation parameters that are considered to generate the oscillation probability plots. These choices of the oscillation parameters are in close agreement with the best-fit values that have been obtained in the recent global fits of the world neutrino oscillation data [5][6][7]. Here we would like to mention that for the 2-3 mixing angle, we have taken the maximal value i.e. sin 2 θ 23 = 0.5, though in the global fit studies, there is a slight hint for a non-maximal value of θ 23 . We have also assumed normal mass ordering, i.e., ∆m 2 31 positive. In Fig. 1, we present the standard and the NSI-modified three-flavor oscillation effective probability as a function of the electron antineutrino energy with a source-detector distance of 1.58 km. In the left panel of Fig   the contribution from the non-oscillatory terms takes the form 4|ε e | + 6|ε e | 2 which always gives an overall positive contribution to the full probability for any choice of |ε e |. We can also see that even for a small value of |ε e | of 0.02, the effective oscillation probability can be more than unity for most of the energies of interest. This is the sign of the nonunitarity effects [64,[68][69][70][71], caused by the presence of neutrino NSI at the source and detector of reactor experiments. In the right panel of Fig. 1, the band shows the changes in the effective probability after varying |ε µ,τ | in the range [0, 0.15] and (δ − φ µ,τ ) over the range [-180 • , 180 • ] simultaneously. Note that in Eq. (2.22), the phases appear in the form of cosine of (δ − φ µ,τ ) and also the NSI terms have a non-trivial L/E dependency. The standard oscillation probability without the NSI terms is shown by the solid black line and the other four lines have been drawn considering particular choices of |ε µ,τ | and (δ − φ µ,τ ) which are mentioned in the figure legends. If (δ − φ µ,τ ) = 0 • (180 • ), then the effective oscillation probability is less (more) compared to the standard value for almost all the choices of neutrino energy as opposed to the case of the NSI parameters associated withν e .    2 shows the impact of the NSI parameters at the probability level with L = 1.58 km for the flavor-universal NSI case where we consider |ε e | = |ε µ | = |ε τ | = |ε| and φ e = φ µ = φ τ = φ. In this plot, the dark salmon region has been generated by varying the NSI parameter |ε| in the range [0, 0.05] and φ in the range [-180 • , 180 • ] simultaneously, keeping the CP phase δ fixed to 0 • . Next, we vary δ in its entire range from -180 • to 180 • along with the NSI parameters |ε| and φ and obtain the extended probability band in the form of the light grey region. In Fig. 2, the solid black line depicts the standard probability without considering the NSI parameters. The other four lines in this plot display the effective oscillation probability for particular combinations of |ε| and φ with δ = 0 • which are mentioned in the figure legends. It is quite clear from Eq. (2.24) that the non-oscillatory terms dominate over the oscillatory terms in the flavor-universal NSI case. Therefore, the dark salmon region of Fig. 2 closely resembles the left panel of Fig. 1 where we consider the NSI parameters which are only associated withν e , namely |ε e | and φ e . Next, we discuss the possible correlations between the NSI parameters and θ 13 with the help of iso-probability plots.

Correlations between NSI parameters and θ 13 : iso-probability plots
We consider the neutrino energy E = 4 MeV and the source-detector distance L = 1.58 km to draw the iso-probability surface plots. Left panel of Fig. 3 shows the iso-probability surface contours in the (sin 2 θ 13 -|ε e |) plane for four different choices of φ e considering sin 2 θ 13 = 0.023 and |ε e | = 0 as best-fit choices. Our best-fit choices correspond to the standard oscillation probability without considering the NSI parameters as given by Eq. (2.13). Now as we consider the finite value of |ε e | with φ e = 0 • (see the solid red line), the NSI terms  .20)). Then, we need to increase the value of sin 2 θ 13 to reduce the SM probability in order to compensate the enhancement due to the NSI contribution. For φ e = 180 • case (see the dot-dashed magenta line), a finite value of |ε e | decreases the overall probability demanding a lower value of sin 2 θ 13 so as to enhance the contribution from the standard probability. In the cases with φ e = 90 • or -90 • (see the dashed blue line), the non-oscillatory NSI term which is linear in |ε e | drops out from the probability expression and the remaining NSI contribution is 2|ε e | 2 (see Eq. (2.20)). Due to this weak quadratic dependence on |ε e |, a very large value of |ε e | is needed to compensate even a very small increment in sin 2 θ 13 . Right panel displays the same in the (sin 2 θ 13 -|ε µ,τ |) plane for different choices of (δ − φ µ,τ ). Here sin 2 θ 13 = 0.023 and |ε µ,τ | = 0 are considered as true choices. In this panel, the combination of phases (δ − φ µ,τ ) shows opposite features for 0 • and 180 • as compared to φ e in the left panel. Note that the impact of |ε µ,τ | on sin 2 θ 13 is weaker compared to |ε e |. If we examine the terms which are linear in |ε e | (see Eq. (2.20)) and |ε µ,τ | (see Eq. (2.22)) then we can see that the contribution coming from |ε µ,τ | is sin θ 13 suppressed even if we work at the first oscillation maximum. For (δ − φ µ,τ ) = 90 • or -90 • , the |ε µ,τ |-dependent terms completely disappear from Eq. (2.22) if sin 2 θ 23 = 0.5 and sin 2 ∆ 31 = 1. Therefore, we do not see any correlation between sin 2 θ 13 and |ε µ,τ | for (δ − φ µ,τ ) = 90 • or -90 • in the right panel of Fig. 3. In Fig. 4, we show the correlation between sin 2 θ 13 and |ε| for the flavor-universal NSI case with the help of iso-probability surface contours. In the left panel, we consider four different choices of the CP phase δ keeping φ fixed to 0 • . In the right panel, we take five different choices of φ assuming δ = 0 • . In both panels, the iso-probability surface contours are the same for 90 • and -90 • choices of phases (see the dashed blue lines) because the phases appear in the form of cosines in Eq. (2.24). For the φ = 0 • case (left panel), the contribution from the non-oscillatory NSI terms is maximum and takes the form 4|ε| + 10|ε| 2 . Thus, it enhances the overall oscillation probability to a great extent even for a small value of |ε|. Now to compensate this enhancement, we need to increase the value of sin 2 θ 13 to reduce the contribution coming from the standard survival oscillation probability. Right panel of Fig. 4 (δ = 0 • case) closely resembles the left panel of Fig. 3 because for the flavor-universal NSI case, the contribution from the non-oscillatory NSI terms dominates which is also true for the NSI parameters associated withν e and also these NSI parameters are δ independent. In the right panel, we present a special case of φ = 108 • to explain the features emerging from the extreme right panel of Fig. 7 (see later in Sec. 4) where we have displayed the allowed region in (sin 2 θ 13 -|ε|) plane using the current data from the Daya Bay experiment allowing the flavor-universal NSI phase φ to vary in the entire range of [-180 • , 180 • ] with the CP phase δ to be fixed to 0 • . Eq. (2.24) suggests that we can always choose some values of φ such that the non-oscillatory terms are canceled. This happens, for instance, for φ = 108 • and |ε| = 0.206. In this case, the dominant oscillatory NSI terms give a negative contribution to the neutrino survival probability that, to be compensated, requires an enhancement of the standard survival probability by reducing the value of sin 2 θ 13 by the same quantity. This is exactly the feature that we can see for φ = 108 • case.

Data analysis of modern reactor experiments
Reactor antineutrinos are produced by the fission of the isotopes 235 U, 239 Pu, 241 Pu and 238 U contributing to the neutrino flux with a certain fission fraction f k . Reactor antineutrinos are detected via inverse β-decay process (IBD),ν e + p → e + + n. The technique used is a delayed coincidence between two gamma rays: one coming from the positron (prompt signal) and the other coming from the neutron capture in the innermost part of the antineutrino detector (AD), containing gadolinium-doped liquid scintillator. The light created is collected by the photo-multipliers (PMTs) located in the outermost mineral oil-region. The antineutrino energy Eν is reconstructed from the positron prompt energy E prompt following the relation: Eν = E prompt +Ē n + 0.78 MeV, whereĒ n is the average neutron recoil energy. The expected number of IBD events at the d-th detector, T d , can be estimated summing up the contributions of all reactors to the detector: where N p is the number of protons in the target volume, P r th is the reactor thermal power, d denotes the efficiency of the detector and E k is the energy release per fission for a given isotope k taken from Ref. [72]. The neutrino survival probability P ee depends also on the distance from r-th reactor to d-th detector, L rd . For the antineutrino flux prediction Φ k (E) we use the parameterization given in Ref. [73] as well as the new normalization for reactor antineutrino fluxes updated in Ref. [74]. The inverse beta decay cross section σ IBD (E ν ) is taken from Ref. [75].

Daya Bay Experiment
Daya Bay is a reactor neutrino experiment with several antineutrino detectors (ADs), arranged in three experimental halls (EHs). Electron antineutrinos are generated in six reactor cores, distributed in pairs, with equal thermal power (P r th =2.9 GW th ) and detected in the EHs. The effective baselines are 512 m and 561 m for the near halls EH1 and EH2 and 1579 m for the far hall EH3 [2]. With this near-far technology Daya Bay has minimized the systematic errors coming from the ADs and thus provided until now the most precise determination of the reactor mixing angle. In the last Neutrino conference, Daya Bay has reported 621 days of data taking combining their results for two different experimental setups [3]: one with six ADs as it was published in Ref. [2] and the other after the installation of two more detectors, eight ADs in total. This new combined data set has four times more statistics in comparison with the previous Daya Bay results. Thus, the precision in the determination of the reactor mixing angle has been improved, and it is now of the order of 6%.
In this work we will consider the most recent data release by the Daya Bay Collaboration described above and we will concentrate on the total observed rates at each detector, that will be analyzed using the following χ 2 expression:

(3.2)
Here T d corresponds to the theoretical prediction in Eq. (3.1), M d is the measured number of events at the d-th AD with its backgrounds (B d ) subtracted and ω d r is the fractional contribution of the r-th reactor to the d-th AD number of events, determined by the baselines L rd and the total thermal power of each reactor. The pull parameters, used to include the systematical errors in the analysis, are given by the set (α r , ξ d , β d ) representing the reactor, detector and background uncertainties with the corresponding set of errors (σ r , σ d , σ B ). Uncertainties in the reactor related quantities are included in σ r (0.8%) while the uncorrelated combined uncertainties in the ADs are included in σ d (0.2%). σ B is the quadratic sum of the background uncertainties taken from Ref. [3]. Finally, we also consider an absolute normalization factor a norm to account for the uncertainty in the total normalization of events at the ADs, given by σ a , and coming for instance from uncertainties in the normalization of reactor antineutrino fluxes. In our analysis we will follow two different approaches concerning this parameter. In Sec. 4 we will take it equal to zero, assuming perfect knowledge of the events normalization. This hypothesis will be relaxed in Sec. 5, where we will allow for a non-zero normalization factor in the statistical analysis, being determined from the fit to the Daya Bay data. As we will see, the results obtained

Bounds on NSI from Daya Bay without normalization error
In this section we will present the bounds on the NSI couplings we have obtained using current Daya Bay reactor data. In all the results, we have assumed maximal 2-3 mixing and we have marginalized over atmospheric splitting with a prior of 3%. For definiteness we will start considering only the couplings relative to electron neutrino: (|ε e |,φ e ), for what we will switch all the other NSI parameters to zero. Next we will do the same for (|ε µ |,φ µ ) and (|ε τ |,φ τ ), that are equivalent for maximal value of θ 23 . Finally, we will consider the possibility of having all NSI couplings with the same value: ε e = ε µ = ε τ = ε. In all cases we will discuss the bounds arising from Daya Bay data in comparison with existing bounds. We will also consider the robustness of the θ 13 measurement by Daya Bay in the presence of NSI.

Constraints on electron-NSI couplings
According to the expression in Eq. (2.20), the effective survival probability in the case when only NSI with electron antineutrinos are considered is independent of the standard CP phase δ. Therefore, in our analysis we consider only two cases, one with the only relevant phase φ e set to zero, and a second case where we allow this phase to vary freely. Our results are presented in Fig. 5. From the left panel in this figure, we can confirm the behaviors shown by the iso-probability curves in the Sec. 2.4, namely, the presence of a non-zero ε e coupling has to be compensated with a slightly larger value of the reactor mixing angle θ 13 . In this panel one also sees how current Daya Bay data constrain very strongly the magnitude of the NSI coupling |ε e |, improving the current bound in Eq. (1.3) by one order of magnitude: However, the situation changes dramatically when the phase φ e is allowed to vary freely, as shown in the right panel of Fig. 5. In this case, the strong bound on |ε e | disappears since the |ε e |-linear term in the neutrino survival probability (see Eq.(2.20)) cancels out for φ e = 90 • and −90 • . Therefore, the sensitivity to |ε e | disappears and no bound can be obtained from reactor data. Concerning the determination of the reactor mixing angle, the presence of the NSIcoupling with electron antineutrinos results in the following allowed range for θ 13 : 0.020 ≤ sin 2 θ 13 ≤ 0.024 (90% C.L.) . (

4.2)
The same interval is obtained for the two panels at Fig. 5 and it also coincides exactly with the allowed range in absence of NSI. In consequence, we can say that the reactor angle determination by Daya Bay is robust in this specific case.

Constraints on muon/tau-NSI couplings
In this subsection we present the results obtained considering only the NSI parameters associated with muon and tau neutrinos. As we have discussed in Sec. 2.2, in this case, the phases δ and φ µ,τ do not appear separately in the expression of the survival probability, see Eq. (2.22). Therefore, it is enough to consider in our calculations the effective phase (δ − φ µ,τ ).
The results corresponding to this particular case are shown in Fig. 6. Here again we can see how the regions presented in the left panel of the figure agree with the behavior shown in the iso-probability plots (right panel of Fig. 3) where there is an anticorrelation between the reactor angle and the NSI coupling |ε µ,τ |. Thus, an increase in |ε µ,τ | is compensated by a shift of the preferred value of the reactor mixing angle toward smaller values, differently to what happens with the NSI coupling |ε e | in Fig. 5. The allowed interval for θ 13 in this case is given by: 0.013 ≤ sin 2 θ 13 ≤ 0.024 (90% C.L.) , while the obtained bound for the NSI coupling is the following In this case reactor data can not improve the present constraints on the NSI couplings at Eq. (1.3), and we get a limit of the same order of magnitude of the ones derived at Ref. [46]. However, in both cases the limits have been derived using different data and assumptions, and therefore, they can be regarded as complementary bounds coming from different data sets. In the right panel of Fig. 6 we show the results obtained when the phase (δ − φ µ,τ ) is allowed to vary. In this case, a wider range in the reactor mixing angle is allowed: 0.013 ≤ sin 2 θ 13 ≤ 0.036 (90% C.L.) . (4.5) The reason is that, in addition to the anticorrelation shown in the left panel, a correlation between the reactor angle and |ε µ,τ | is also possible when the cosine function in Eq. (2.23) is negative. Note, however, that even though there is a wider allowed region in the reactor angle, the bound on the |ε µ,τ | NSI coupling is nearly the same as the one obtained when the phase (δ − φ µ,τ ) is set to zero, namely: |ε µ,τ | ≤ 5.2 × 10 −2 (90% C.L.) . (4.6)

Constraints for the flavor-universal NSI case
Here we present the results obtained under the hypothesis of flavor-universal NSI, that is, we assume all NSI couplings are present and they take the same value. Therefore, we consider only two NSI parameters: |ε| and φ, in addition to the standard model parameters entering in the calculations. In this case, the effective survival probability is given by the expression in Eq. (2.24), with separate dependence on the phases δ and φ. Therefore we have considered four different cases in our analysis: one with all the phases set to zero, two cases varying only one of the phases with the other set to zero and a last case varying the two phases simultaneously. Our results are presented at Fig. 7. The left panel in Fig. 7 shows the tight constraint obtained for the magnitude of the flavor-universal NSI coupling |ε| when the phases are set to zero: |ε| ≤ 1.2 × 10 −3 (90% C.L.) . This result follows directly from the tendency already observed for |ε e | at Fig. 5. This happens because, even though the NSI couplings |ε µ,τ | are also present in the flavor-universal case, the dominant contribution comes from the non-oscillatory term depending on |ε e |.
The same behavior is also present in the middle panel of Fig. 7, where the Dirac phase δ is allowed to freely vary. In this case, the effect of the Dirac phase is just a scaling in the last term in Eq. (2.25), driven by cos δ, but again the total effect on θ 13 is dominated by the non-oscillatory term in Eq. (2.25). Under this assumption, we obtain the following bound on the magnitude of the flavor-universal coupling: Note that in the two former cases the allowed range for the mixing angle θ 13 is very close to the one obtained in the standard case, and therefore the Daya Bay determination is barely affected by the presence of NSI. However, when δ = 0 and the NSI phase φ is allowed to take different values, a completely different behavior results, as it is shown in the right panel of Fig. 7. In this case, in the fit of Daya Bay data, the phase φ takes a value such that the non-oscillatory term in Eq. (2.24) is cancelled up to order |ε| 4 . Most importantly, with a preferred value of cos φ −1.5|ε|, new terms of second order in |ε| appear at Eq. (2.24) and therefore the first order expression in Eq. (2.25) can not satisfactory explain the degeneracy between θ 13 and |ε|. Actually, the shift in the effective reactor angle is given now by this other expression (up to order |ε| 2 ): Then, it is possible to see that, even with s 13 = 0, one can reproduce the measured value of θ 13 in Daya Bay with |ε| ∼ 0.11, as shown in the corresponding plot. As a consequence of the degeneracy in the plane s 2 13 -|ε|, this specific case is not very restrictive and we get a loose bound on |ε|: |ε| ≤ 1.1 × 10 −1 (90% C.L.) . (4.10) phases sin 2 θ 13 |ε| electron-type NSI coupling φ e = 0 0.020 ≤ sin 2 θ 13 ≤ 0.024 |ε e | ≤ 0.0012 φ e free 0.020 ≤ sin 2 θ 13 ≤ 0.024 |ε e | unbound muon or tau-type NSI couplings (δ − φ µ,τ ) = 0 0.013 ≤ sin 2 θ 13 ≤ 0.024 |ε µ,τ | ≤ 0.051 (δ − φ µ,τ ) free 0.013 ≤ sin 2 θ 13 ≤ 0.036 |ε µ,τ | ≤ 0.052 universal NSI couplings δ = φ = 0 0.020 ≤ sin 2 θ 13 ≤ 0.024 |ε| ≤ 0.0012 δ free, φ = 0 0.020 ≤ sin 2 θ 13 ≤ 0.025 |ε| ≤ 0.0013 δ = 0, φ free sin 2 θ 13 ≤ 0.024 |ε| ≤ 0.110  As commented above, the presence of flavor universal NSI implies that the reactor mixing angle may be compatible with zero. Nevertheless, the degeneracy between the mixing angle and the new physics parameter |ε| observed here may be lifted by a combined analysis with accelerator long-baseline neutrino experiments. A global analysis of neutrino data assuming the simultaneous presence of NSI in reactor and accelerator neutrino data would be very useful for this purpose. Besides solving the degeneracy, the combined analysis might provide further constraints on the NSI couplings as well as improve the agreement between the preferred θ 13 value from reactors and long-baseline experiments, as discussed in Ref. [76]. However, since the production, detection and propagation of neutrinos is quite different in both kind of experiments, a global analysis would require a very detailed study with many new physics parameters involved, besides the consideration of a specific model for NSI. In any case, this point is out of the scope of the present analysis and it will be considered elsewhere.
Finally, let us comment that we have also considered the case of flavor-universal NSI with all the phases different from zero. However, we have not presented the results obtained for this setup because, in this case, the confusion between θ 13 and |ε| is complete, and therefore no information on any of the parameters can be extracted from the analysis of Daya Bay data. All the results obtained in this section are summarized in Table 2. Needless to mention, to obtain all the limits on sin 2 θ 13 presented along this section as well as in the table, we have marginalized over the NSI couplings over a wide range. Similarly, to place bounds on the NSI parameters, sin 2 θ 13 has also been allowed to float over a wide range.

Bounds on NSI from Daya Bay with 5% normalization error
In the previous section, we have not considered any normalization error in the statistical analysis of Daya Bay reactor data. This means that we have assumed a perfect knowledge of the event normalization at the experiment, disregarding the presence of uncertainties in the flux reactor normalization or in the detection cross section, among others. This procedure has been followed in most of the previous phenomenological analyses of Daya Bay data in presence of NSI [44,45]. However, a more detailed analysis of reactor data can not ignore the presence of such normalization errors. Therefore, in this section we present a χ 2 analysis of Daya Bay data using the expression defined at Eq. (3.2), where a free normalization factor is considered in order to account for the uncertainties in the total event number normalization. This point is very relevant in the study of NSI with reactor experiments, since the uncertainty in the event normalization presents a degeneracy with the zero-distance effect due to NSI. In consequence, the far over near technique exploited by Daya Bay in order to reduce the dependence upon total normalization does not work equally fine in the presence of NSI, where the non-oscillatory zero-distance effect, simultaneously present at near and far detectors, does not totally cancel. Actually, in the standard model case without NSI, the number of events expected at the near detector is given by: while, in the presence of NSI, the event number at the near detector is calculated as follows: Here a norm controls the normalization of far and near detector events in the fit and, together with the NSI couplings ε, it fixes the total zero-distance effect. As a result, if we set a norm to zero, we artificially increase the power of Daya Bay data to constrain the zero-distance effect due to NSI, getting non-realistic strong bounds on the NSI couplings. On the other hand, we can not leave the factor a norm totally free in our statistical analysis, as it is usually done in the standard Daya Bay analysis, where the factor is kept small thanks to the far over near technique. Actually, we have found that leaving the normalization factor totally free, and due to the degeneracy with the NSI couplings, it could achieve very large values, of the order of 10-20%. For this reason it is necessary the use of a prior on this magnitude. Recent reevaluations of the reactor antineutrino flux indicate an uncertainty on the total flux of about 3% [73,77]. However, an independent analysis in Ref. [78] claims that this uncertainty may have been underestimated due to the treatment of forbidden transitions in the antineutrino flux evaluation, and proposes a total uncertainty of 4%. Since the total normalization errors may also include uncertainties coming from other sources, we follow the conservative approach of taking a total uncertainty on the reactor event normalization of 5%. This is the value we have assumed for σ a in Eq. (3.2), so that a norm is constrained to be smaller than 5%.
To illustrate the differences with respect to the results obtained in the previous section, assuming no uncertainties in the event rate normalization, here we have considered only Case sin 2 θ 13 |ε| φ e = 0 0.020 ≤ sin 2 θ 13 ≤ 0.025 |ε e | ≤ 0.015 (δ − φ µ,τ ) = 0 sin 2 θ 13 ≤ 0.024 |ε µ,τ | ≤ 0.176 δ = φ = 0 0.017 ≤ sin 2 θ 13 ≤ 0.024 |ε| ≤ 0.015 the cases where all phases are set to zero. The results obtained with these assumptions are presented in Fig. 8 and Table 3. In the left panel of Fig. 8 we present the allowed region in the plane sin 2 θ 13 -|ε e | when only NSI with electron antineutrinos are present. In this case, the range for the reactor mixing angle is rather similar to the one shown in the left panel of Fig. 5: while the bound on |ε e |, however, is much weaker than the one given at Eq. (4.1) (although still slightly better than the one at Eq. (1.3)): The same bound is also obtained in the flavor-universal case for the NSI parameter |ε|, see the right panel of Fig. 8. In this case, the allowed range for the reactor mixing angle is a bit enlarged with respect to the previous one: 0.017 ≤ sin 2 θ 13 ≤ 0.024 (90% C.L.) , (5.5) due to the presence of NSI oscillation terms driven by the new physics couplings with muon and tau antineutrinos. As commented above, the loss of sensitivity to the NSI couplings is due to the degeneracy between the normalization uncertainty and the zero-distance terms induced by the presence of NSI. In this way, a larger value of the NSI parameters can be compensated with a non-zero normalization factor a norm , without spoiling the good agreement with experimental reactor data. Finally, the middle panel of Fig. 8 shows the results obtained when only NSI with muon or tau antineutrinos are considered. In this case, the cancellation between the normalization term and the zero-distance effect due to terms of second order in |ε µ,τ | results in an extended region in the sin 2 θ 13 -|ε µ,τ | plane with an upper bound of: sin 2 θ 13 ≤ 0.024 (90% C.L.) . (5.6) The bound on the NSI coupling is given by:

Comparing NSI constraints from 217 and 621 days of Daya Bay run
The new high-precision data from Daya Bay with 621 days of running time [3] has improved the measurement of sin 2 2θ 13 dramatically with a relative 1σ precision of ∼ 6% as compared to ∼ 11% obtained using the previously published 217 days of data [2]. This clearly shows the impact of the four time more statistics that the Daya Bay experiment has accumulated with the help of eight ADs in comparison with the previously released Daya Bay data set with six ADs. Now it would be quite interesting to see how much we can further constrain the allowed ranges for these NSI parameters under consideration using the new 621 days of Daya Bay data sin 2 θ 13 |ε| electron-type NSI parameters Current (621 days) 0.020 ≤ sin 2 θ 13 ≤ 0.024 |ε e | ≤ 0.0012 Previous (217 days) 0.019 ≤ sin 2 θ 13 ≤ 0.027 |ε e | ≤ 0.0024 muon or tau-type NSI parameters Current (621 days) 0.013 ≤ sin 2 θ 13 ≤ 0.024 |ε µ,τ | ≤ 0.051 Previous (217 days) 0.011 ≤ sin 2 θ 13 ≤ 0.026 |ε µ,τ | ≤ 0.070 universal NSI parameters Current (621 days) 0.020 ≤ sin 2 θ 13 ≤ 0.024 |ε| ≤ 0.0012 Previous (217 days) 0.019 ≤ sin 2 θ 13 ≤ 0.026 |ε| ≤ 0.0024 Table 4: 90% C.L. (1 d.o.f) bounds on sin 2 θ 13 and the NSI parameters using the old (new) 217 (621) days of Daya Bay data. Here all the phases are considered to be zero. We also do not consider the uncertainty on the normalization of reactor events and set a norm = 0 in the statistical analysis.
Daya Bay data in comparison with the old 217 days of Daya Bay run. In Fig. 9, we compare the performance of the current and the previous data sets of Daya Bay in constraining the allowed regions in sin 2 θ 13 -|ε e | plane (left panel) and in sin 2 θ 13 -|ε µ,τ | plane (right panel). In both the panels, the solid (dashed) lines portray the results with the new (old) 621 (217) days of Daya Bay data. For the sake of illustration, we have only chosen the cases of NSI parameters which are associated withν e (left panel) andν µ (right panel) assuming all the phases to be zero. In particular, we have focused on the situations which are presented in sections 4.1 and 4.2, where we do not consider the normalization uncertainty on the reactor events and set a norm = 0 in the statistical analysis. As compared to the old data set, with the current data, the improvement in constraining the allowed parameter space between sin 2 θ 13 and the NSI parameters is quite significant as can be readily seen from Fig. 9.
In Table 4, we present the 90% C.L. (1 d.o.f) constraints on sin 2 θ 13 and the NSI parameters using the previous (current) 217 (621) days of Daya Bay data. Here all the phases are considered to be zero. We do not consider the normalization uncertainty on the reactor events and take a norm = 0 in the statistical analysis. We can see from Table 4 that in case of electron-type and universal NSI parameters, the bounds on |ε e | and |ε| are the same while analyzing 621 days of Daya Bay data. This feature is also there in the case of 217 days of Daya Bay run. The constraints on |ε e | and |ε| get improved by factor of two when we consider the current 621 days of Daya Bay data compared to its previous 217 days of data. We also get better bounds on |ε µ,τ | using the current Daya Bay data as compared to the old data set. These results suggest that the future data from the Daya Bay experiment with more statistics is going to play an important role to further constrain the allowed parameter space for the NSI parameters. There are also marginal improvements on the bounds to sin 2 θ 13 with the current 621 days of Daya Bay data.

Summary and conclusions
The success of the currently running Daya Bay, RENO, and Double Chooz reactor antineutrino experiments in measuring the smallest lepton mixing angle θ 13 with impressive accuracy signifies an important advancement in the field of modern neutrino physics with nonzero mass and three-flavor mixing. With this remarkable discovery, the neutrino oscillation physics has entered into a high-precision era opening up the possibility of observing sub-dominant effects due to possible new physics beyond the Standard Model of particle physics. At present, undoubtedly the Daya Bay experiment in China is playing a leading and an important role in this direction. The recent high-precision and unprecedentedly copious data from the Daya Bay experiment has provided us an opportunity to probe the existence of the non-standard interaction effects which might crop up at the production point or at the detection stage of the reactor antineutrinos.
In this paper for the first time, we have reported the new constraints on the flavor nonuniversal and also flavor universal NSI parameters obtained using the currently released 621 days of Daya Bay data. While placing the bounds on these NSI parameters, we have assumed that the new physics effects are just inverse of each other in the production and detection processes of the reactor antineutrino experiment ,i.e., ε s eγ = ε d * γe . Considering this special case, we have discussed in detail the impact of the NSI parameters on the effective antineutrino survival probability expressions which we ultimately use to analyze the Daya Bay data. With this special choice of the NSI parameters, we have observed a global shift of the oscillation amplitude without distorting the shape of the oscillation probability. This shift in the depth of the oscillation dip can be caused due to the NSI parameters and as well as θ 13 , making it quite difficult to disentangle the NSI effects from the standard oscillations. Before presenting the final results, we have studied the correlations between the NSI parameters and θ 13 with the help of iso-probability surface plots in Sec. 2. This study has been quite useful to understand the final bounds on the NSI parameters that we have obtained from the fit. Since the shape of the oscillation probability is not distorted with the special choice of the NSI parameters considered in this paper, an analysis based on the total event rate at the Daya Bay experiment is sufficient to obtain the limits on the NSI parameters.
As far as the flavor non-universal NSI parameters are concerned, first we have considered the NSI parameters |ε e | and φ e which are associated withν e . Assuming φ e = 0 • and a perfect knowledge of the normalization of the event rates, the current Daya Bay data places a strong constrain on |ε e | ≤ 1.2 × 10 −3 (90% C.L.) improving the present bound on |ε e | by one order of magnitude. Now if we consider the uncertainty in the normalization of event rates with a prior of 5%, then this limit changes to 15 × 10 −3 , diluting the constraint by almost one order of magnitude. No limits can be placed on |ε e | if we allow φ e to vary freely in the fit. We have also observed that the determination of the 1-3 mixing angle is quite robust in this specific case and it is almost independent of the issue of uncertainty in the normalization of event rates and the choice of φ e . In fact, the allowed range of sin 2 θ 13 coincides exactly with the allowed range in absence of NSI. Next we turn our attention to the NSI parameters |ε µ | and φ µ which are associated withν µ . With (δ − φ µ ) = 0 • and a perfect knowledge of the event rates normalization, the current Daya Bay data sets a limit of |ε µ | ≤ 5.1 × 10 −2 (90% C.L.) which is comparable and complementary to the existing bound obtained using a different data set under different assumptions. This limit on |ε µ | becomes 17.6 × 10 −2 once we consider the uncertainty in the normalization of event rates with a prior of 5%. These limits on |ε µ | remain almost unchanged even if we allow (δ − φ µ ) to vary freely in the fit. This is not true while placing the constraint on sin 2 θ 13 . For an example, in the case when we do not consider any uncertainty in the normalization of event rates, the upper bound on sin 2 θ 13 increases from 2.4 × 10 −2 to 3.6 × 10 −2 due to freely varying (δ − φ µ ) in the fit instead of setting it to zero. We cannot place a lower bound on sin 2 θ 13 when we consider the uncertainty in the normalization of event rates with a prior of 5% even if (δ − φ µ ) is taken to be zero in the fit. The above mentioned limits on |ε µ | and sin 2 θ 13 are also valid for the NSI parameters |ε τ | and φ τ as long as the 2-3 mixing angle is maximal ,i.e., sin 2 θ 23 = 0.5.
In the case of flavor universal NSI parameters, we have placed limits on |ε| and sin 2 θ 13 under certain assumptions on δ and φ. With δ = φ = 0 • and perfect knowledge of the normalization of the event rates, the upper bound on |ε| at 90% C.L. (1 d.o.f) is 1.2 × 10 −3 . Though this limit does not change much when we marginalize over δ in the fit, it deteriorates by almost ninety times when we allow φ to vary freely in the fit providing a limit of |ε| ≤ 110 × 10 −3 (90% C.L.). With a 5% uncertainty in the normalization of total event rates and δ = φ = 0 • , the constraint on |ε| becomes ≤ 15 × 10 −3 at 90% C.L. With the same assumptions, we can restrict θ 13 within a range of 0.017 ≤ sin 2 θ 13 ≤ 0.024 (90% C.L.).
One of the interesting studies that we have performed in this paper is the comparison of the constraints on the NSI parameters placed with the current 621 days of Daya Bay data with the limits obtained using the previously released 217 days of Daya Bay run. In this analysis, all the phases are considered to be zero and the normalization of events is also kept fixed in the statistical analysis with a norm = 0. We have observed that the constraints on |ε e | and |ε| get improved by factor of two when we analyze the current 621 days of Daya Bay data compared to its previous 217 days data. This comparative study reveals the merit of the huge statistics that Daya Bay has already accumulated. It also suggests that the future high-precision data from the Daya Bay experiment with enhanced statistics is inevitable to further probe the sub-leading effects in neutrino flavor conversion due to the presence of the possible NSI parameters beyond the standard three-flavor oscillation paradigm.
A Effective Survival Probability with the NSI parameters: ε s eγ = ε d * γe In this appendix, we present the effective probability expressions for the physical scenarios where ε s eγ = ε d * γe . In such cases, the spectral study of the reactor data plays an important role since the NSI parameters are not only responsible for a shift in θ 13 i.e. the change of the depth of the first oscillation maximum but they also modify the shape of the probability due to the shift in its energy. We have already mentioned that a detailed analysis of the Daya Bay data considering such interesting physical cases will be performed in [65].