Neutrino oscillations and Non-Standard Interactions

Current neutrino experiments measure the neutrino mixing parameters with an unprecedented accuracy. The upcoming generation of experiments will be sensitive to subdominant effects that can give information on the unknown neutrino parameters: the Dirac CP-violating phase, the mass ordering and the $\theta_{23}$ octant. Determining the exact values of neutrino mass and mixing parameters is crucial to test neutrino models and flavor symmetries. In the first part of this review, we summarize the current status of neutrino oscillation parameters. We consider the most recent data from solar experiments and the atmospheric data from Super-Kamiokande, IceCube and ANTARES. We implement the data from the reactor experiments KamLAND, Daya Bay, RENO and Double Chooz as well as the long baseline data from MINOS, T2K and NOvA. If in addition to the standard interactions, neutrinos have subdominant Non-Standard Interactions (NSI) with matter, extracting the values of these parameters will suffer from new degeneracies. We review such effects and formulate the conditions on the NSI parameters under which the precision measurement of neutrino oscillation parameters can be distorted. Like standard weak interactions, NSI can be categorized into Charged and Neutral Current NSI. Our focus will be on NC NSI since it is possible to build a class of models giving rise to sizeable NC NSI with effects on neutrino oscillations. These models are based on new $U(1)$ gauge symmetry with a boson of mass $\lesssim 10$ MeV. The UV complete model should be electroweak invariant which implies that along with neutrinos, charged fermions acquire new interactions on which there are strong bounds. We enumerate the bounds that exist on such models and show that it is possible to build viable models avoiding all the bounds. We review methods to test these models and suggest approaches to break the degeneracies caused by NSI.


