Shifts of neutrino oscillation parameters in reactor antineutrino experiments with non-standard interactions

We discuss reactor antineutrino oscillations with non-standard interactions (NSIs) at the neutrino production and detection processes. The neutrino oscillation probability is calculated with a parametrization of the NSI parameters by splitting them into the averages and differences of the production and detection processes respectively. The average parts induce constant shifts of the neutrino mixing angles from their true values, and the difference parts can generate the energy (and baseline) dependent corrections to the initial mass-squared differences. We stress that only the shifts of mass-squared differences are measurable in reactor antineutrino experiments. Taking Jiangmen Underground Neutrino Observatory (JUNO) as an example, we analyze how NSIs influence the standard neutrino measurements and to what extent we can constrain the NSI parameters.


Introduction
After the observation of non-zero θ 13 from recent reactor [1,2,3,4,5] and accelerator [6,7] neutrino experiments, we have established a standard picture of three active neutrino oscillations with three mixing angles and two independent mass-squared differences [8]. Therefore the remaining neutrino mass ordering and CP-violating phase, which manifest themselves as the generic properties of three neutrino oscillations, constitute the main focus of future neutrino oscillation experiments. On the other hand, the probe of new physics beyond the Standard Model (SM) is another motivation for future precision oscillation measurements.
Experiments using reactor antineutrinos have played important roles in the history of neutrino physics, which can be traced back to the discovery of neutrinos [9], to establishment [10] of the Large Mixing Angle (LMA) Mikheyev-Smirnov-Wolfenstein (MSW) solution of the long-standing solar neutrino problem, and more recently to the discovery of non-zero θ 13 [1,2,4,5]. Moreover, future reactor experiments would keep their competitive roles in the determination of the neutrino mass hierarchy, precision measurement of oscillation parameters, and search for additional neutrino types and interactions. The survival probability for the reactor antineutrino ν e → ν e oscillation in the three neutrino framework can be written as P ee = 1 − c 4 13 sin 2 2θ 12 sin 2 ∆ 21 − c 2 12 sin 2 2θ 13 sin 2 ∆ 31 − s 2 12 sin 2 2θ 13 sin 2 ∆ 32 , with c ij = cos θ ij , s ij = sin θ ij , and ∆ ji = ∆m 2 ji L/(4E) where L is the baseline distance between the source and detector, E is the antineutrino energy, and ∆m 2 ji = m 2 j − m 2 i is the mass-squared difference between the ith and jth mass eigenstates. Because there is a large hierarchy between different masssquared differences, 30∆m 2 21 ∼ |∆m 2 31 | ∼ |∆m 2 32 | , different reactor antineutrino experiments may measure different oscillation terms of ∆ 12 or (∆ 31 , ∆ 32 ), which can be categorized into three different groups: • Long baseline reactor antineutrino experiments, such as KamLAND [10,11]. The aim of these experiments is to observe the slow oscillation with ∆ 21 and measure the corresponding oscillation parameters ∆m 2 21 and θ 12 .
• Medium baseline reactor antineutrino experiments. They stand for the next generation experiments of reactor antineutrinos, with typical representatives of Jiangmen Underground Neutrino Observatory (JUNO) [12] and RENO-50 [13]. They can determine the neutrino mass ordering (m 1 < m 2 < m 3 or m 3 < m 1 < m 2 ). In addition, they are expected to provide the precise measurement for both the fast and slow oscillations and become a bridge between short baseline and long baseline reactor antineutrino experiments.
High-dimensional operators originating from new physics can contribute to the neutrino oscillation in the form of non-standard interactions (NSIs) [14,15]. They induce effective four-fermion interactions after integrating out some heavy particles beyond the SM, where the heavy particles can be scalars, pseudo-scalars, vectors, axial-vectors, or tensors [16]. For reactor antineutrino experiments NSIs may appear in the antineutrino production and detection processes, and can modify the neutrino oscillation probability. Therefore, the neutrino mixing angles and mass-squared differences can be shifted and the mass ordering (MO) measurement will be affected. There are some previous discussions on NSIs in reactor antineutrino experiments [17,18,19] and other types of oscillation experiments [20]. In this work, we study the NSI effect in reactor antineutrino oscillations in both specific models and also the most general case. Taking JUNO as an example, we apply our general framework to the medium baseline reactor antineutrino experiment. We discuss how NSIs influence the standard 3-generation neutrino oscillation measurements and to what extent we can constrain the NSI parameters.
The remaining part of this work is organized as follows. Section 2 is to derive the analytical formalism. We develop a general framework on the NSI effect in reaction antineutrino oscillations, and calculate the neutrino survival probability in the presence of NSIs. In section 3, we give the numerical analysis for the JUNO experiment. We analyze the NSI impacts on the precision measurement of mass-squared differences and the determination of the neutrino mass ordering, and present the JUNO sensitivity of the relevant NSI parameters. Finally, we conclude in section 4.