I. INTRODUCTION
In the framework of "old" electroweak theory, formulated by Glashow, Weinberg and Salam, lepton flavor is conserved and neutrinos are massless. As a result, a neutrino of flavor α (α ∈ {e, µ, τ }) created in charged current weak interactions in association with a charged lepton of flavor α will maintain its flavor. Various observations have however shown that the flavor of neutrinos change upon propagating long distances. Historically, solar neutrino anomaly (deficit of the solar neutrino flux relative to standard solar model predictions) [1] and atmospheric neutrino anomaly (deviation of the ratio of muon neutrino flux to the electron neutrino flux from 2 for atmospheric neutrinos that cross the Earth before reaching the detector) [2] were two main observations that showed the lepton flavor was violated in nature. This conclusion was further confirmed by observation of flavor violation of man-made neutrinos after propagating sizable distances in various reactor [3][4][5] and long baseline experiments [6,7]. The established paradigm for flavor violation which impressively explain all these anomalies is the three neutrino mass and mixing scheme. According to this scheme, each neutrino flavor is a mixture of different mass eigenstates. As neutrinos propagate, each component mass eigenstate acquires a different phase so neutrino of definite flavor will convert to a mixture of different flavors; hence, lepton flavor violation takes place.
Within this scheme, the probability of conversion of ν α to ν β (as well as that ofν α toν β ) in vacuum or in matter with constant density has a oscillatory dependence on time or equivalently on the distance traveled by neutrinos 1 . For this reason, the phenomenon of flavor conversion in neutrino sector is generally known as neutrino oscillation. Neutrino flavor eigenstates are usually denoted by ν α . That is ν α is defined as a state which appears in W boson vertex along with charged lepton l α (α ∈ {e, µ, τ }). The latter corresponds to charged lepton mass eigenstates. Neutrino mass eigenstates are denoted by ν i with mass m i where i ∈ {1, 2, 3}. The flavor eigenstates are related to mass eigenstate by a 3 × 3 unitary matrix, U , known as PMNS after Pontecorvo-Maki-Nakagawa-Sakata: ν α = i U αi ν i . The unitary mixing matrix can be decomposed as follows where the mixing angles θ ij are defined to be in the [0, π/2] range while the phase δ can vary in [0, 2π). In this way, the whole physical parameter space is covered. Historically, ν 1 , ν 2 and ν 3 have been defined according to their contribution to ν e . In other words, they are ordered such that |U e1 | > |U e2 | > |U e3 | so ν 1 (ν 3 ) provides largest (smallest) contribution to ν e . Notice that with this definition, θ 12 , θ 13 ≤ π/4. It is then of course a meaningful question to ask which ν i is the lightest and which one is the heaviest; or equivalently, what is the sign of ∆m 2 ij = m 2 i − m 2 j . The answer to this question comes from observation. Time evolution of ultra-relativistic neutrino state is governed by the following Hamiltonian: H vac + H m , where the effective Hamiltonian in vacuum is given by Within the Standard Model (SM) of particles, the effective Hamiltonian in matter H m in the framework of the medium in which neutrinos are propagating can be written as where it is assumed that the medium is electrically neutral (N e = N p ), unpolarized and composed of non-relativistic particles . In vacuum, H m = 0 and we can write By adding or subtracting a matrix proportional to the identity I 3×3 to the Hamiltonian, neutrinos obtain an overall phase with no observable physical consequences. That is why neutrino oscillation probabilities (both in vacuum and in matter) are sensitive only to ∆m 2 ij rather than to m 2 i . As a result, it is impossible to derive the mass of the lightest neutrino from oscillation data alone. Similarly, neutrino oscillation pattern within the SM only depends on N e not on N n . Similar arguments can be repeated for antineutrinos by replacing U with U * (or equivalently δ → −δ) and replacing H m → −H m . The phase δ, similarly to its counterpart in the CKM mixing matrix of quark sector, violates CP. Just like in the quark sector, CP violation in neutrino sector is given by the Jarlskog invariant: J = sin θ 13 cos 2 θ 13 sin θ 12 cos θ 12 cos θ 23 sin θ 23 sin δ.
As we will see in detail in section II, the mixing angles θ 12 , θ 13 and θ 23 are derived from observations with remarkable precision. The mixing angle θ 23 has turned out to be close 45 • but it is not clear within present uncertainties whether θ 23 < π/4 or θ 23 > π/4. This uncertainty is known as the octant degeneracy. The value of δ is also unknown for the time being, although experimental data start indicating a preferred value close to 3π/2. The absolute value and the sign of ∆m 2 21 are however determined. While |∆m 2 31 | is measured, the sign of ∆m 2 31 is not yet determined. If ∆m 2 31 > 0 (∆m 2 31 < 0), the scheme is called normal (inverted) ordering or normal (inverted) mass spectrum. The main goals of current and upcoming neutrino oscillation experiments are determining sgn(cos 2θ 23 ), sgn(∆m 2 31 ) and the value of the CP-violating phase δ.
The neutrino oscillation program is entering a precision era, where the known parameters are being measured with an ever increasing accuracy. Next generation of long-baseline neutrino experiments will resolve the subdominant effects in oscillation data sensitive to the yet unknown oscillation parameters (e.g. δ). Of course, all these derivations are within 3 × 3 neutrino mass and mixing scheme under the assumption that neutrinos interact with matter only through the SM weak interactions (plus gravity which is too weak to be relevant). Allowing for Non-Standard Interaction (NSI) can change the whole picture. Non-Standard Interaction of neutrinos can be divided into two groups: Neutral Current (NC) NSI and Charged Current (CC) NSI. While the CC NSI of neutrinos with the matter fields (e, u, d) affects in general the production and detection of neutrinos, the NC NSI may also affect the neutrino propagation in matter. As a result, both types of interaction may show up at various neutrino experiments. In recent years, the effects of both types of NSI on neutrino experiments have been extensively studied in the literature, formulating the lower limit on the values of couplings in order to have a resolvable impact on the oscillation pattern in upcoming experiments. On the other hand, non-standard interaction of neutrinos can crucially affect the interpretation of the experimental data in terms of the relevant neutrino mass parameters. Indeed, as it will be discussed in this work, the presence of NSI in the neutrino propagation may give rise, among other effects, to a degeneracy in the measurement of the solar mixing angle θ 12 [9][10][11]. Likewise, CC NSI at the production and detection of reactor antineutrinos can affect the very precise measurement of the mixing angle θ 13 in Daya Bay [12,13]. Moreover, it has been shown that NSI can cause degeneracies in deriving the CP-violating phase δ [14][15][16][17], as well as the correct octant of the atmospheric mixing angle θ 23 [18] at current and future long-baseline neutrino experiments. Along this review, we will discuss possible ways to resolve the parameter degeneracies due to NSI, by exploiting the capabilities of some of the planned experiments such as the intermediate baseline reactor neutrino experiments JUNO and RENO50 [19].
Most of the analyses involving NSI in neutrino experiments parameterize such interactions in terms of effective four-Fermi couplings. However, one may ask whether it is possible to build a viable renormalizable electroweak symmetric UV complete model that underlay this effective interaction with coupling large enough to be discernible at neutrino oscillation experiments. Generally speaking if the effective coupling comes from integrating out a new state (X) of mass m X and of coupling g X , we expect the strength of the effective four-Fermi interaction to be given by g 2 X /m 2 X . We should then justify why X has not been so far directly produced at labs. As far as NC NSI is concerned, two solutions exist: (i) X is too heavy; i.e., m X m EW . Recent bounds from the LHC imply m X > 4 − 5 TeV [20] which for even g X ∼ 1 implies g 2 X /m 2 X G F [21][22][23]. (ii) Second approach is to take m X m EW and g X 1 such that g 2 X /m 2 X ∼ G F . In this approach, the null result for direct production of X is justified with its very small coupling. In [24][25][26], this approach has been evoked to build viable models for NC NSI with large effective couplings. For CC NSI, the intermediate state, being charged, cannot be light. That is, although its Yukawa couplings to neutrinos and matter fields can be set to arbitrarily small values, the gauge coupling to the photon is set by its charge so the production at LEP and other experiments cannot be avoided. We are not aware of any viable model that can lead to a sizable CC NSI. Interested reader may see Refs. [13,[27][28][29][30][31][32][33].
This review is organized as follows. In section II, we review the different neutrino oscillation experiments and discuss how neutrino oscillation parameters within the standard three neutrino scheme can be derived. We then discuss the prospect of measuring yet unknown parameters: δ, sgn(∆m 2 31 ) and sgn(cos 2θ 23 ). In section III, we discuss how NSI can affect this picture and review the bounds that present neutrino data sets on the effective parameters. We then discuss the potential effects of NSI on future neutrino experiments and suggest strategies to solve the degeneracies. In section IV, we introduce models that can lead to effective NSI of interest and briefly discuss their potential effects on various observables. In section V, we review methods suggested to test these models. Results will be summarized in section VI.

II. NEUTRINO OSCILLATIONS
In this section, we will present the current status of neutrino oscillation data in the standard three-neutrino framework. Most recent global neutrino fits to neutrino oscillations can be found in Refs. [34][35][36]. Here we will focus on the results of Ref. [34], commenting also on the comparison with the other two analysis. First, we will describe the different experiments entering in the global neutrino analysis, grouped in the solar, reactor, atmospheric and longbaseline sectors. For each of them we will also discuss their main contribution to the determination of the oscillation parameters.
A. The solar neutrino sector: (sin 2 θ12, ∆m 2 21 ) Under the denomination of solar neutrino sector, one finds traditionally not only all the solar neutrino experiments, but also the reactor KamLAND experiment, sensitive to the same oscillation channel, under the assumption of CPT conservation. Solar neutrino analysis include the historical radiochemical experiments Homestake [37], Gallex/GNO [38], SAGE [39], sensitive only to the interaction rate of electron neutrinos, but not to their energy or arrival time to the detector. This more detailed information was available with the start-up of the real-time solar neutrino experiment Kamiokande [40], that confirmed the solar neutrino deficit already observed by the previous experiments. Its successor Super-Kamiokande, with a volume 10 times larger, has provided very precise observations in almost 20 years of operation. Super-Kamiokande is a Cherenkov detector that uses 50 kton of ultra pure water as target for solar neutrino interactions, that are detected through elastic neutrino-electron scattering. This process is sensitive to all neutrino flavors, with a larger cross section for ν e due to the extra contribution from the chargedcurrent neutrino-electron interaction. The correlation between the incident neutrino and the recoil electron in the observed elastic scattering makes possible the reconstruction of the incoming neutrino energy and arrival direction. After its first three solar phases [41][42][43], Super-Kamiokande is already in its fourth phase, where a very low energy detection threshold of 3.5 MeV has been achieved [44]. Moreover, during this last period, Super-Kamiokande has reported a 3σ indication of Earth matter effects in the solar neutrino flux, with the following measured value of the day-night asymmetry [45,46] Likewise, they have reconstructed a neutrino survival probability consistent with the MSW prediction at 1σ [47,48]. The Sudbury Neutrino Observatory (SNO) used a similar detection technique with 1 kton of pure heavy water as neutrino target. The use of the heavy water allowed the neutrino detection through three different processes: charged-current ν e interactions with the deuterons in the heavy water (CC), neutral-current ν α with the deuterons (NC), and as well as elastic scattering of all neutrino flavors with electrons (ES) The measurement of the neutrino rate for each of the three reactions allows the determination of the ν e flux and the total active ν α flux of 8 B neutrinos from the Sun. SNO took data during three phases, each of them characterized by a different way of detecting the neutrinos produced in the neutrino NC interaction with the heavy water [49,50].
Apart from Super-Kamiokande, the only solar neutrino detector at work nowadays is the Borexino experiment. Borexino is a liquid-scintillator experiment sensitive to solar neutrinos through the elastic neutrino-electron scattering, with a design optimized to measure the lower energy part of the spectrum. During its first detection phase, Borexino has reported precise observations of the 7 Be solar neutrino flux, as well as the first direct observation of the mono-energetic pep solar neutrinos and the strongest upper bound on the CN O component of the solar neutrino flux [51]. Moreover, Borexino has also measured the solar 8 B rate with a very low energy threshold of 3 MeV [52] and it has also provided the first real-time observation of the very low energy pp neutrinos [53].
The simulation of the production and propagation of solar neutrinos requires the knowledge of the neutrino fluxes produced in the Sun's interior. This information is provided by the Standard Solar Model (SSM), originally built by John Bahcall [54]. The more recent versions of the SSM offer at least two different versions according to the solar metallicity assumed [55,56]. Ref. [34] uses the low metallicity model while [36] reports its main results for the high-metallicity model. For a discussion on the impact of the choice of a particular SSM over the neutrino oscillation analysis see for instance Refs. [36,57].
KamLAND is a reactor neutrino experiment designed to probe the existence of neutrino oscillations in the so-called LMA region, with ∆m 2 21 ∼ 10 −5 eV 2 . KamLAND detected reactor antineutrinos produced at an average distance of 180 kilometers, providing the first evidence for the disappearance of neutrinos traveling to a detector from a power FIG. 1. Allowed regions at 90% and 99% C.L. from the analysis of solar data (black lines), KamLAND (blue lines) and the global fit (colored regions). θ13 has been marginalized according to the latest reactor measurements [34].
reactor [58]. In KamLAND, neutrinos are observed through the inverse beta decay processν e + p → e + + n, with a delayed coincidence between the positron annihilation and the neutron capture in the medium that allows the efficient reduction of the background. The final data sample released by KamLAND contains a total live time of 2135 days, with a total of 2106 reactor antineutrino events observed to be compared with 2879±118 reactor antineutrino events plus 325.9±26.1 background events expected in absence of neutrino oscillations [59]. Fig. 1 reports the allowed region in the sin 2 θ 12 -∆m 2 21 plane from the analysis of all solar neutrino data (black lines), from the analysis of the KamLAND reactor experiment (blue lines) and from the combined analysis of solar + KamLAND data (colored regions). Here the value of the θ 13 has been marginalized following the most recent short-baseline reactor experiments which will be described in the next subsection. From the figure, one can see that the determination of θ 12 is mostly due to solar neutrino experiments, while the very accurate measurement of ∆m 2 21 is obtained thanks to the spectral information from KamLAND. There is also a mild but noticeable tension between the preferred values of ∆m 2 21 by KamLAND and by solar experiments. While the first one shows a preference for ∆m 2 21 = 4.96 × 10 −5 eV 2 , the combination of all solar experiments prefer a lower value: ∆m 2 21 = 7.6 × 10 −5 eV 2 . This discrepancy appears at the 2σ level. As we will see in the next section, non-standard neutrino interactions have been proposed as a way to solve the tension between solar and KamLAND data.
The best fit point for the global analysis corresponds to: Maximal mixing is excluded at more than 7σ.

B. Short-baseline reactor neutrino experiments and θ13
Until recently, the mixing angle θ 13 was pretty much unknown. Indeed, the only available information about the reactor angle was an upper-bound obtained from the non-observation of antineutrino disappearance at the CHOOZ and Palo Verde reactor experiments [60,61], sin 2 θ 13 < 0.039 at 90% C.L. for ∆m 2 31 = 2.5 × 10 −3 eV 2 . Later on, the interplay between different data samples in the global neutrino oscillation analyses started showing some sensitivity to the reactor mixing angle θ 13 . In particular, from the combined analysis of solar and KamLAND neutrino data, a non-zero θ 13 value was preferred [57,62,63]. The non-trivial constraint on θ 13 mainly appeared as a result of the different correlation between sin 2 θ 12 and sin 2 θ 13 present in the solar and KamLAND neutrino data samples [64,65]. Moreover, a value of θ 13 different from zero helped to reconcile the tension between the ∆m 2 21 best fit points for solar and KamLAND separately. Another evidence for a non-zero value of θ 13 was obtained from the combination of atmospheric and long-baseline neutrino data [66,67]. Due to a small tension between the preferred values of |∆m 3 31 | at θ 13 = 0 for MINOS Super-Kamiokande atmospheric data, the combined analysis of both experiments showed a preference for θ 13 > 0 [68][69][70].
Nevertheless, the precise determination of θ 13 was possible thanks to the new generation of reactor neutrino experiments, Daya Bay, RENO and Double Chooz. The main features of these new reactor experiments are, on one side, their increased reactor power compared to their predecessors and, on the other side, the use of several identical antineutrino detectors located at different distances from the reactor cores. Combining these two features results in an impressive increase on the number of detected events. Moreover, the observed event rate at the closest detectors is used to predict the expected number of events at the more distant detectors, without relying on the theoretical predictions of the antineutrino flux at reactors. Several years ago, in the period between 2011 and 2012, the three experiments found evidence for the disappearance of reactor antineutrinos over short distances, providing the first measurement of the angle θ 13 [3][4][5]. We will now briefly discuss the main details of each experiment as well as their latest results.
The Daya Bay reactor experiment [4] in China is a multi-detector and multi-core reactor experiment. Electron antineutrinos produced at six reactor cores with 2.9 GW thermal power are observed at eight antineutrino 20 ton Gadolinium-doped liquid scintillator detectors, located at distances between 350 and 2000 m from the cores. The latest data release from Daya Bay has reported the detection of more than 2.5 millions of reactor antineutrino events, after 1230 days of data taking [71]. This enormous sample of antineutrino events, together with a significant reduction of systematical errors has made possible the most precise determination of the reactor mixing angle to date [71] sin 2 2θ 13 = 0.0841 ± 0.0027 (stat.) ± 0.0019 (syst.) .
Likewise, the sensitivity to the effective mass splitting ∆m 2 ee has been substantially improved, 2 reaching the accuracy of the long-baseline accelerator experiments, originally designed to measure this parameter.
The RENO experiment [5] in South Korea consists of six aligned reactor cores, equally distributed over a distance of 1.3 km. Reactor antineutrinos are observed by two identical 16 ton Gadolinium-doped Liquid Scintillator detectors, located at approximately 300 (near) and 1400 m (far detector) from the reactor array center. The RENO Collaboration has recently reported their 500 live days observation of the reactor neutrino spectrum [73,74], showing an improved sensitivity to the atmospheric mass splitting, |∆m 2 ee | = 2.62 +0.24 −0.26 ×10 −3 eV 2 . Their determination for θ 13 is consistent with the results of Daya Bay: sin 2 2θ 13 = 0.082 ± 0.009 (stat.) ± 0.006 (syst.) (9) The Double Chooz experiment in France detects antineutrinos produced at two reactor cores in a near and far detectors located at distances of 0.4 km and 1 km from the neutrino source, respectively [3]. The latest results presented by the Double Chooz collaboration correspond to a period of 818 days of data at far detector plus 258 days of observations with the near detector. From the spectral analysis of the multi-detector neutrino data, the following best fit value for θ 13 is obtained [75] sin 2 2θ 13 = 0.119 ± 0.016 (stat. + syst.).
(10) Fig. 2 illustrates the sensitivity to sin 2 θ 13 obtained from the analysis of reactor and global neutrino data for normal and inverted mass ordering. The black line corresponds to the result obtained from the combination of all the reactor neutrino data samples while the others correspond to the individual reactor data samples, as indicated. One can see from the figure that the more constraining results come from Daya Bay and RENO experiments, while Double Chooz shows a more limited sensitivity to θ 13 . Moreover, the global constraint on θ 13 is totally dominated by the Daya Bay measurements, with some contributions from RENO to its lower bound. Notice also that global analyses of neutrino data do not show relevant differences between the preferred value of θ 13 for normal or inverted mass ordering, as we will discuss later. For more details on the analysis of reactor data presented in Fig. 2, see Ref. [34]. The atmospheric neutrino flux was originally studied as the main source of background for the nucleon-decay experiments [76][77][78]. For several years, most of the dedicated experiments observed a deficit in the detected number of atmospheric neutrinos with respect to the predictions. The solution to this puzzling situation arrived in 1998, when the observation of the zenith angle dependence of the µ-like atmospheric neutrino data in Super-Kamiokande indicated an evidence for neutrino oscillations [2]. Some years later, the Super-Kamiokande Collaboration reported a L/E distribution of the atmospheric ν µ data sample characteristic of neutrino oscillations [79]. Super-Kamiokande has been taking data almost continuously since 1996, being now in its fourth phase. Super-Kamiokande is sensitive to the atmospheric neutrino flux in the range from 100 MeV to TeV. The observed neutrino events are classified in three types, fully contained, partially contained and upward-going muons, based on the topology of the event. The subsequent data releases by the Super-Kamiokande Collaboration have increased in complexity. Currently it is very complicated to analyze the latest results by independent groups [34,36]. From the analysis of the latest Super-Kamiokande atmospheric data, the following best fit values have been obtained for the oscillation parameters [80]: Thus, a slight preference for θ 23 > π/4 is reported. Likewise, the normal mass ordering (i.e., ∆m 2 31 > 0) is preferred over the inverted one (i.e., ∆m 2 31 < 0) .
In recent years, atmospheric neutrinos are also being detected by neutrino telescope experiments. IceCube and ANTARES, originally designed to detect higher energy neutrino fluxes, have reduced their energy threshold in such a way that they can measure the most energetic part of the atmospheric neutrino flux.
The ANTARES telescope [81], located under the Mediterranean Sea, observes atmospheric neutrinos with energies as low as 20 GeV. Neutrinos are detected via the Cherenkov light emitted after the neutrino interaction with the medium in the vicinity of the detector. In Ref. [82], the ANTARES Collaboration has analyzed the atmospheric neutrino data collected during a period of 863 days. Their result for the oscillation parameters are in good agreement with current world data. For the first time, ANTARES results have also been included in a global neutrino oscillation fit [34].
The IceCube DeepCore detector is a sub-array of the IceCube neutrino observatory, operating at the South Pole [83]. DeepCore was designed with a denser instrumentation compared to IceCube, with the goal of lowering the energy threshold for the detection of atmospheric muon neutrino events down to 10 GeV. Neutrinos are identified trough the Cherenkov radiation emitted by the secondary particles produced after their interaction in the ice. The most recent data published by DeepCore correspond to a live time of three years [84]. A total of 5174 atmospheric neutrino events were observed, compared to a total 6830 events expected in absence of neutrino oscillations. The obtained best fit values for the atmospheric neutrino parameters sin 2 θ 23 = 0.53 +0.09  2 31 parameter. Both plots correspond to the normal ordered neutrino mass spectrum [34].
also compatible with the atmospheric results of the Super-Kamiokande experiment.
The left panel of Fig. 3 shows the allowed regions at 90% and 99% C.L. in the atmospheric neutrino oscillation parameters sin 2 θ 23 and ∆m 2 31 obtained from ANTARES, DeepCore and Super-Kamiokande phases I to III [34]. From the combination one sees that DeepCore results start being competitive with the determination of the atmospheric oscillation parameters by long-baseline experiments. Indeed, a recent reanalysis of DeepCore atmospheric data [85] shows an improved sensitivity with respect to the region plotted in Fig. 3. The sensitivity of ANTARES shown in Fig. 3 is not yet competitive with the other experiments. However, it is expected that the ANTARES collaboration will publish an updated analysis that will certainly improve their sensitivity to the atmospheric neutrino parameters.

D. Long-baseline accelerator experiments
After the discovery of neutrino oscillations in the atmospheric neutrino flux, several long-baseline accelerator experiments were planned to confirm the oscillation phenomenon with a man-made neutrino source. The first two experiments trying to probe the ν µ disappearance oscillation channel in the same region of ∆m 2 as explored by atmospheric neutrinos were K2K and MINOS. Their successors, T2K and NOνA are still at work today.
The KEK to Kamioka (K2K) experiment used a neutrino beam produced by a 12 GeV proton beam from the KEK proton synchrotron. The neutrino beam was detected by a near detector 300 m away from the proton target and by the Super-Kamiokande detector, at a distance of 250 km. The number of detected neutrino events, as well as the spectral distortion of the neutrino flux observed by K2K was fully consistent with the hypothesis of neutrino oscillation [86].
The Main Injector Neutrino Oscillation Search (MINOS) experiment observed neutrino oscillations from a beam produced by the NuMI (Neutrinos at Main Injector) beamline at Fermilab in an underground detector located at the Soudan Mine, in Minnesota, 735 km away. MINOS searched for oscillations in the disappearance (ν µ → ν µ ) and appearance channels (ν µ → ν e ), for neutrinos and antineutrinos as well. After a period of 9 years, the MINOS experiment collected a data sample corresponding to an exposure of 10.71 × 10 20 protons on target (POT) in the neutrino mode, and 3.36 × 10 20 POT in the antineutrino mode [87,88]. The combined analysis of all MINOS data shows a slight preference for inverted mass ordering and θ 23 below maximal as well as a disfavored status for maximal mixing with ∆χ 2 = 1.54 [89]. The allowed ranges for the atmospheric parameters from the joint analysis of all MINOS data are the following The Tokai to Kamioka (T2K) experiment uses a neutrino beam consisting mainly of muon neutrinos, produced at the J-PARC accelerator facility and observed at a distance of 295 km and an off-axis angle of 2.5 • by the Super-Kamiokande detector. The most recent results of the T2K collaboration for the neutrino and antineutrino channel have been published in Refs. [90,91]. A separate analysis of the disappearance data in the neutrino and antineutrino channels has provided the determination of the best fit oscillation parameters for neutrinos and antineutrinos [90]. The obtained results are consistent so, no hint for CPT violation in the neutrino sector has been obtained [92,93]. In both cases, the preferred value of the atmospheric angle is compatible with maximal mixing. The combined analysis of the neutrino and antineutrino appearance and disappearance searches in T2K, that corresponds to a total sample of 7.482 × 10 20 POT in the neutrino mode, and 7.471 × 10 20 POT in the antineutrino mode, results in the best determination of the atmospheric oscillation parameters to date [91] sin 2 θ 23 = 0.532 (0.534) , |∆m 2 32 | = 2.545 (2.510) × 10 −3 eV 2 , for normal (inverted) mass ordering spectrum. Furthermore, thanks to the combination of neutrino and antineutrino data, T2K has already achieved a mild sensitivity to the CP violating phase, reducing the allowed 90% C.L. range of δ in radians to [−3.13, −0.39] for normal and [−2.09, −0.74] for inverted mass ordering [91].
In the NOνA (NuMI Off-Axis ν e Appearance) experiment, neutrinos produced at the Fermilab's NuMI beam are detected in Ash River, Minnesota, after traveling 810 km through the Earth. In the same way as the T2K experiment, the NOνA far detector is located slightly off the centerline of the neutrino beam coming from Fermilab. Thanks to this configuration, a large neutrino flux is obtained at energies close to 2 GeV, where the maximum of the muon to electron neutrino oscillations is expected. The most recent data release from the NOνA collaboration corresponds to an accumulated statistics of 6.05×10 20 POT in the neutrino run [94,95]. For the muon antineutrino disappearance channel, 78 events have been observed, to be compared with 82 events expected for oscillation and 473 ± 30 events predicted under the no-oscillation hypothesis. The searches for ν µ → ν e transitions in the accelerator neutrino flux have reported the observation of 33 electron neutrino events, with an expected background of 8.2±0.8 ν e events. The analysis of the NOνA Collaboration disfavors maximal values of θ 23 at the 2.6σ level [94]. On the other hand, from the analysis of the appearance channel it is found that the inverted mass ordering is disfavored at 0.46σ, due to the small number of event predicted for this ordering in comparison to the observed results [95]. Furthermore, the combination of appearance and disappearance NOνA data with the θ 13 measurement at the reactor experiments results disfavors the scenario with inverted neutrino mass ordering and θ 23 < π/2 at 93% C.L., regardless of the value of δ [95].
The right panel of Fig. 3 shows the 90% and 99% C.L. allowed region in the atmospheric neutrino oscillation parameters sin 2 θ 23 and ∆m 2 31 according to the MINOS, T2K and NOνA data for normal mass ordering [34]. Note the different scale for ∆m 2 31 with respect to the left panel. The three long-baseline experiments provide similar constraints for this parameter, while the constraint on θ 23 obtained from T2K is a bit stronger. One can also see some small differences between the preferred values of θ 23 by the three experiments. While T2K prefers maximal mixing, MINOS and NOνA show a slight preference for non-maximal θ 23 . In any case, these differences are not significant and the agreement among the three experiments is quite good. Although not shown here, the agreement for inverted mass ordering is a bit worse, since in that case the rejection of NOνA against maximal mixing is stronger, whereas the preference for θ 23 ∼ π/4 in T2K remains the same as for normal ordering.

E. Global fit to neutrino oscillations
In the previous subsections, we have reviewed the different experimental neutrino data samples, discussing their dominant sensitivity to one or two oscillation parameters. However, every data sample offers subleading sensitivities to other parameters as well. Although the information they can provide about such parameters may be limited, in combination with the rest of data samples, relevant information can emerge. This constitutes the main philosophy behind global analyses of neutrino oscillation data: joint analyses trying to exploit the complementarity of the different experiments to improve our knowledge on the neutrino oscillation parameters. Here, we will show the results of a combined analysis of neutrino oscillation data in the framework of the three-flavor neutrino oscillation scheme. Fig. 4 reports the 90%, 95% and 99% C.L. allowed regions in the parameters sin 2 θ 23 , sin 2 θ 13 , |∆m 2 31 | and δ from the global fit in Ref. [34] for normal and inverted mass ordering. For the allowed regions in the solar plane sin 2 θ 12 -∆m 2 21 see Fig. 1. The best fit points, along with the corresponding 1σ uncertainties and 90% C.L. ranges for each parameter, are quoted in Table. I. The relative uncertainties on the oscillation parameters at 1σ range from . Allowed regions at 90%, 95% and 99% C.L. in the planes sin 2 θ23-|∆m 2 31 |, sin 2 θ13-|∆m 2 31 |, sin 2 θ23-δ and sin 2 θ13-δ for normal (lines) and inverted mass ordering (colored regions). The star indicates the global best fit point, corresponding to normal ordering, while the circle indicates the local minimum in inverted ordering. Adapted from Ref. [34].
parameter best fit ± 1σ 90% C.L. range  around 2% for the mass splittings to 7-10% (depending on the mass ordering) for sin 2 θ 23 . In case of the CP phase, the 1σ uncertainties are of the order of 15-20%. Note also that, at the 3σ level, the full range of δ is still allowed for normal ordering. For the case of inverted ordering, a third of the total range is now excluded at the 3σ level. These results are in good agreement with Refs. [35,36].
Despite the remarkable sensitivity reached in the determination of most of the neutrino oscillation parameters, there are still three unknown parameters in the oscillation of standard three neutrino scheme: the octant of θ 23 , the value of the CP phase δ and the neutrino mass ordering. The current status of these still unknown parameters will be discussed next.
Let us now comment on the maximality/non-maximality and octant preference for the atmospheric mixing angle. So far, experimental neutrino data have not shown a conclusive preference for values of θ 23 smaller, equal or larger than π/4. Different experiments may show a limited preference for one of the choices, but for the moment all the results are consistent at the 3σ level. On the other hand, one finds that the available global analyses of neutrino data [34][35][36], using very similar data samples show slightly different results for the octant preference. For this particular case, one can find the origin of the possible discrepancies in the different treatment of the Super-Kamiokande atmospheric data. See the previous references for more details on the chosen approach at each work. The results in Fig. 4 and Table I, corresponding to the analysis in Ref. [34], show a preference for θ 23 in the first octant. This global best fit point corresponds to normal mass ordering, but a local minimum can also be found with θ 23 > π/4 and inverted mass ordering with a ∆χ 2 = 4.3. In the same way, additional local minima can be found with θ 23 in the second octant and inverted mass spectrum and the other way around. All these possibilities are allowed at 90% C.L. as can be seen in the right panels of Fig. 4. With current data, the status of the maximal atmospheric mixing is a bit delicate, being allowed only at 99% C.L. However, this result may change after the implementation of the partially published data release of T2K [96] in the global fit.
In the same way, the current neutrino oscillation data do not offer a definitive determination for the neutrino mass ordering. Individual neutrino experiments show in general a limited sensitivity to the mass ordering, with the exception of the latest atmospheric data from Super-Kamiokande, that prefer normal mass ordering with a significance of ∆χ 2 = 4.3. Note however that this data sample is not included in some of the global analyses of neutrino oscillations [34,36]. The sensitivity to the mass ordering in the global analysis arises instead from the interplay of the different neutrino data, as a result of the existing correlations and tensions among the other neutrino parameters. Indeed, the three global analysis discussed in this review show a preference for normal mass ordering, although the significance may be different in each case, depending on the particular details of the specific global fit. In the work in Ref. [34], discussed in a bit more details here, a preference for normal ordering over inverted is obtained, with a significance of ∆χ 2 = 4.3. In any case, the results reported are not conclusive yet, and we will have to wait for the next generation of experiments devoted to this purpose (among others), such as DUNE [97], PINGU [98], ORCA [99], JUNO [100] or RENO-50 [101].
Finally, we comment on the sensitivity to the CP-violating phase δ. Prior to the publication of the antineutrino run data from T2K, combined analyses were already showing a weak preference for δ = 3π/2, while δ = π/2 was disfavored above the 2σ level [102][103][104]. This sensitivity, absent in all the individual data samples, emerged from the tension between the value of θ 13 measured at the reactor experiments and the preferred value of θ 13 for δ = π/2 in T2K. This scenario has changed after the release of T2K results from its antineutrino run and now the sensitivity to δ comes mainly from the combined analysis of the neutrino and antineutrino channel in T2K. The remaining experiments contribute only marginally to the determination of the CP-violating phase.

III. CURRENT BOUNDS ON NON-STANDARD INTERACTIONS
New neutrino interactions beyond the Standard Model are natural features in most neutrino mass models [105,106]. As commented in the introduction, these Non-Standard Interactions (NSI) may be of Charged-Current (CC) or of Neutral-Current (NC) type. In the low energy regime, neutrino NSI with matter fields can be formulated in terms of the effective four-fermion Lagrangian terms as follows: where G F is the Fermi constant and P X denote the left and right chirality projection operators P R,L = (1 ± γ 5 )/2.
The dimensionless coefficients f f X αβ and f X αβ quantify the strength of the NSI between leptons of α and β flavour and the matter field f ∈ {e, u, d} (for NC-NSI) and f = f ∈ {u, d} (for CC-NSI). At the limit f X αβ → 0, we recover the standard interactions, while αβ ∼ 1 corresponds to new interactions with strength comparable to that of SM weak interactions. If αβ is non-zero for α = β, the NSIs violate lepton flavor. If αα − ββ = 0, the lepton flavor universality is violated by NSI.
The presence of neutrino NSI may affect the neutrino production and detection at experiments as well as their propagation in a medium through modified matter effects [47,48]. In the literature, it is common denoting the CC-NSI couplings as s αβ or d αβ since they may often affect the source (s) and detector (d) interactions at neutrino experiments. On the other hand, mf αβ is used to refer to the NC-NSI couplings with the fundamental fermion f generally affecting the neutrino propagation in matter (m). In this case, what is relevant for neutrino propagation in a medium is the vector part of interaction f V αβ = f L αβ + f R αβ . In fact, the neutrino propagation in a medium is sensitive  [11]. See this reference for further details.
to the following combinations 3 so most of the bounds from oscillation experiments are presented in the literature in terms of αβ rather than in terms of f V αβ . Inside the Sun, N u /N e 2N d /N e 1 [55] and inside the Earth, N u /N e N d /N e 3 [107]. When studying the effect of NSI at neutrino detection, there will be independent sensitivity for the left and right chirality coefficients f L αβ and f R αβ . Although this kind of interactions has not been confirmed experimentally, their potential effects have been extensively studied in a large variety of physical scenarios. As a result, stringent bounds on their strength have been derived [11,105,106]. Moreover, it has been shown that NSI may interfere with neutrino oscillations in different contexts, giving rise to parameter degeneracies that can affect the robustness of the neutrino parameter determination. In this section, we will review these results.

A. NSI in solar experiments
NSI may affect the propagation of solar neutrinos within the Sun and the Earth as well as the detection, depending on the type of NSI considered. Before the confirmation of neutrino oscillation as the phenomenon responsible for the solar neutrino anomaly by the KamLAND experiment, NSI with massless neutrinos was also proposed as the mechanism behind this anomaly [108][109][110][111][112]. After KamLAND confirmed the phenomenon of mass-induced electron neutrino (antineutrino) oscillation, NSI was excluded as the main mechanism behind the solar neutrino oscillations, although its presence has been considered at subleading level in solar neutrino experiments, see for instance [9-11, 113, 114]. These analyses have found that a small amount of NSI, dV ee 0.3, is in better agreement with data than the standard solution at the level of 2σ. On one hand, this result is due to the non-observation of the up-turn of the solar neutrino spectrum predicted by the standard LMA-MSW solution at around 3 MeV. On the other hand, there exists a small tension between the preferred value of ∆m 2 21 by KamLAND and by solar experiments that can be eased by introducing NSI. More surprisingly, these studies revealed an alternative solution to the standard LMA-MSW, known as LMA-Dark or LMA-D solution [9][10][11], requiring NSI with strength dV τ τ − dV ee 1. The presence of this new degenerate solution to the solar neutrino anomaly, shown in Fig. 5, can be understood in the framework of two-neutrino mixing as follows. Under this approximation (justified by the fact that sin 2 θ 13 1), the two by two Hamiltonian matrix can be diagonalized with an effective mixing angle given by 3 Note that some references use the It is very relevant to distinguish between both notations, since the reported bounds will be different by a factor of 3. For this reason, we prefer to quote directly the results in terms of the effective lagrangian coefficients f V αβ . In any case, for the analysis including Earth matter effects, the relation between the corresponding NSI couplings is straightforward.
The splitting between two eigenvalues is given by Under the simultaneous transformations cos 2θ 12 → − cos 2θ 12 and (H m ) ee → −(H m ) ee , we find θ m 12 → π/2−θ m 12 and ∆ m → ∆ m which means that the off-diagonal elements of the 2 × 2 Hamiltonian remains the same but the diagonal elements (the 11 and 22 elements) flip. That is, P (ν e → ν e ) changes to P (ν c → ν c ) where ν c ≡ c 23 ν µ − s 23 ν τ is the combination that ν e converts to (that is ν c |ν 3 = ν c |ν e = 0). Since in two neutrino approximation, we can write P (ν e → ν e ) + P (ν e → ν c ) = 1, P (ν c → ν e ) + P (ν c → ν c ) = 1 and P (ν e → ν c ) = P (ν c → ν e ). We therefore conclude P (ν e → ν e ) = P (ν c → ν c ). As a result, under the transformation described above, P (ν e → ν e ) remains invariant. This transformation is not possible for the case of standard matter effects, where the value of (H m ) ee is fixed to √ 2G F N e . However, if one considers the presence of neutrino NSI with the matter field f , the effective Hamiltonian in the medium is modified to: Allowing for sufficiently large values of the f V ee coupling in the effective Hamiltonian in matter H m , it is now possible to apply the transformation described above, obtaining a degenerate solution to the solar neutrino anomaly with cos 2θ 12 < 0. Notice that, letting cos 2θ 12 to change sign, we are violating the historical choice of keeping θ 12 in the first octant and, therefore, ν 1 will not be anymore the state giving the largest contribution to ν e , but the lighter one between those two eigenstates that give the main contribution to ν e . This change of definition is in fact equivalent to maintain the same convention regarding the allowed range for the mixing angle, but allowing ∆m 2 21 to be negative. Indeed, changing ∆m 2 21 → −∆m 2 21 instead of cos 2θ 12 → − cos 2θ 12 , we would find the same degeneracy. In other words, for a given H m , solar neutrino data only determine the sign of the product ∆m 2 21 cos 2θ 12 , not the signs of ∆m 2 21 and cos 2θ 12 separately, and therefore there is a freedom in definition. Since the LMA-D solution was introduced in the literature keeping ∆m 2 21 positive while allowing θ 12 to vary in the range (0, π/2) [9] and this convention has become popular in the literature since then, we will use it along this review. Note also that the degeneracy found at the neutrino oscillation probability is exact only for a given composition of matter (i.e., for a given N n /N p = N n /N e ). The composition slightly varies across the Sun radius and of course is quite different for the Sun and the Earth. Because of this, the allowed regions in the neutrino oscillation parameter space for the LMA and LMA-D solutions are not completely degenerate. A small χ 2 difference between the best fit point of the LMA solution and the local minimum of LMA-Dark solution appears because the relevant data analyses take into account the varying composition of the Sun and the day-night asymmetry due to propagation in the Earth.
Unfortunately, the degeneracy between the LMA and LMA-Dark solutions could not be lifted by the KamLAND reactor experiment because KamLAND was not sensitive to the octant of the solar mixing angle due to the lack of matter effects. A possible way to solve this problem was proposed in Ref. [10]. There, it was found that the combination of solar experiments, KamLAND and neutrino neutral-current scattering experiments, such as CHARM [115], may help to probe the LMA-D solution. The relevance of the degeneracy in the solar neutrino parameter determination has been explored recently in Ref. [116]. As discussed in this analysis, the ambiguity of LMA-D does not affect only the octant of the solar mixing angle but it also makes impossible the determination of the neutrino mass ordering at oscillation experiments. More recently, a global analysis of neutrino scattering and solar neutrino experiments was performed to further investigate the situation of the LMA-D solution [117]. Besides the accelerator experiment CHARM, the authors also considered the NuTeV experiment [118]. They found that the degenerate LMA-D solution may be lifted for NSI with down quarks, although it does not disappear for the case of neutrino NSI with up quarks. As discussed in that work, constraints from CHARM and NuTeV experiments can be however directly applied only for NSI with relatively heavy mediators. For the case of NSI mediated by lighter particles (above 10 MeV), constraints coming from coherent neutrino-nucleus scattering experiments may be used to resolve the degeneracy. Indeed, after the recent observation of such process at the COHERENT experiment [119], a combined analysis of neutrino oscillation data including the observed number of events in this experiment has excluded the LMA-D solution (for up and down quarks) at the 3σ level [120]. One should, however, bear in mind that the analysis [120] assumes the mediator of interaction in Eq. (16) is heavier than 50 MeV. As we shall see in sect V, for light mediator, their conclusion should be revised. Besides that, COHERENT data along with neutrino oscillation data has been used to improve the current bounds on the flavor-diagonal NSI parameters [120] These limits on NC vector interactions of ν τ improve previous bounds by one order of magnitude [10,11,121]. For the flavor-changing NC NSI couplings, however, the improvement is much smaller: The spectrum of coherent elastic neutrino-nucleus scattering events at COHERENT has also been analyzed to constrain the amplitude of NSI in Ref. [122]. Besides their impact on solar neutrino propagation, NSI can also affect the detection processes at solar neutrino experiments. In experiments like Super-Kamiokande and Borexino, for instance, the presence of NSI may modify the cross section of neutrino elastic scattering on electrons, used to observe solar neutrinos. Analyzing data from solar neutrino experiments, and in particular the effect of NSI on neutrino detection in Super-Kamiokande, in combination with KamLAND, Ref. [123] reported limits on the NSI parameters which are competitive and complementary to the ones obtained from laboratory experiments. For the case of ν e NSI interaction with electrons, the reported bounds (taking one parameter at a time) are: while for the case of ν τ NSI interaction with electrons, looser constraints are obtained: The sensitivity of the Borexino solar experiment to NSI has also been investigated in Refs. [124,125]. Using 7 Be neutrino data from Borexino Phase I, the following 90% C.L. bounds have been derived [125] − 0.046 < eL ee < 0.053 , −0.21 < eR ee < 0.16 , As can be seen, the NSI constraints obtained from Borexino and the combined analysis of solar (mainly Super-Kamiokande) and KamLAND data are comparable. It is expected that future results from Borexino Phase II, as well as the combination of all solar data, including Borexino, plus KamLAND data would allow a significant improvement on the current knowledge of neutrino NSI with matter.

B. NSI in atmospheric neutrino experiments
The impact of non-standard neutrino interactions on atmospheric neutrinos was originally considered in Refs. [126][127][128][129]. Assuming a two-flavor neutrino system, it was shown [127] that the presence of large NSI couplings together with the standard mechanism of neutrino oscillation can spoil the excellent description of the atmospheric neutrino anomaly given by neutrino oscillations. Thus, quite strong bounds on the magnitude of the non-standard interactions were derived. Using atmospheric neutrino data from the first and second phase of the Super-Kamiokande experiment, the following constraints were obtained, under the two-flavor neutrino approach [130]: However, Ref. [128,129] showed that a three-family analysis significantly relaxes the previous bounds in such a way that the values of the NSI couplings with quarks comparable to the standard neutral current couplings can be still compatible with the Super-Kamiokande atmospheric data. A more recent three-neutrino analysis of NSI in the atmospheric neutrino flux can be found in Ref. [131], where the following limits on the effective NSI couplings with electrons have been obtained: − 0.035 (−0.035) < eV µτ < 0.018 (0.035) , | eV τ τ − eV µµ | < 0.097 (0.11) (90% C.L.) (28) for the case of real (complex) eV µτ coupling.
The IceCube extension to lower energies, DeepCore, has made possible the observation of atmospheric neutrinos down to 5 GeV with unprecedented statistics. Indeed, with only 3 years of data, DeepCore allows the determination of neutrino oscillation parameters with similar precision as the one obtained from the long-lived Super-Kamiokande or the long-baseline accelerator experiments [85]. Using the most recent data release from DeepCore, the IceCube collaboration has reported the following constraints on the flavor-changing NSI coupling [132]: From a different data sample containing higher energy neutrino data from IceCube, the authors of Ref. [133] have derived somewhat more restrictive bounds on the same NSI interactions: Both results are fully compatible and constitute the current best limits on NSI in the ν µ -ν τ sector. Note, however, that unlike in Ref. [131], these results have been obtained under the assumption of no flavor diagonal NSI couplings, Future prospects on NSI searches in atmospheric neutrino experiments have been considered in the context of PINGU, the future project to further lower the energy threshold at the IceCube observatory. Ref. [134] shows that, after three years of data taking in PINGU in the energy range between 2 and 100 GeV, the Super-Kamiokande constraints on the NSI couplings may be improved by one order of magnitude: Likewise, the impact of NSI interactions on atmospheric neutrinos on the future India-based Neutrino Observatory (INO) has been analyzed in Ref. [135]. Besides discussing its constraining potential towards NSI, this work studies how the sensitivity to the neutrino mass hierarchy of INO, one of the main goals of the experiment, may change in the presence of NSI.

C. NSI in reactor experiments
Modern reactor neutrino experiments, like Daya Bay, RENO and Double Chooz, provide a very accurate determination of the reactor mixing angle θ 13 [73,136,137]. Being at the precision era of the neutrino parameter determination, it is imperative to investigate the robustness of this successful measurement in the presence of NSI. Refs. [12,13,[138][139][140] have addressed this point. In principle, short-baseline reactor experiments may be affected by the presence of new neutrino interactions in β and inverse-β decay processes, relevant for the production and detection of reactor antineutrinos [141]. The NSI parameters relevant for these experiments are the CC NSI couplings between up and down quarks, positrons and antineutrinos of flavor α, ud eα . Considering unitarity constraints on the CKM matrix as well as the non-observation of neutrino oscillations in the NOMAD experiment, one may find the following constraints on these CC NSI couplings [33]: Ref. [13] explored the correlations between the NSI parameters and the reactor mixing angle determination, showing that the presence of NSI may lead to relatively large deviations in the measured value of θ 13 in Daya Bay, as it can be seen in Fig. 6. Conversely, the total number of events observed in Daya Bay was used to constrain the corresponding NSI couplings under two assumptions: i) perfect theoretical knowledge of the reactor neutrino flux in absence of NSI and ii) assuming a conservative error on its total normalization. In the latter case, it was shown that assuming an uncertainty of 5% on the reactor flux can relax the bounds by one order of magnitude, obtaining the following conservative limits on the NSI strengths 5 with P=L,R,V,A. Note that these results improve the existing bounds on the ud ee coupling reported above. On the other hand, one finds that an improved knowledge of the standard absolute neutrino flux from nuclear reactors together with a larger data sample from Daya Bay will result in a more stringent bound on the other two couplings in the near future. Notice also that previous results have been obtained assuming that the NSI couplings at neutrino production and detection satisfy s αβ = d * αβ . In this case, the presence of NSI would only produce a shift in the oscillation amplitude without altering the L/E pattern of the oscillation probability, and therefore, the analysis of the total neutrino rate in Daya Bay provides enough information. The investigation of more exotic scenarios where s αβ = d * αβ will require the spectral analysis of the Daya Bay data [12].
NSI at future intermediate baseline reactor experiments like JUNO and RENO-50 (see, for instance, Refs. [142] and [140]) are discussed at Section V.

D. NSI in long-baseline neutrino experiments
Besides neutrino production and detection, NSI can also modify the neutrino propagation through the Earth in long-baseline accelerator experiments 6 . This effect will be larger for experiments with larger baselines such as MINOS or NOνA. Using their neutrino and antineutrino data sample, the MINOS Collaboration reported the following bounds on the flavor-changing NC NSI with electrons [143]: − 0.20 < eV µτ < 0.07 (90% C.L.) .
MINOS appearance data were also used to constrain NSI interactions between the first and third family [144], although the reported bound, | eV eτ | < 3.0 (90% C.L.) does not improve the previous limits on that parameter [33].
Regarding the long-baseline experiment NOνA, the presence of NSI in the neutrino propagation has been proposed as a way to solve the mild tension between the measured values of the atmospheric mixing angle in T2K and NOνA [145]. Under this hypothesis, the deviation from maximal mixing in the NOνA preferred value for θ 23 would be explained through the NSI-modified matter effects. The T2K experiment, with a shorter baseline, has a limited sensitivity to matter effects in the neutrino propagation so, its θ 23 measurement would be unaffected by NSI. Note, however, that the size of the NSI required to reproduce the observed results is of the same order as the standard neutrino interaction (to be more precise eτ , ( τ τ − µµ ) ( τ τ − ee ) ∼ O(1)).
The presence of NSI has also been considered to reconcile the measured value of θ 13 in reactor experiments and T2K [146]. In that case, it is suggested that CC-NSI in the neutrino production and detection processes may be responsible for the different values of the reactor mixing angle measured in Daya Bay and T2K. Finally, it has been shown that long-baseline neutrino facilities can also suffer from degeneracies in the reconstruction of some parameters due to the existence of new neutrino interactions with matter. For instance, Ref. [14] states that NC NSI may affect the sensitivity to the CP-violating phase δ in experiments like T2K and NOνA. According to this analysis, it would be possible confusing signals of NSI with a discovery of CP violation, even if CP is conserved in nature. This result is illustrated in Fig. 7, where it is shown how the standard CP-violating scenario may be confused with an hybrid standard plus NSI CP-conserving scenario.
Future sensitivities to NSI as well as the presence of new degeneracies due to NSI in future long-baseline experiments such as DUNE, T2HK and T2HKK are analyzed in more detail in Sect. V. It is worth mentioning that CC-NSI, affecting the production and detection of neutrinos can show up also in short baseline experiments [147][148][149].