NSI-induced neutrino oscillations 2.1 Basic formalism
NSIs may occur in the neutrino production, detection and propagation processes in neutrino oscillation experiments. The neutrino and antineutrino states produced in the source and observed in the detector are superpositions of flavor states, in which the superscripts 's' and 'd' denote the source and detector, respectively, and are normalization factors. In general, NSIs in different physical processes may have distinct contributions. For a certain type of neutrino experiments, the same set of effective NSI parameters can be introduced to describe the NSI effect. But when one turns to anther type of neutrino experiments, neutrinos can have totally different origins and one should use another set of NSI effective parameters to parametrize the NSI effect. The parameters s,d αβ used here are strongly experiment-dependent, and in principle also energy-dependent. However, they are usually considered as the averaged effects and treated as constant values.
In order to measure the average and difference between neutrino production and detection processes, we introduce two sets of NSI parameters as to rewrite the NSI effect. One should note that the NSI parameters δ αβ are negligibly small compared with˜ αβ when neutrinos are purely left-handed particles [16,20]. However, if neutrinos are right-handed particles and the four-fermion interactions are mediated by heavy scalars beyond SM, δ eα can reach the percent level [16]. The effective Hamiltonian that describes the vacuum neutrino oscillation is given by where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [21] and can be expressed in the form [22] where s ij = sin θ ij , P l = diag{e iφ 1 , e iφ 2 , e iφ 3 }, and P ν = diag{1, e iρ , e iσ }. P l is unphysical since charged leptons are Dirac particles. On the other hand, P ν can not be neglected for Majorana neutrinos but does not contribute to neutrino oscillations. We write the PMNS matrix in the form as in Eq. (7) to keep the first row and the third column of P l U P † ν real. The mixing angles θ 12 , θ 13 , θ 23 and CP phase δ take the same values as those in the PDG parametrization [8], respectively.
In the presence of NSI effects at the source and detector, the amplitude of the ν α → ν α transition is where L is the baseline, and is the amplitude of ν β → ν γ without NSIs. It is useful to definẽ and i |Ũ αi | 2 = 1 is required. Thus we can obtainÃ αα as and the survival probability for ν α → ν α is expressed as Because of the smallness of δU αi /Ũ αi , we can rewrite the above equation as with∆ α ji = ∆ ji + Im Note that∆ α 31 =∆ α 32 +∆ α 21 holds, just as the relation in the standard three neutrino case with ∆ 31 = ∆ 32 + ∆ 21 . For the effective mass-squared differences we have ∆m 2α ji (E/L) = ∆m 2 ji + Im which is an energy/baseline-and flavor-dependent effective quantity. With Eq. (14) we have obtained a standard-like expression for the antineutrino survival probability in vacuum in the presence of NSIs. The corresponding NSI effects are encoded in the effective mass and mixing parameters. The average parts induce constant shifts for the neutrino mixing elements, and the difference parts generate energy and baseline dependent corrections to the mass-squared differences. Without prior information on the true mass and mixing parameters, we cannot distinguish between the true and effective parameters. Only the energy-and baseline-dependent feature of the effective parameters can tell us the NSI effects, and correspondingly the difference parts of NSIs will be constrained with precise spectral measurements of reactor antineutrino oscillations.