E. NSI in non-oscillation neutrino experiments
Neutrino scattering experiments constitute a very precise tool towards the understanding of neutrino interactions with matter. Indeed, this kind of experiments has been often used to measure the electroweak mixing angle θ W [150]. Non-standard neutrino interactions may contribute significantly to the neutrino-electron elastic scattering cross section and therefore they cannot be ignored when studying this process. Ref. [151] compiled most of the neutrino scattering experiments potentially modified by the presence of NSI, from the neutrino accelerator-based experiments LSND and CHARM to the short-baseline neutrino reactor experiments Irvine, Rovno and MUNU, including as well as the measurement of the process e + e − → ννγ at LEP. From a combined analysis of all experimental data, allowed ranges on the e αβ were obtained. Some of these results are among the current strongest constraints on NSI couplings, and are reported in Table II. The antineutrino-electron scattering data collected by the TEXONO Collaboration has been also used to constrain the presence of neutrino NC NSI with electrons [152] as well as CC NSI at neutrino production and detection [153].
In order to constrain the NSI between neutrinos and quarks, one may use data from the neutrino-nucleus experiments NuTeV, CHARM and CDHS. From the combination of atmospheric and accelerator data from NuTeV, CHARM and CDHS, the following limits on the non-universal vectorial and axial NSI parameters were derived [154]: For the case of the flavor changing NSI couplings (with q = u, d) | qV µτ | < 0.007 , | qA µτ | < 0.039 (90% C.L.).
Under this category we include also the first observation of coherent neutrino-nucleus scattering observed at the COHERENT experiment recently [119]. As discussed above, the COHERENT data have been used to constrain neutrino NSI with quarks in Refs. [120,122]. The combination of solar neutrino oscillation data with COHERENT has been exploited to investigate the status of the solar degenerate solution LMA-D.  Here we summarize the current constraints on the NSI couplings from different experiments discussed throughout this section. For more details about the assumptions considered in each case, we refer the reader to the previous subsections as well as to the original references where the constraints have been calculated. The limits summarized in Tables II, III and IV have been obtained assuming only one nonzero NSI coupling at a time. Table II contains the limits on the flavor diagonal NC NSI couplings between neutrinos and electrons eP αα and neutrinos and quarks qP αα , with P = L, R, V, A being the chirality index and q = u, d. The table indicates the origin of the reported bound as well as the reference where it has been obtained as well. Most of the limits have been derived from the combination of neutrino oscillation and detection or production experimental results. For instance, the joint analysis of atmospheric neutrino data and accelerator measurements in NuTeV, CHARM and CDHS [154], or solar and KamLAND data together with the recent bounds of COHERENT [120]. 7 In other cases the constraints reported in the table come just from one type of experiment, as the limits derived only from CHARM [121], TEXONO [152] or atmospheric data [131]. Note that, for the latter case, we have adapted the bound on eV τ τ reported in Ref. [131] to the corresponding bound for quarks, qV τ τ . Table III collects the limits of the flavor changing NC NSI couplings between neutrinos and electrons eP αβ and neutrinos and quarks qP αβ , with the same conventions indicated above for P and q. As discussed before, in this case most of the bounds also emerge from the complementarity of different types of experiments, as the combination of reactor and accelerator non-oscillation experiments in Ref. [151]. On the other hand, the first analyses on NSI obtained from IceCube data [132,133] offer very strong bounds on qV µτ . This last constraint has also been adapted to get the equivalent bound for NSI with electrons, eV µτ .
Finally, Table IV contains the limits on the neutrino CC NSI with quarks and electrons (semileptonic CC NSI) and the CC NSI with leptons only (purely-leptonic CC NSI) in terms of the couplings udP αβ and ll P αβ , respectively. The former ones, have been discussed in the context of the neutrino production and detection in the Daya Bay reactor experiment, as analyzed in Ref. [13]. Previous bounds on this type of NSI have been derived using the negative searches for neutrino oscillations at short distances in the NOMAD experiment [155,156], as reported in the table [33]. Constraints on leptonic CC NSI using the results of the KARMEN experiment [157] as well as the deviations of Fermi's constant G F in the presence of these interactions, have also been obtained in Ref. [33]. We refer the reader to that work for further details on the derivation of these constraints. As we saw in the previous section, neutral current NSI of neutrinos with matter fields can lead to observable effect on neutrino oscillation provided that the NSI parameters αβ are large enough. As briefly discussed in the introduction, it is possible to build viable models for NSI by invoking an intermediate state of relatively light mass (∼ 10 MeV) which has escaped detection so far because of its very small coupling. In this chapter, we review the models that give rise to sizeable NSI through integrating out a new gauge boson Z with a mass smaller than ∼ 100 MeV. We however note that an alternative model has been suggested [158] in which NSI are obtained from SU (2) L scalar doublet-singlet mixing. We shall not cover this possibility in the present review. The models described in this chapter introduce a new U (1) gauge interaction which is responsible for NSI between neutrinos and quarks.
In section IV A, we describe the general features of the model gauging a linear combination of lepton flavors and Baryon number with a light O(10 MeV) gauge boson. We then outline general phenomenological consequences. We show how a simple economic model can be reconstructed to reproduce the NSI pattern that gives the best fit to neutrino data, solving the small tension between KamLAND and solar neutrino by explaining the suppression of the upturn in the low energy part of the solar neutrino spectrum. In section IV B, we describe another model which can provide arbitrary flavor structure u αβ = d αβ (both lepton flavor violating and lepton flavor conserving) without introducing new interactions for charged leptons. In sect. IV C, the impact of the recent results from the COHERENT experiment is outlined.

A. NSI from new U (1)
In this section, we show how we can build a model based on U (1) × SU (2) L × U (1) Y gauge symmetry which gives rise to NSI for neutrinos. Notice that the NSI of interest for neutrino oscillation involves only neutrinos and quarks of first generation which make up the matter. However, to embed the scenario within a gauge symmetric theory free from anomalies, the interaction should involve other fermions.
Let us first concentrate on quark sector and discuss the various possibilities of U (1) charge assignment. Remember that, in the flavor basis by definition, the interaction of W µ boson with quarks is diagonal: W µ On the other hand, as discussed in the introduction, in order to maintain the SM prediction for the total neutrino flux measured at the SNO experiment via NC interactions, the coupling to (at least the first generation of) quarks should be non-chiral. Thus, the U (1) charges of u 1L , u 1R , d 1L and d 1R should be all equal. In principle, different generations of quarks can have different U (1) charges. Such a freedom opens up abundant possibilities for anomaly cancelation. However, if the coupling of the new gauge bosons to different quark generations is non-universal, in the quark mass basis, off-diagonal couplings of form Z µqi γ µ q j | i =j appear which can lead to q i → Z q j with a rate enhanced by m 3 qi /m 2 Z due to longitudinal component of Z . To avoid these decays, we assume the quarks couple to Z universally. In other words, the U (1) charges of quarks are taken to be proportional to baryon number, B. Yukawa couplings of quarks to the SM Higgs will be then automatically invariant under U (1) .
Let us now discuss the couplings of leptons to the new gauge boson. There are two possibilities: 1) U (1) charges are assigned to a combination of lepton numbers of different flavors. In this case, the U (1) charges of charged leptons and neutrinos will be equal. (2) Neutrinos couple to Z through mixing with a new fermion with mass larger than m Z . In this case charged leptons do not couple to Z at tree level. We shall return to the second case in sect. IV B. In the present section, we focus on the first case. As discussed in [26], it is possible to assign U (1) charge to linear combinations of leptons which do not even correspond to charged lepton mass eigenstates. However, let us for the time being study the charge assignment as follows Denoting the new gauge coupling by g , the coupling of each generation of leptons and quarks to Z are, respectively, g a α and g /3. There are strong bounds on the new couplings of the electrons. If a e = 0, Z with a mass of ∼ 10 MeV will dominantly decay into e − e + so strong bounds from beam dump experiments apply. These bounds combined with supernova cooling study yield g a e < 3 × 10 −11 (see Fig. 4 of [159]). On the other hand, for m Z < m π , the bound from π 0 → γZ is g < 3 × 10 −3 [160]. (See Fig. 8 which is taken from [25].) These bounds are too stringent to lead to a discernible ee . We therefore set a e = 0 which means at tree level, neither electron nor ν e couple to Z . With such charge assignment, we obtain Notice that, with this technique, we only obtain lepton flavor conserving NSI. Notice that for neutrino oscillation not only the absolute value of αα − ββ but also its sign is important. In fact, neutrino oscillation data favor positive value of ee − µµ ee − τ τ ∼ 0.3. If a µ + a τ = −3, the anomalies cancel without any need for new generations of leptons and/or quarks. However, just like in B − L and L µ − L τ gauge theories, the presence of right-handed neutrinos is necessary to cancel the U (1) − U (1) − U (1) anomaly. Let us take a µ = a τ = −3/2 so that anomalies cancel; moreover, we obtain µµ = τ τ . We can then accommodate the best fit with g = 4 × 10 −5 m Z 10 MeV For the LMA-Dark solution ee − µµ < 0 is required, so the value of a µ a τ > 0. As a result, more chiral fermions are needed to be added to cancel anomalies. We will return to this point later.
Since the U (1) charges of the left-handed and right-handed charged leptons are equal, their Yukawa coupling (and therefore their mass terms) preserve U (1) automatically. We should however consider the mass matrix of neutrinos with more care. While the flavor diagonal elements of neutrino mass matrix can be produced without any need for U (1) breaking, if a α is not universal, obtaining the neutrino mass mixing requires symmetry breaking. As mentioned above, right-handed neutrinos are also required to cancel anomalies. If the masses of neutrinos are of Dirac type, right-handed neutrinos will be as light as left-handed neutrinos. They can be produced in the early universe via U (1) coupling so, if they are light, they can contribute to the relativistic degrees of freedom. To solve both problems at one shot, we can invoke the seesaw mechanism. For simplicity, we take a µ = a τ so that the mixing between the second and third generation does not break U (1) . Generalization to a µ = a τ will be straightforward. Let us call the right-handed neutrino of generation "i" with N i . Under U (1) , Dirac mass terms come from By changing the basis, either of λ 4 and λ 5 can be set to zero, but the nonzero one will mix the second and the third generations. Moreover, we add electroweak singlet scalars S 1 and S 2 with U (1) charges −2a µ and −a µ , respectively. We can then write the following potential  π→eνZ' FIG. 10. 90% C.L. constraints on i g 2 ei versus m Z from constraints on π −→ eνZ [162] and K + −→ e + ννν [163] branching ratios, from current and projected Rπ measurement by PIENU [164] and from the RK measurement by NA62 [165]. gei is the coupling of Z µνe γ µ νi where νi can be any neutrino state much lighter than ∼ 100 MeV.
Once S 1 and S 2 develop a vacuum expectation value (VEV), U (1) will be broken leading to the desired neutrino mass and mixing scheme. The VEVs of S 1 and S 2 induce a mass of for the Z boson. Taking g a µ ∼ 10 −5 , we find that as long as S 1 ∼ S 2 ∼TeV, the contribution to the Z mass will be ∼ 10 MeV as desired. In case that more scalars charged under U (1) are added to the model (we shall see examples in sect IV B), the Z mass receives further contributions. For m Z < m π , the Z can decay only to neutrino pair at tree level with a lifetime of As a result, Z evades the bounds from the beam dump experiments. In the following, we go through possible experiments that can search for the Z boson.  11. 90% C.L. constraints on i g 2 µi versus m Z from K + −→ µ + ννν branching ratio [166]. The band shows the parameter space within Lµ gauge models (giving rise to equal couplings to µ and νµ) that can explain the (g − 2)µ anomaly [167]. gµi is the coupling of Z µνµ γ µ νi where νi can be any neutrino state much lighter than 100 MeV. [168][169][170][171] which are too weak to be relevant for our models; see Figs. 10 and 11, which are taken from [168]. As shown in [168], the bound on the ν e coupling to the Z boson can be dramatically improved by customized searches for three body decays (K + → e + + missing energy) and (π + → e + + missing energy).
• The Z − γ mixing: In principle, Z can kinetically mix with the hypercharge gauge boson which gives rise to Z mixings both with the photon and the Z bosons. Even if we set the kinetic mixing to zero at tree level, it can be produced at loop level as long as there are particles charged under both U (1) gauge symmetries. Going to a basis where the kinetic terms of gauge bosons is canonical, the Z boson obtains a coupling to the electron given by e where is the kinetic mixing between Z and the photon. This coupling can affect neutrino interaction with the electron on which there are strong bounds from solar experiments (mostly Super-Kamiokande and Borexino) [123,125]. Ref [172], setting the tree level kinetic mixing equal to zero, has calculated the kinetic mixing for the L µ − L τ models and has found it to be finite and of order of eg /8π 2 . The Borexino bound [125] can then be translated into g e < ∼ 10 −4 which can be readily satisfied for g < 10 −4 . The loop contribution to the photon Z mixing from a charged particle is very similar to its contribution to the vacuum polarization (photon field renormalization) replacing (qe) 2 with (qe)(a α g ). In case of the L µ − L τ gauge symmetry, a e = 0, a µ = −a τ and since the electric charges of µ and τ are the same, the infinite parts of their contribution to the mixing cancels out. In general, we do not however expect such a cancelation and counter terms are therefore required. Once we open up the possibility of tree level kinetic mixing, the sum of tree level and loop level mixing can be set to arbitrarily small value satisfying any bound.
• Z −Z mixing: The above discussion on the Z −γ kinetic mixing also applies to the Z −Z kinetic mixing. Here, we should also check the Z − Z mass mixing [173]. It is straightforward to show that, since the Z couplings are taken to be non-chiral, there is no contribution to the Z − Z mass mixing at one loop level. If the model contains scalars that are charged both under electroweak and U (1) and develop VEV, mass mixing between Z and Z appears even at tree level. In the minimal version of the model that is described above there is no such scalar but we shall come back to this point in sect. IV B.
• Big Bang Nucleosynthesis (BBN): Decay of Z to neutrino pairs can warm up the neutrino background during and right after the BBN era. The effect can be described by the contribution to the effective extra relativistic degrees of freedom ∆N ef f . As shown in Ref. [172], BBN bounds rule out m Z < 5 MeV. Of course, this lower bound on m Z applies only if the coupling is large enough to bring Z to thermal equilibrium with neutrinos before they decouple from the plasma at T ∼ 1 MeV. That is, for 1 MeV < m Z < 5 MeV the coupling should be smaller than ∼ 3×10 −10 (m Z /5 MeV) while for m Z < 1 MeV it should be smaller than ∼ 7×10 −11 (MeV/m Z ).
• Supernova cooling: NSI can leave its imprint on the flavor composition of supernova neutrino flux [174][175][176]. Moreover, Z particle can be produced and decay back to neutrinos within the supernova core. This leads to a shortening of the mean free path of neutrinos inside the supernova core [172]. This, in turn, results in prolonging the duration of neutrino emission from supernova. To draw a quantitative conclusion and bound, a full simulation is required.
• Dips in the spectrum of high energy neutrinos: Once we introduce the new interaction for neutrinos, high energy neutrinos (or antineutrinos) traveling across the universe resonantly interact with cosmic background antineutrinos (or neutrinos) producing Z which decays back to νν pair with energies lower than that of initial neutrino (or antineutrino). This will result in a dip in the spectrum of high energy neutrinos. Taking the cosmic background neutrinos as non-relativistic, we expect the position of the dip to be given by E ν ∼ PeV(m Z /10 MeV) 2 (0.05 eV/m ν ). The value is tantalizingly close to the observed (but by no means established) gap in the high energy IceCube data. Moreover, as shown in [172] with g ∼ 10 −5 − 10 −4 (the range of interest to us), the optical depth is larger than one. Thus, this rather robust prediction can be eventually tested by looking for the dip in the high energy neutrino data.
• Magnetic dipole moment of the muon: The contribution from the Z loop to (g − 2) µ can be estimated as g 2 /(8π) up to corrections of order O(m 2 Z /m 2 µ ) ∼ 0.01. For g < 10 −4 , the contribution is too small to explain the claimed discrepancy [177].
• Neutrino scattering experiments: The amplitude of the contribution from t-channel Z exchange to neutrino quark scattering is suppressed relative to that from SM by a factor of m 2 Z /(t − m 2 Z ), where t is the Mandelstam variable. At CHARM and NuTeV experiments, the energy momentum exchange was about 10 GeV (t m 2 Z ), so the new effects were suppressed. As a result, the bound found in [10,117,121] does not apply to the model with a light gauge boson. However, as discussed in [25,117,178], low energy scattering experiments can be sensitive to low mass gauge interactions. Three categories of scattering experiments have been studied in this regard: (1) Scattering of solar neutrino at direct dark matter search experiments [25,161,178]. As shown in [25], the upcoming Xenon based experiments such as LUX-Zeplin and the future Germanium based experiments such as superCDMS at SNOLAB can test most of the parameter space of our interest (see Fig. 9, adapted from [25]). (2) As shown in detail in [117,179], the running COHERENT experiment [180] is an ideal setup to probe NSI with a light mediator. At this experiment, low energy ν µ and ν e fluxes are produced via pion and muon decay at rest. The LMA-Dark solution can be entirely probed by this experiment [117]. The COHERENT experiment has recently released its preliminary results, ruling out a significant part of the parameter space. We shall discuss the new results in sect. IV C. (3) Scattering of reactorν e flux off nuclei can also probe NSI of type we are interested in [181? ? -186].
• Neutrino trident production: The Z gauge boson coupled to ν and µ can contribute to the so-called neutrino trident production ν + A → ν + A + µ + + µ − , where A is a nucleus. The rate of such interaction was measured by the CCFR [187] and CHARM II [188] collaborations, and is found to be consistent with the SM prediction. This observation sets the bound g a < 9 × 10 −4 for m Z ∼ 10 MeV [189].
As we saw earlier, taking a µ = a τ = −3/2, the contributions from the field content of the SM to anomalies cancel out. We can then obtain any negative values of µµ − ee = τ τ − ee ∼ −1 by choosing g ∼ 10 −4 (| µµ − ee |) 1/2 (see Eq. (39)). Let us now discuss if with this mechanism we can reconstruct a model that embeds the LMA-Dark solution with positive µµ − ee = τ τ − ee ∼ 1. The condition µµ − ee = τ τ − ee = 1 can be satisfied if a e = 0, a µ = a τ > 0 and We should however notice that with a µ = a τ > 0, the cancelation of U (1) − SU (2) − SU (2) and U (1) − U (1) − U (1) anomalies require new chiral fermions. A new generation of leptons with U (1) charge equal to −(3 + a µ + a τ ) can cancel the anomalies but in order for these new fermions to acquire masses large enough to escape bounds from direct production at colliders, their Yukawa couplings enter the non-perturbative regime. Similar argument holds if we add a new generation of quarks instead of leptons. Another option is to add a pair of new generations of leptons (or quarks) with opposite U (1) Y charges but equal U (1) charge of −(3 + a µ + a τ )/2. Let us denote the field content of the fourth generation with ν R4 , e R4 and L 4 , and similarly that of the fifth generation with ν R5 , e R5 and L 5 . As pointed out, the hypercharges of fourth and fifth generation are opposite so we can write Yukawa terms of type where S is singlet of the electroweak symmetry group, SU (2) × U (1) with a U (1) charge of 3 + a µ + a τ . Even for Y 1 ∼ Y 2 ∼ 1, in order to obtain heavy enough mass, S should be of order of TeV. On the other hand, S contributes to Z mass so Masses of 4th and 5th generation S < 5 TeV m Z 10 MeV 2 × 10 −6 g (3 + a µ + a τ ) .
In other words, g In general, the cancelation of U (1) −SU (2)−SU (2) and U (1) −U (1)−U (1) anomalies requires new chiral fermions charged under U (1) and SU (2) × U (1) (or both). In the former case, we need new U (1) charged scalars whose VEV contribute to the Z mass. The lower bounds on the masses of new particles set a lower bound on the VEV of new scalars which, in turn, can be translated into an upper bound on g /m Z which leads to µµ < ∼ 0.01. In the second case, large masses of the 4th and 5th generations requires non-perturbative Yukawa coupling to the Higgs. If masses of the new fermions could be about few hundred GeV, none of these obstacles would exist. Fortunately, there is a trick to relax the strong lower bounds from colliders on the masses of new particles. Let us suppose the charged particles are just slightly heavier than their neutral counterparts. Their decay modes can be then of type e − 4(5) → ν 4(5) lν, ν 4(5) q q with a final charged lepton or jet too soft to be detected at colliders. In this case, the new generation can be as light as few 100 GeV so their mass can come from a perturbative Yukawa coupling to the SM Higgs or new scalars charged under U (1) and VEV of ∼ 100 GeV opening up a hope for g /m Z ∼ 10 −5 MeV −1 and therefore for µµ = τ τ ∼ 1.