Parametrization of the effective PMNS matrixŨ
We have defined the matrixŨ in Eq. (10), which is considered as an effective PMNS matrix compared with the original one U in the standard three neutrino framework. However, the unitary conditions for U do not hold any more. Instead, we have We plan to perform a similar parametrization forŨ as in Eq. (7). However,Ũ is not a unitary matrix and to perform such a parametrization is not a trivial task. Taking account of the relation in Eq. (17) and the detailed configuration of neutrino oscillation experiments, we parametrizeŨ as follows, with s ij = sinθ ij andc ij = cosθ ij for simplicity, 1) We first extract two diagonal phase matricesP l = diag{e iφ 1 , e iφ 2 , e iφ 3 } andP ν = diag{1, e iρ , e iσ } to make the first row and third column ofP lŨP † ν real and positive.P l includes unphysical phases to be redefined by rephasing the charged lepton fields, butρ andσ are the effective Majorana CP phases modified by NSIs, which do not contribute to the neutrino oscillation but will essentially influence the process of the neutrinoless double beta decay.

5)
Finally, there are 4 additional CP-violating phasesα µ1 ,α µ2 ,α τ 1 andα τ 2 which are defined as follows: To summarize, we have obtained a parametrization similar to Eq. (7): Fifteen parameters are included inŨ and ten of them (θ 12 ,θ 13 ,θ 23 ,θ 23 ,δ,δ ,α µ1 ,α µ2 ,α τ 1 ,α τ 2 ) are related to neutrino oscillations. In the case of˜ → 0, we can restore the original PMNS matrix of three active neutrinos as One should notice that although the effective PMNS matrix takes a simple form as in Eq. (19), the corresponding parameters are in general dependent on the type of experiments. We should use different U to characterize different realizations of the effective PMNS matrix for reactor and accelerator neutrino experiments. This is different from the non-unitary effect of the PMNS matrix, where˜ is used to parametrize the universal mixing between the active and sterile neutrinos and δ = 0 by definition. In this case,Ũ is an effective PMNS matrix for all neutrino oscillation experiments.