B. A model both for LF conserving and LFV NSI
As mentioned in section IV A, the coupling of Z to neutrinos can be achieved with two mechanisms: (i) The (ν α Lα ) doublet is assigned a charge under U (1) , so neutrinos directly obtain a gauge coupling to the Z boson. This route was discussed in section IV A. (ii) Active neutrinos mix with a new fermion singlet under electroweak group symmetry but charged under new U (1) . In this section, we focus on the second route. Using the notation in Eq.(37), in the present scenario one has a e = a µ = a τ = 0 so to cancel the U (1) − SU (2) − SU (2) and U (1) − U (1) − U (1) anomalies, we should add new fermions. As discussed in the previous section, in order to make these new fermions heavier than ∼ 1 TeV, we need new scalars charged under U (1) with a VEV of 1 TeV. To keep the contribution from the new VEV to Z mass under control, Notice that this tentative bound is stronger than the bound from π → Z γ (see Fig. 8). Let us introduce a new Dirac fermion Ψ which is neutral under electroweak symmetry but charged under U (1) . Its U (1) charge denoted by a Ψ can be much larger than one. Since we take equal U (1) charges for Ψ L and Ψ R , no anomaly is induced by this Dirac fermion. Let us denote the mixing of Ψ with neutrino of flavor α with κ α . Such a mixing of course breaks U (1) . Mixing can be obtained in two ways: • We add a sterile Dirac N (neutral both under electroweak and under U (1) ) and a scalar (S) to break U (1) . The U (1) charge of the S is taken to be equal to that of Ψ. We can then add terms like the following to the Lagrangian: Notice that we were allowed to add a term of λ R SΨ L N R too, but this term is not relevant for our discussion.
Taking Y α H , λ L S , m Ψ m N , we can integrate out N and obtain Since we take λ L S < m N , in order to have sizeable κ α , the mass of Ψ cannot be much larger than Y α H . On the other hand, Y α determines the new decay mode of H → νN which is observationally constrained [190]. We therefore find an upper bound on m Ψ of few GeV. For example, taking m Ψ = 2 GeV, m N = 20 GeV, Y α H = 0.1 GeV, λ L = 1 and S ∼ 4 GeV, we obtain κ α 0.01. With such small Y α , the rate of the Higgs decay into N and ν will be as small as Γ(H → µµ) and therefore negligible. With g < 10 −4 , the contribution from S to m Z will also be negligible.
• Another scenario which has been proposed in [25] invokes a new Higgs doublet H which has U (1) charge equal to that of Ψ. The Yukawa coupling will be then equal to L = − α y αLα H T cΨ which leads to where tan β = H / H . The VEV of H can contribute to the Z mass so we obtain Thus, to obtain sizeable κ α (e.g., κ α > 0.03), we find Moreover, the VEV of H can induce Z − Z mixing on which there are strong bounds [191]. These bounds can be translated into cos β < 10 −4 (m Z /10 MeV)(1/g Ψ ) which is slightly weaker than the bound in (48). The smallness of H , despite its relatively large mass, can be explained by adding a singlet scalar(s) of charge a Ψ with L = −µS † 1 H † H which induces H = −µ S 1 /(2M 2 H ). Taking S 1 µ M 2 H , we find cos β 1. The components of H can be pair produced at colliders via electroweak interactions. They will then decay to Ψ and leptons. In particular, the charged component H − can decay into charged lepton plus Ψ which appears as missing energy. Its signature will be similar to that of a charged slepton [25]. According to the present bounds [25,192], m H 300 GeV.
Regardless of the mechanism behind the mixing between Ψ and ν, it will lead to the coupling of Z to active neutrinos as follows which leads to u αβ = d αβ = Notice that if the mixing of Ψ with more than one flavor is nonzero, we can have lepton flavor violating NSI with αβ | 2 . If more than one Ψ is added, we may label the mixing of ith Ψ to ν α with κ iα . The Schwartz inequality ( i κ * iα κ iβ ) 2 < ( i |κ iα | 2 )( i |κ iβ | 2 ) then still applies Taking κ iα = δ iα , meaning that each Ψ i mixes with only one ν α , only diagonal elements of αβ will be nonzero, preserving lepton flavor.
Under certain assumptions, Ref. [193] also derives independent bound on κ α κ * β | α =β from lepton flavor violation (LFV) processes l − α → l − β γ, but these bounds are valid only for m Ψ m W . For our case with m Ψ m W , a GIM mechanism is at work and suppresses the contribution to l − α → l − β γ from ν − Ψ mixing. In the case that the mixing comes from Yukawa coupling to H , because of the LFV induced by H and Ψ coupling to more than one flavor, a new contribution to l − α → l − β γ appears. As shown in [25], from Br(τ → eγ) < 3.3 × 10 −8 , Br(τ → µγ) < 4.4 × 10 −8 [195] and Br(µ → eγ) < 4.2 × 10 −13 [196], one can respectively derive |y e y τ | < 0.46(m H − /(400 GeV)) 2 , |y µ y τ | < 0.53(m H − /(400 GeV)) 2 and |y e y µ | < 7 × 10 −4 (m H − /(400 GeV)) 2 . As demonstrated in [25], except for eµ which is strongly constrained by the bound from µ → eγ within the model described in [25], all components of αβ can be within the reach of current and upcoming long baseline neutrino experiments. If mixing is achieved via the mechanism described in Eq. (47), no new bound from LFV rare decay applies and we can obtain all αβ (including eµ ) of the order of Notice that in this model, the coupling of Z to neutrino pairs can be much larger than the coupling to quarks: |g a Ψ κ α κ β | g , see Eq. (50). The bounds from meson decays on Z coupling to neutrino pairs have been studied in [168]. The results are shown in Figs. 10 and 11. The strongest bound for m Z ∼ 10 MeV is of order of 0.001, which can be readily satisfied if κ α κ β < 10 −3 . However, further data on [π + → e + (µ + ) + missing energy] or on [K + → e + (µ + ) + missing energy] can probe parts of the parameter space of interest to us.

C. Impact of recent results from COHERENT experiment
Recently, the COHERENT experiment has released its first results confirming the SM prediction of elastic scattering of neutrinos off nuclei at 6.7σ, studying the interaction of ν µ ,ν µ and ν e flux from Spallation Neutron Source (SNS) at the Oak Ridge National Laboratory on a 14.6 kg CsI[Na] scintillator detector [119]. The preliminary results already set strong bounds on NSI.
Assuming the validity of the contact interaction approximation (i.e., assuming the mass of the mediator is heavier than ∼ 10 MeV), Ref. [120] shows that the recent COHERENT data rules out LMA-Dark solution. Ref. [122], taking a universal coupling of Z to SM fermions finds that Let us discuss how this bound can constrain our model(s) for NSI. Regardless of the details of the underlying theory, we can write where g q = g B /3 is the coupling of Z boson to quarks and (g ν ) αβ is its coupling to ν α and ν β . In the model described in sect IV A, (g ν ) αβ = δ αβ a α g and in the model of sect IV B, (g ν ) αβ = g a Ψ κ α κ β . Remember that, in order for αβ to show up in neutrino oscillation experiments, it should be non-universal. For example, the LMA-Dark solution requires µµ − ee = τ τ − ee ∼ 1. Setting ee = 0 and µµ = 0, we expect the bound on | µµ | from the COHERENT experiment to be slightly weaker than that found in Ref. [122] taking µµ = ee = 0. Thus, the LMA-Dark solution is may still survive with the present COHERENT data, although further data from COHERENT as well as the from an upcoming reactor neutrino-nucleus coherent scattering experiments such as CONUS [186] can probe the most interesting part of the parameter space. In recent years, rich literature has been developed on the possibility of detecting the NSI effects in upcoming long baseline neutrino experiments. In particular, degeneracies induced by the presence of NSI in the DUNE experiment have been scrutinized [15][16][17][18][197][198][199][200][201][202][203][204]. In sect. V A, we review the effects of NSI at DUNE, T2HK and T2HKK experiments. In sect V B, we show how intermediate baseline reactor experiments such as JUNO and RENO-50 can help to determine sign(cos 2θ 12 ) and therefore test the LMA-Dark solution. In sect V C, we show how the MOMENT experiment can help to determine the octant of θ 23 and the true value of δ despite the presence of NSI. Throughout this section, we set µµ = 0 for definiteness and consistency with the majority of our references.

A. NSI at upcoming long baseline neutrino experiments
Let us first briefly review the setups of the three upcoming state-of-the-art long baseline neutrino experiments which are designed to measure the yet unknown neutrino parameters with special focus on the Dirac CP-violating phase of the PMNS matrix.
• DUNE: The source of the DUNE experiment will be at the Fermilab and the detector will be located at Sanford Underground Research Facility (SURF) at Homestake mine in South Dakota [97]. The baseline will be 1300 km. The far detector will be a 40 kton liquid Argon detector sitting on axis with the beam so the spectrum will be broad band. The energy of the neutrino beam will be around 3 GeV which comes from an 80 GeV proton beam with 1.47 × 10 21 POT per year. A reasonable assumption for data taking is 3.5 years in each neutrino and antineutrinos modes.
• T2HK: The source of T2HK [205] will be upgraded 30 GeV JPARC beam with 2.7 × 10 21 POT per year. The Hyper-Kamiokande detector with fiducial volume of 0.56 Mton [206] (25 times that of Super-Kamiokande) will be located in Kamiokande, 2.5 • off axis so the spectrum will be narrow band. The energy of neutrinos will be around 0.6 GeV. The baseline of this experiment is 295 km. A reasonable assumption for data taking is the 2TankHK-staged configuration 8 for which the data taking time is 6 years for one tank plus another 4 years with second tank [204]. The ratio of running time in neutrino mode to that in the antineutrino mode is 1:3.
• T2HKK: This project is an extension of T2HK [201] with an extra detector in Korea with a baseline of 1100 km. Two options with 2.5 • off-axis-angle and 1.5 • off-axis-angle have been discussed which respectively correspond to neutrino energies of 0.6 GeV and 0.8 GeV.
Notice that at T2HK, both energy of the neutrino beam and the baseline are lower than those at DUNE. We therefore expect the DUNE experiment to be more sensitive to both standard and non-standard matter effects than T2HK and T2HKK. Although the baseline for the Korean detector of T2HKK is comparable to the DUNE baseline, the DUNE experiment will be more sensitive to matter effects than T2HKK, because the energy of the beam at T2HKK is lower. Detailed simulation confirms this expectation [204]. In the presence of NSI, new degeneracies will appear in long baseline neutrino experiments for determination of the value of δ, mass ordering and the octant of θ 23 . One of the famous degeneracies is the so called generalized mass ordering degeneracy [11,13,19,116,204]. The oscillation probability remains invariant under the following simultaneous transformations αβ depends on the composition of medium. Notice that the LMA-Dark solution with θ 12 > π/4 and f ∼ −1 [13] is related to the generalized mass ordering transformation from the standard LMA solution with = 0. For the Earth (with N p N n ), we can write N u /N e N d /N e = 3. Notice however that the transformation in Eq. (55) does not depend on the beam energy or baseline. As a result, by carrying out long baseline neutrino experiments on Earth with different baseline and beam energy configurations, this degeneracy cannot be resolved. Resolving this degeneracy requires media with different N n /N e composition. Notice that although the nuclear compositions of the Earth core and mantle are quite different, N n /N e is uniformly close to 1 across the Earth radius [107]. However, N n /N e in the Sun considerably differs from that in the Earth. Moreover, it varies from the Sun center (with N n /N e 1/2) to its outer region (with N n /N e 1/6) [55,56]. As a result, the solar neutrino data can in principle help to solve this degeneracy. In fact, [11] by analyzing solar data shows that the LMA solution with u ee 0.3 is slightly favored over the LMA-Dark solution. The global analysis of solar, atmospheric and (very) long baseline data can in principle help to solve degeneracies. For the time being, however, since the terrestrial experiments are not precise enough to resolve the effects of sign(∆m 2 31 ) and/or sign(cos 2θ 12 ), the generalized mass ordering degeneracy cannot be resolved.
At relatively low energy long baseline experiments such as T2HK and T2HKK for which the contribution to the oscillation probability from higher orders of O(V ef f /|∆m 2 31 |) can be neglected, the appearance oscillation probability along the direction eµ / eτ = tan θ 23 will be equal to that for standard = 0 [204]. The DUNE experiment being sensitive to higher orders of (V ef f /|∆m 2 31 |) can solve this degeneracy [204]. At the DUNE experiment, another degeneracy appears when ee and τ e are simultaneously turned on and the phase of eτ is allowed to be nonzero. As shown in Fig. 12, (which corresponds to Fig 4 of [198]) and confirmed in Fig 10 of [204], in the presence of cancelation due to the phase of eτ for | eτ | ∼ 0.2 − 0.3, | ee | as large as 2 cannot be disentangled from standard case with ee = eτ = 0 at the DUNE experiment. However, this figure also demonstrates that when information on from already existent data is used as prior, the degeneracy can be considerably solved, ruling out the ee < 0 wing of solutions. That is because solar data rules out ee < 0 for θ 12 < π/4 [11]. NSI can induce degeneracies in deriving sign(cos 2θ 23 ). In principle, even with eµ ( eτ ) as small as O(0.01), the degeneracy due to the phase of eµ ( eτ ) makes the determination of the octant of θ 23 problematic [18]. To determine the sign of ∆m 2 31 , two reactor neutrino experiments with baseline of ∼ 50 km are proposed: The JUNO experiment in China which is planned to start data taking in 2020 and RENO-50 which is going to be an upgrade of the RENO experiment in South Korea. At reactor experiments, since the energy is low, |∆m 2 31 |/E √ 2G F N e . Thus, the matter effects can be neglected and the survival probability can be written as where ∆ ij = ∆m 2 ij L/(2E ν ) in which L is the baseline. Notice that the first parenthesis (which could be resolved at KamLAND) is only sensitive to sin 2 2θ 12 and cannot therefore resolve the octant of θ 12 . The terms in the last parenthesis, however, are sensitive to the octant of θ 12 . To solve these terms two main challenges have to be overcome: (i) These terms are suppressed by s 2 13 ∼ 0.02 so high statistics is required in order to resolve them. (ii) In the limit, ∆ 12 → 0, we can write P (ν e →ν e ) = c 4 13 + s 4 13 + 2s 2 13 c 2 13 cos ∆ 31 so the sensitivity to θ 12 is lost. To determine θ 12 baseline should be large enough (i.e., L 10 km). (iii) Condition ∆ 12 1 naturally implies ∆ 13 1 so the terms sensitive to the octant of θ 12 (and sign of ∆m 2 31 ) oscillate rapidly. To resolve these terms, the energy resolution and accuracy of reconstruction of the total energy scale must be high. Notice that reactor experiments such as Daya Bay satisfy the first condition and resolve the terms proportional to s 2 13 , but cannot overcome the second challenge because at these experiments, ∆ 12 1. At KamLAND, ∆ 12 > 1 but the statistics was too low to resolve the s 2 13 terms. JUNO and RENO-50, being designed to be sensitive to these terms to determine sign(∆m 2 31 ), can overcome all these Moreover, the energy calibration error can be as low as 3 %. Using the GLoBES software [207,208], Ref [19] shows how JUNO and RENO-50 experiments can test LMA-Dark solution with θ 12 > π 4 . Results are shown in Fig. 15 and Fig. 16. The star denotes the true value of ∆m 2 31 and θ 12 . In Fig. 15 and Fig. 16, normal and inverted mass orderings are respectively assumed. Ellipses show 3σ C.L. contours, after 5 years of data taking. As seen from these figures, these upcoming experiments will be able to determine |∆m 2 31 | with much better accuracy than the present global data analysis so no prior on |∆m 2 31 | is assumed. The uncertainties of other relevant neutrino parameters are taken from [209] and are treated by pull-method.
For JUNO experiment, the uncertainties in the flux normalization and the initial energy spectrum at the source are taken respectively equal to 5 % and 3 %. RENO-50 enjoys having a near detector (the detectors of present RENO) which can measure the flux with down to O(0.3%) uncertainty. To perform the analysis, the energy range of 1.8 MeV to 8 MeV is divided to 350 bins of 17.7 keV size. The pull-method is applied by defining where N i is the number of events at bin i. α i is the pull parameter that accounts for the uncertainty in the initial spectrum at bin i. Pull parameters taking care of the other uncertainties are collectively denoted by θ pull . As seen from these figures, JUNO and RENO-50 can determine the octant of θ 12 for a given mass ordering. This result is relatively robust against varying the calibration error but as expected, is extremely sensitive to the energy resolution. Increasing the uncertainty in energy resolution from 3% to 3.5%, Ref. [19] finds that JUNO and RENO-50 which stems from the generalized mass ordering degeneracy that we discussed in subsect. V A.