Reactor antineutrino oscillation probabilities
In reactor antineutrino oscillations, only the electron antineutrino survival probability is relevant because of the high threshold of the µ/τ production. With the parametrization ofŨ in Eq. (19), we can rewrite P ee with these effective mixing parameters as with∆ where the superscript α = e in∆ e ji has been ignored. The average part˜ can be treated as constant shifts to mixing angles θ 12 and θ 13 , and the difference part δ leads to energy-and baseline-dependent shifts to the mass-squared differences ∆m 2 ji as However, only two combinations of the three parameters δε i contribute to the oscillation probability thanks to the relation ( . It is notable that one cannot distinguish the effect of mixing angle shifts from the scenario of three neutrino mixing using reactor antineutrino oscillations. This degeneracy can only be resolved by including different types of neutrino oscillation experiments, where the NSI parameters and their roles in neutrino oscillations are totally distinct. On the other hand, the shifts of mass-squared differences are clearly observable due to the baseline-and energy-dependent corrections in the reactor antineutrino spectrum. Different kinds of reactor antineutrino oscillation experiments have their own advantages in measuring the NSI parameters. The long baseline reactor antineutrino experiment (e.g., KamLAND) can measure the slow oscillation term ∆ 21 and thus are sensitive to the measurement of δε 1 − δε 2 . Since the fast oscillation terms ∆ 31 and ∆ 32 are averaged out, the oscillation probabilityP ee is reduced tõ The short baseline reactor antineutrino experiments (e.g., Daya Bay) are designed to measure the fast oscillation terms ∆ 31 and ∆ 32 (or equivalently, ∆ ee ), and are effective to constrain δε 1 − δε 3 . Since the slow oscillation term ∆ 21 is negligible, the oscillation probabilityP ee can be simplified as P ee = 1 − sin 2 2θ 13 sin 2 ∆ ee − η c 2 12 (δε 1 − δε 3 ) +s 2 12 (δε 2 − δε 3 ) sin 2 2θ 13 cos 2∆ ee .
For the reactor antineutrino oscillation at the medium baseline (e.g., 52.5 km), we can generalize the formalism given in Ref. [23] to see how the mass ordering measurement can be influenced by the NSI effect. The oscillation probability in Eq. (21) can be rewritten as where with η = ±1 for the normal mass ordering (NMO) and inverted mass ordering (IMO), respectively. As a function of δε 1 − δε 3 and δε 2 − δε 3 , S varies with the NSI parameters, which can further alter the difference of oscillation probabilities between NMO and IMO. Since NSIs are constrained to the percent level [24], we expect that NSIs will not significantly affect the mass ordering measurement. Finally, one should keep in mind that there is an additional correction from terrestrial matter effects during the neutrino propagation inside the Earth. In addition to the Hamiltonian in Eq. (6), there is a matter potential term written as where −A CC characterizes the contribution of charged-current interactions between antineutrinos and electrons in matter and A CC = 2 √ 2G F N e E, with N e being the electron number density. For reactor antineutrino experiments, A CC is sufficiently small compared with the kinetic term ∆m 2 ji where the matter corrections to oscillation parameters of the solar sector are given by and the magnitude of these corrections is estimated as with ρ being the matter density along the antineutrino trajectory in the Earth. In comparison, the correction to parameters of the atmospheric sector is only the order of 10 −5 . In our following studies, matter effects during the neutrino propagation will be neglected in the analytical calculation but effectively included in our numerical analysis.
The N e E suppression leads to negligible NSI effects in the propagation process. The scenario is very different from other types of neutrino oscillation experiments, where NSIs during propagation lead to much larger corrections to oscillation probabilities than those at the source and detector. In reactor antineutrino oscillations NSIs during propagation can be safely neglected.

Numerical simulation
In this section, following the JUNO nominal setup in Ref. [12], we will numerically show the NSI effect in the medium baseline reactor antineutrino experiment. We will illustrate how the NSI effect shifts the mass-squared differences, how it influences the mass ordering measurement and to what extent we can constrain the NSI parameters.

Nominal setups
We employ the true power and baseline distribution in Table 1 of Ref. [12]. The weighted average of the baseline is around 52.5 km with baseline differences less than 500 m. We use the nominal running time of 1800 effective days for six years, and detector energy resolution 3%/ E(MeV) as a benchmark. A NMO is assumed to be true otherwise mentioned explicitly. The relevant oscillation parameters arẽ θ 12 ,θ 13 , ∆m 2 21 and ∆m 2 ee and the NSI parameters are δε 1 − δε 2 , δε 1 − δε 3 .
The measured mixing angles θ D 13 and θ K 12 are from Daya Bay [25] and KamLAND [26], respectively. For some recent discussions on how the true mixing angles are modified by the average of the NSI effect in a simplified version, see [18,27]. However, the measured mass-squared differences ∆m 2 D ee from Daya Bay [25] and ∆m 2 K 21 from KamLAND [26] can be viewed as the true parameters rather than effective oscillation parameters, i.e., ∆m 2 ee = ∆m 2 D ee = 2.44 × 10 −3 eV 2 , ∆m 2 21 = ∆m 2 K 21 = 7.54 × 10 −5 eV 2 .
The reason is that although the mass-squared differences have energy/baseline-dependent corrections as shown in Eq. (23), these shifts are sufficiently small compared with the current level of uncertainties. In details, δε i may be in the percent level and the ratio E/L for Daya Bay is around (E/L) D 2 MeV/km 4 × 10 −4 eV 2 , and for KamLAND is around (E/L) K 0.02 MeV/km 4 × 10 −6 eV 2 . Therefore the shifts are still below the sensitivity of Daya Bay [25] and KamLAND [26], and the measured parameters can be approximate to the true values for ∆m 2 ee and ∆m 2 21 , respectively. In our following numerical analysis, two different treatments for the NSI parameters will be explored.
• The second one is the general treatment for the NSI parameters. As shown in Eq. (21), only two combinations of the NSI parameters δε i are independent, thus we can treat (δε 1 − δε 2 ) and (δε 1 − δε 3 ) as free parameters to cover the full parameter space. Note that the non-unitary effect is identical with NSIs for the case of δ αβ = 0. If we can observe the splittings of∆ α ji compared with ∆ ji , we may distinguish NSIs from the non-unitary effect.

Antineutrino spectrum
To calculate the expected reactor ν e spectrum in the presence of NSIs, we need first deal with the standard case without the NSI effect. The energy spectrum of detected events S(E vis ), as a function of the visible energy E vis of the inverse β-decay ν e + p → e + + n (IBD), is parametrized as where Φ i (E) is the antineutrino flux with i standing for different reactor cores and E the antineutrino energy, N i is the corresponding normalization and conversion factor, P ee (E/L i ) is the oscillation probability of ν e → ν e with different baseline L i from the ν e source i to the detector, dσ(E, E e )/dE e is the IBD differential cross section with E e being the true positron energy, r(E e + m e , E vis ) is the Gaussian energy resolution function with a standard deviation σ E defined as (E e + m e )/MeV .
In the presence of NSIs, we can use the following replacements to show the NSI effect during the antineutrino oscillation, production and detection processes, respectively Therefore, we can obtain the NSI-modified reactor ν e spectrum as where N s e , N d e andÑ e are defined in Eqs. (4) and (11), and can be absorbed by redefining the couplings of nuclear matrix elements in the reactor ν e production and detection processes. From Eq. (11),Ñ e is related to the average parts of NSIs, and contributes to the normalization factor of the reactor antineutrino flux. In this work we takeÑ e = 1.0 for simplicity, as we mainly stress on the difference parts of NSIs and the experimental spectral measurements.
We show the effect of NSIs in the reactor ν e spectra at a baseline of 52.5 km in Fig. 1, where the influences of δε 1 − δε 2 and δε 1 − δε 3 are presented in the upper panel and lower panel, respectively. The oscillation parameters are taken as in Eqs. (32) and (33). The scenario of 3-generation neutrino oscillations with δU = 0 is also shown for comparison. In the upper panel, we fix δε 1 − δε 3 = 0 and find that non-zero δε 1 − δε 2 introduces the spectral distortion to the slow oscillation term ∆ 21 . For δε 1 − δε 2 = 0.02, the spectrum is suppressed in the high energy region with E > 3 MeV and enhanced for the low energy range 2 MeV < E < 3 MeV. In comparison, negative δε 1 − δε 2 gives the opposite effect on the spectrum distortion. In the lower panel, we set δε 1 − δε 2 = 0 and observe that δε 1 − δε 3 can affect the spectral distribution for the fast oscillation term ∆ 31 . Non-trivial NSI effect will contribute a small phase advancement or retardance to the fast oscillation depending upon the sign of δε 1 − δε 3 .

Statistical analysis
In this part, we shall implement the standard χ 2 statistical method to do the numerical analysis with the above setup. A general χ 2 function using the spectrum calculated in Eq. (38) can be defined as whereM i andT i are the measured and predicted (with oscillation) reactor ν e fluxes in the i-th energy bin respectively. The definition of bin sizes is identical to that assumed in Ref. [12]. The systematic uncertainties σ k together with the corresponding pull parameters k for the nominal setups are also the same as those in Ref. [12], which include the correlated (absolute) reactor uncertainty (2%), the uncorrelated (relative) reactor uncertainty (0.8%), the flux spectrum uncertainty (1%) and the detectorrelated uncertainty (1%). The sets of p and δε are for the standard oscillation parameters and NSI parameters respectively with p = {θ 12 ,θ 13 , ∆m 2 12 , ∆m 2 ee } and, δε = {δU/c 12c13 , ±δU/s 12c13 , ±δU/s 13 } for specified models, or δε = {δε 1 − δε 2 , δε 1 − δε 3 } for the general model defined in section 3.1.

Shifts of mass-squared differences
Neglecting the existing non-zero NSIs, we may get biased best-fit oscillation parameters. In this part we shall evaluate the sizes and properties on the shifts of mass-squared differences due to the NSI effect. In the numerical simulation, we use the spectrum with non-zero NSIs as the true spectrum, and the spectrum of the standard neutrino oscillation without NSIs as the predicted spectrum. In other word, Figure 2: The NSI-induced shifts for ∆m 2 21 and ∆m 2 31 in the four specific models defined in Eq. (34) with δU = 0.01. The red stars stand for the true values of ∆m 2 21 and ∆m 2 31 , and contours are the 68.3% (1σ), 95.5% (2σ), 99.7% (3σ) allowed regions for ∆m 2 21 and ∆m 2 31 when the NSI effect is neglected. The NMO is assumed for illustration. the true spectrum is defined in Eq. (38) with the oscillation probability given in Eq. (21), and the predicted spectrum is given in Eq. (35) with the oscillation probability in Eq. (1). Then we minimize the χ 2 function and find out the best-fit mixing angles and mass-squared differences.
The effects of mass shifts are shown in Figs. 2 and 3, corresponding to the first and second treatments, respectively, where NMO has been assumed. In the first treatment, the mass-squared differences have different shift sizes and directions in each specific models, dependent upon the sign of δU . The best-fit value of ∆m 2 21 decreases from its true value in the models of S1, S2 or increases for S3, S4, and the best-fit value of ∆m 2 31 decreases in S1, S3 or increases in S2, S4. A simple explanation can be found in the following estimation. The NSI parameters in these models are numerically given by S1 : δε 1 − δε 2 = −0.54 δU, δε 1 − δε 3 = −5.60 δU, S2 : δε 1 − δε 2 = −0.54 δU, δε 1 − δε 3 = +8.06 δU, S3 : δε 1 − δε 2 = +3.00 δU, δε 1 − δε 3 = −5.60 δU, where δU is fixed at 0.01 in Fig. 2. The sign of δε 1 − δε 2 is "−" in S1, S2 and "+" in S3, S4, which reduces or increases the measured value of ∆m 2 21 in the LHS of Eq. (23), respectively. Similar analysis is also valid for explaining the shift of ∆m 2 31 as shown in Fig. 2. Moreover, due to the smallness of the coefficients of δε 1 − δε 2 in models S1, S2, the shift of ∆m 2 21 in S1, S2 is much smaller than that in S3, S4. The relative mass shift for ∆m 2 21 is around 0.4% in S1, S2 and 2% in S3, S4. Although the magnitude of δε 1 − δε 3 is in general larger than δε 1 − δε 2 , the absolute value of ∆m 2 31 is much larger than ∆m 2 21 , and thus the relative shift of ∆m 2 31 is not significant, just roughly around 0.2% for the best-fit data in four models.
The effects of NSI-induced mass shifts in the general case are presented in Fig. 3, without any assumptions on the relation of NSI parameters. Our simulation results can be understood using the relation in Eq. (23) that the fitted mass-squared differences ∆m 2 J 21 and ∆m 2 J 31 in JUNO are linearly dependent upon δε 1 − δε 2 and δε 1 − δε 3 , respectively. With the JUNO nominal setup, we can simplify Eq. (23) into the following formulae where 4(E/L) J 5 × 10 −5 eV 2 roughly holds. The relative mass shift of ∆m 2 21 is about 2 3 (δε 1 − δε 2 ), in the same level of δε 1 − δε 2 , i.e., in the same order as Im(δU ei ) and δ αβ . The relative mass shift of ∆m 2 31 is about 0.02(δε 1 − δε 3 ). Keeping in mind δε 3 = Im(δU ei )/s 13 , we obtain that the relative mass shift of ∆m 2 31 is one order smaller than Im(δU ei ) and δ αβ . In the case of the IMO, as ∆m 2 31 = −|∆m 2 31 | holds, the shift of |∆m 2 31 | will go in the opposite direction to the NMO.  The iso-∆χ 2 contours for the MO sensitivity in the generic NSI model as a function of two effective NSI parameters δε 1 − δε 2 and δε 1 − δε 3 . The NMO is assumed for illustration.

Impacts on the MO measurement
When fitting the χ 2 function in Eq. (39) with both NMO and IMO, we can take the difference of the minima to measure the sensitivity of neutrino mass ordering, where the discriminator is defined as For these specific models defined in Eq. (34), we illustrate in Fig. 4 the NSI effect on the MO measurement by showing the dependence of the MO sensitivity on the true value of the NSI parameter δU , where the NMO is assumed for illustration. The NSI effect with a negative δU in S1, S3 or positive δU in S2, S4 will decrease the ∆χ 2 (MO) value and thus degrade the sensitivity of the MO determination. However in the other half possibilities, the NSI effect can increase the ∆χ 2 (MO) value and enhance the MO sensitivity. Moreover, the NSI effect shows stronger influence on the MO measurement in models S2, S3 than S1, S4. On the other hand, we illustrate in Fig. 5 the iso-∆χ 2 contours for the MO sensitivity in the generic NSI model as a function of two effective NSI parameters δε 1 − δε 2 and δε 1 − δε 3 . The NMO is assumed for illustration. We can learn from the figure that the smaller δε 1 − δε 2 and larger δε 1 − δε 3 will reduce the possibility of the MO measurement. If δε 1 − δε 2 decreases by 0.03 or δε 1 − δε 3 increases by 0.05, ∆χ 2 will be suppressed by 2 units.

Constraints on NSI parameters
In this part we shall discuss the constraints on the NSI parameters with the JUNO nominal setup. In our numerical calculation, the true oscillation parameters are taken as in Eqs. (32) and (33), and the true NSI parameters are taken as δε 1 − δε 2 = δε 1 − δε 3 = 0. In the fitting process, we fix the oscillation parameters but take the NSI parameters as free. With the above simplification, we can obtain the constraints on the considered NSI parameters. In Fig. 6, we show the limit on these two parameters at the 1, 2, 3σ confidence levels. For δε 1 − δε 2 , the precision is much better than 1%. However, the precision for δε 1 − δε 3 is around the 10% level. JUNO is designed for a precision spectral measurement at the oscillation maximum of ∆m 2 21 . From Eq. (21), the precision for δε 1 − δε 2 can be compatible with that of sin 2 2θ 12 , where a sub-percent level can be achieved [28]. On the other hand, the precision for sin 2 2θ 13 is also at the 10% level, also consistent with that of δε 1 − δε 3 in our numerical simulation. Because δε 1 − δε 3 is suppressed by sinθ 13 , the above two constraints are actually compatible if we consider the physical NSI parameters δ αβ defined in Eq. (5). Notice that different assumptions (e.g., the uncertainties of oscillation parameters) on the experimental systematics may alter the quantitative precision of the NSI parameters, but our qualitative conclusion is reasonable in any realistic systematical assumptions.

Conclusion
In this work we have presented a complete and new derivation on the generic NSI effects in reactor antineutrino oscillations, where the NSI parameters are divided into the average and difference parts of the antineutrino production and detection processes. The average part can induce an effective nonunitary PMNS matrix and shift the true values of the mixing angles. On the other hand, the difference part of the NSI effect can be parametrized with only two independent parameters (i.e., δε 1 − δε 2 , δε 1 − δε 3 ), and give the energy-and baseline-dependent corrections to the mass-squared differences. Eq. (21) is our key formula for the reactor ν e → ν e survival probability, where Figure 6: The experimental constraints on the generic NSI parameters δε 1 − δε 2 and δε 1 − δε 3 , where the true values are fixed at δε 1 − δε 2 = δε 1 − δε 3 = 0, and the contours are the 68.3% (1σ), 95.5% (2σ), 99.7% (3σ) allowed regions. The NMO is assumed for illustration.
• we define the mixing angle shifts as the deviations of measured mixing anglesθ 12 ,θ 13 from their true values θ 12 , θ 13 . However, we stress that these constant shifts are undetectable in reactor antineutrino experiments.
• the two NSI parameters δε 1 − δε 2 and δε 1 − δε 3 can be absorbed into the mass-squared differences and the corresponding E/L-dependent effective parameters ∆m 2 12 and ∆m 2 31 can be defined as the shifts of mass-squared differences. These shifts are detectable in the spectral measurement of reactor antineutrino oscillations.
Our analytical formalism is applied to the future medium baseline reactor antineutrino experiment JUNO. Two different treatments (a class of specific models and the most general case with the full parameter space) of the NSI parameters are employed in our numerical analysis. We analyze the NSI impact on the precision measurement of mass-squared differences and the determination of the neutrino mass ordering, and present the JUNO sensitivity of the relevant NSI parameters. Numerically, • we find that the relative mass shift of ∆m 2 21 is around 2 3 (δε 1 − δε 2 ), in the same order of the original NSI parameters δ αβ ; and the relative shift of ∆m 2 31 is around 0.02(δε 1 − δε 3 ), one order smaller than the magnitude of δ αβ . However, cancelations may appear in δε 1 − δε 2 and suppress the mass shift of ∆m 2 21 (see the models S1 and S2).
• a positive δε 1 − δε 2 or negative δε 1 − δε 3 may enhance the sensitivity of the neutrino MO measurement at JUNO.
• due to the specific configuration of JUNO, the constraint on δε 1 − δε 2 can be better than 1%, but δε 1 − δε 3 can only be measured at the 10% precision level.
Compared with long baseline and short baseline reactor antineutrino experiments (e.g., KamLAND and Daya Bay), the medium baseline reactor antineutrino experiment JUNO is more suitable for constraining the NSI effect because both the slow and fast oscillation terms are measurable in the reactor antineutrino spectral measurement. Taking into account the complementary properties of reactor antineutrino experiments at different baselines, it is desirable to present a sophisticated global analysis of all the reactor antineutrino experiments and therefore, we may obtain the most complete and precision testing of the NSIs or other new physics beyond the SM.