C. NSI at the MOMENT
The MOMENT experiment is a setup which has been proposed to measure the value of CP-violating phase, δ [210,211]. MOMENT stands for MuOn-decay MEdium baseline NeuTrino beam. This experiment will be located in China. The neutrino beam in this experiment is provided by the muon decay. Beam can switch between muon decay (µ → eν e ν µ ) and antimuon decay (μ → e + ν eνµ ). The energies of neutrinos will be relatively low with a maximum energy at 700 MeV and peak energy at 150 MeV. The detector is going to be Gd doped water Cherenkov with fiducial mass of 500 kton, located at a distance of 150 km from the source. The detection modes are Gd at the detector can capture the final neutron so although the detector lacks magnetic field, it can distinguish between neutrino and antineutrino with Charge Identification (CI) of 80 % [212]. The MOMENT experiment, with a baseline of 150 km and relatively low energy is not very sensitive to matter effects so it enjoys an ideal setup to determine δ and octant of θ 23 without ambiguity induced by degeneracies with NSI. The potential of this experiment for determining δ and the octant of θ 23 is studied in [213] using GLoBES [207,208]. The unoscillated flux of each neutrino mode is taken to be 4.7 × 10 11 m −2 year −1 and 5 years of data taking in each muon and antimuon modes is assumed. Uncertainties in flux normalization ofν e and ν µ modes are taken to be correlated and equal to 5 % but the uncertainties of fluxes from muon and antimuon modes are uncorrelated. One of the main sources of background is atmospheric neutrinos. Since the neutrino beam at the MOMENT experiment will be sent in bunches, this source of background can be dramatically reduced. Reduction of background   Fig. (a), are taken to be ∆m 2 31 = 2.417 × 10 −3 eV 2 , θ12 = 33.57 • , ∆m 2 21 = (7.45 ± 0.45) × 10 −5 eV 2 and θ13 = (8.75 ± 0.5) • . Plots are taken from [19].
is parameterized by Suppression Factor (SF). Results of Ref [213] are shown in Fig. 17. The assumed true value of δ and θ 23 are shown with a star. The mass ordering is taken to be normal and assumed to be known. All the appearance and disappearance modes are taken into account. In all these figures, the true values of are taken to be zero. In Figs 17-b to 17-d, pull method is applied on = f (N f /N e ) f , taking 1σ uncertainties on as follows [11] | eµ | < 0.16, | eτ | < 0.26 and | µτ | < 0.02 (58) and − 0.018 < τ τ − µµ < 0.054 and 0.35 < ee − µµ < 0.93.
Results shown in Fig 17-c and 17-d assume that T2K (NOνA) takes data in neutrino mode for 2 (3) years and in antineutrino mode for 6 (3) years. For more details on the assumptions, see Ref [213]. As seen from Fig 17- FIG. 16. The same as Fig. 15 except that the true values are taken to be ∆m 2 31 = 2.417 × 10 −3 eV 2 and θ12 = 56.43 • . In other words, the LMA-dark solution is assumed to be true. Plots are taken from [19].
these experiments cannot determine the octant of θ 23 . Moreover, in the presence of NSI, these experiments cannot rule out the wrong octant even at 1 σ C.L. But, comparing Fig a and b, we observe that turning on NSI within range (58,59) does not considerably reduce the power of MOMENT to measure the CP-violating phase and rule out the wrong octant solution. In this figure, SF is taken to be 0.1 % which is rather an optimistic assumption. Fig 18 shows that increasing SF up to 10 %, the power of octant determination is significantly reduced but the determination of δ is not dramatically affected.
Similar result holds valid when instead of normal mass ordering, inverted mass ordering is taken (and again assumed that the ordering is known) [213]. Ref [211] shows that the MOMENT experiment itself can determine the mass ordering. According to [213], as long as can vary in the range shown in Eqs. (58,59), MOMENT maintains its power to determine the mass ordering. Of course, once we allow to vary in a wide range such that transformation in Eq. (55) can be made, the power of mass ordering determination is lost due to the generalized mass ordering degeneracy. which are taken to be their present best fit values [104]. Both appearance and disappearance modes are taken into account. For MOMENT, SF= 0.1% . Fig (a) shows the sensitivity of MOMENT for standard scenario without NSI. In Figs (b), (c) and (d), pull method is used to treat the uncertainties of shown in Eqs. (58,59) . Fig (b) displays the sensitivity of the MOMENT experiment alone . Fig (c) shows the sensitivity of the NOνA and T2K experiments combined and Fig (d) demonstrates the combined sensitivity of all three experiments. Plots are taken from [213].

VI. SUMMARY
After multiple decades of experimental progress in the area of neutrino physics, the phenomenon of neutrino oscillations has been observed in a wide variety of experiments. The discovery of neutrino oscillations implies the existence of neutrino masses and therefore a need for an extension of the SM to include them. Many possibilities have been proposed so far, see for instance Refs. [214][215][216][217]. Motivated by the two original anomalies in the solar and atmospheric neutrino sector, various other experiments have been proposed to search for neutrino oscillations in the solar and atmospheric neutrino flux as well as in man-made neutrino beams such as reactors or accelerators. The large amount of experimental data collected over more than 20 years has allowed a very precise determination of some of the parameters responsible for the oscillations. These include the solar mass splitting, ∆m 2 21 , the absolute value of the atmospheric mass splitting, |∆m 2 31 |, as well as the solar (θ 12 ) and reactor (θ 13 ) mixing angles, mea-MOMENT, Standard Osc. sured with relative accuracies below 5%. Nevertheless, the current precision of atmospheric angle θ 23 and the CP phase δ is not at that level. The sign of cos(2θ 23 ) or in other words the octant of θ 23 is yet unknown. Moreover, although there are some hints for CP phase (δ) to be close to 3π/2, its value is not yet established. The sign of ∆m 2 31 or equivalently the scheme of mass ordering (normal versus inverted) is also still unknown. In Section II, we have discussed the most relevant experimental information used in the global fits of neutrino oscillations [34][35][36] to obtain precise measurements of the oscillation parameters, exploiting the complementarity of the different data sets. The main results of these analysis have also been commented, with an emphasis on the still unknown parameters.
Since their discovery, neutrinos have always surprised us by showing unexpected characteristics. In the dawn of the neutrino precision era, it is intriguing to ask whether neutrinos have new interactions beyond those expected within the standard model of particles. Such new interactions can give a signal in different neutrino oscillation as well as non-oscillation experiments. No evidence for the presence of NSI has been reported so far. As a consequence of these negative searches, upper bounds on the magnitude of the new interactions can be set. In Section III, we have discussed the constraints on the NSI interactions, parameterized in terms of the αβ couplings introduced in Eqs. (15) and (16). The presence of NSI has been extensively analyzed in the literature, at the level of the production, detection and propagation of neutrinos in matter. The most restrictive limits on NSI are summarized in Tables II, III and IV. In principle, adding any new particle which couples both to neutrinos and to quarks will induce Non-Standard Interaction (NSI) for neutrinos. However, it is very challenging to build an electroweak symmetric model that leads to large enough NSI to be discernible at neutrino oscillation experiments without violating various bounds. We have discussed a class of models in which the new particle responsible for NSI is a light U (1) gauge boson Z with mass 5 MeV − few 10 MeV with a coupling of order of 10 −5 − 10 −4 to quarks and neutrinos. Within this range of parameter space, the NSI effective coupling can be as large as the standard effective Fermi coupling, G F .
The total flux of solar neutrino has been measured by SNO experiment via dissociation of Deuteron through axial part of neutral current interaction and has been found to be consistent with the standard model prediction. To avoid a deviation from this prediction, the coupling of quarks to the new gauge boson is taken to be non-chiral with equal U (1) charges for left-handed and right-handed quarks. Moreover, the U (1) charges of up and down quarks are taken to be equal to make the charged current weak interaction term invariant under U (1) . The U (1) charges of quarks is taken to be universal; otherwise, in the mass basis of quarks, we would have off-diagonal interactions leading to huge q i → q j Z rate enhanced by (m qi /m Z ) 2 . In summary, U (1) charges of quarks is taken to be proportional to their baryon number. We have discussed two different scenarios for U (1) charge assignment to leptons: (i) assigning U (1) charge to the SM fermion as a e L e + a µ L µ + a τ L τ + B where L α denotes lepton flavor α and B denotes Baryon number. With this assignment, lepton flavor will be preserved and both charged leptons and neutrinos obtain lepton flavor conserving NSI. A particularly interesting scenario is a e = 0, a µ = a τ = −3/2 for which gauge symmetry anomalies automatically cancel without a need to add new specious. Choosing appropriate value of coupling [i.e., 4 × 10 −5 (m Z /10 MeV)], the best fit to solar neutrino data with µµ − ee = τ τ − ee = −0.3 can be reproduced. ii) In the second scenario, the leptons are not charged under U (1) . A new Dirac fermion, denoted by Ψ, with a mass of 1 GeV which is singlet under SM gauge group but charged under U (1) is introduced which mixes with neutrinos. As a result, neutrinos obtain coupling to Z through mixing with the new fermion but charged leptons do not couple to Z at the tree level. If the new fermion mixes with more than one flavor, both LFV and LFC NSI will be induced. Within this scenario, new fermions are needed to cancel the gauge anomalies. We have discussed different possibilities. To give masses to these new fermions, new scalars charged under U (1) are required whose VEV also gives a significant contribution to the mass of Z boson.
We have suggested two mechanism for inducing a mixing between Ψ and neutrinos: (1) Introducing a new Higgs doublet, H , with U (1) charge equal to that of Ψ which couples to left-handed lepton doublets and Ψ. H obtains a VEV of few MeV which induces mixing. (2) Introducing a sterile neutrino, N (singlet both under SM gauge group and U (1) ) and a new scalar singlet with a U (1) charge equal to that of Ψ which couples to N and Ψ. Its VEV then induces the coupling.
Even though the mass of Z particle is taken to be low (i.e., of order of solar neutrino energies and much smaller than the typical energies of atmospheric neutrinos or the energies of the neutrinos of long baseline experiments), the effect of new interaction on propagation of neutrinos in matter can be described by an effective four-Fermi Lagrangian integrating out Z because at forward scattering of neutrinos off the background matter, the energy momentum transfer is zero. At high energy scattering experiments, such as NuTeV and CHARM, the energy momentum transfer, q 2 , is much higher than m 2 Z so the effective four-Fermi coupling loses its viability. The amplitude of new effects will be suppressed by a factor of m 2 Z /q 2 1 relative to SM amplitude and will be negligible. Thus, unlike the case that the intermediate state responsible for NSI is heavy, these experiments cannot constrain ∼ 1. However, by studying scattering of low energy neutrinos (E ν ∼few 10 MeV) off matter, these models can be tested. The current COHERENT experiment and the upcoming CONUS experiment are ideal set-ups to eventually test this model. An alternative way to test such models is to search for a dip in the energy spectrum of high energy neutrinos around few hundred TeV.