Running of oscillation parameters in matter with flavor-diagonal non-standard interactions of the neutrino

In this article we unravel the role of matter effect in neutrino oscillation in the presence of lepton-flavor-conserving, non-universal non-standard interactions (NSI’s) of the neutrino. Employing the Jacobi method, we derive approximate analytical expressions for the effective mass-squared differences and mixing angles in matter. It is shown that, within the effective mixing matrix, the Standard Model (SM) W -exchange interaction only affects θ12 and θ13, while the flavor-diagonal NSI’s only affect θ23. The CP-violating phase δ remains unaffected. Using our simple and compact analytical approximation, we study the impact of the flavor-diagonal NSI’s on the neutrino oscillation probabilities for various appearance and disappearance channels. At higher energies and longer baselines, it is found that the impact of the NSI’s can be significant in the νμ → νμ channel, which can probed in future atmospheric neutrino experiments, if the NSI’s are of the order of their current upper bounds. Our analysis also enables us to explore the possible degeneracy between the octant of θ23 and the sign of the NSI parameter for a given choice of mass hierarchy in a simple manner.

While the three-flavor oscillation probabilities in matter can be calculated numerically on a computer, with or without NSI's, and properly taking into account the changing mass-density along the baseline, the computer program is a blackbox which does not offer any deep understanding as to why the probabilities depend on the input parameters in a particular way. If the mass-density along the baseline is approximated by an average constant value, then exact analytical expressions for the three-flavor oscillation probabilities can, in principle, be derived as was done for the SM case in refs. [56][57][58][59]. But even for the SM case, the expressions are extremely lengthy and too complicated to yield much physical insight. Thus, for the SM case, further approximations which simplify the analytic expressions while maintaining the essential physics have been developed by various authors [60][61][62][63][64][65][66][67][68][69][70][71] to help us in this direction. Indeed, approximate analytical expressions of the SM neutrino oscillation probabilities in constant-density matter have played important roles in understanding the nature of the flavor transitions as functions of baseline L and/or neutrino energy E [65][66][67]69]. To obtain similar insights for the NSI case, approximate analytical expressions for the three-flavor oscillation probabilities in constant-density matter in the presence of NSI's are called for. Once they have provided us with the intuition we seek, on how and why the oscillation probabilities behave in a particular way, we can then resort to numerical techniques to further refine the analysis, e.g. taking into account the non-constant mass-density, if the need arises. 1 For a recent review, see ref. [20]. 2 There are two possibilities: it can be either 'normal' if δm 2 31 ≡ m 2 3 − m 2 1 > 0, or 'inverted' if δm 2 31 < 0. 3 Present status and future prospects of NSI's are discussed in recent reviews [38,39].

JHEP11(2015)035
In previous works [68,70], we looked at the matter effect on neutrino oscillation due to the SM W -exchange interaction between the matter electrons and the propagating electron neutrinos. Employing the Jacobi method [72], we showed that the said matter effects could be absorbed into the 'running' of the effective mass-squared-differences, and the effective mixing angles θ 12 and θ 13 in matter as functions of the parameter a = 2 √ 2G F N e E, while the effective values of θ 23 and the CP-violating phase δ remained unaffected. Here, G F is the Fermi muon decay constant, N e is the electron density, and E is the energy of the neutrino. The approximate neutrino oscillation probabilities were obtained by simply replacing the oscillation parameters in the vacuum expressions for the probabilities with their running in-matter counterparts.
This running-effective-parameter approach has several advantages over other approaches which approximate the neutrino oscillation probabilities directly. First, the resulting expressions for the probabilities are strictly positive, which is not always the case when the probabilities are directly expanded in some small parameter, and the series truncated after a few terms. Second, the behavior of the oscillation probabilities as functions of L and E can be understood easily as due to the running of the oscillation parameters with a, as was shown in several examples in refs. [68,70]. Third, as will be shown later, it is very convenient in exploring possible correlations and degeneracies among the mass-mixing parameters that may appear in matter in a non-trivial fashion.
In this paper, we extend our previous analysis and investigate how our conclusions are modified in the presence of neutrino NSI's of the form where subscripts α, β = e, µ, τ label the neutrino flavor, f = e, u, d mark the matter fermions, C = L, R denotes the chirality of the f f current, and ε f C αβ are dimensionless quantities which parametrize the strengths of the interactions. The hermiticity of the interaction demands ε f C βα = (ε f C αβ ) * . (1.2) For neutrino propagation through matter, the relevant combinations are (1.3) where N f denotes the density of fermion f . In this current work, we limit our investigation to flavor-diagonal NSI's, that is, we only allow the ε αβ 's with α = β to be non-zero. The case of flavor non-diagonal NSI's will be considered in a separate work [73]. In the following, we will show how the presence of such flavor-diagonal NSI's affect the running of the effective neutrino oscillation parameters (the mass-squared differences, mixing angles, and CP-violating phase), and ultimately how they alter the oscillation probabilities. We find that due to the expected smallness of the ε αα 's as compared to the SM W -exchange interaction, there is a clear separation in the ranges of a at which the NSI's and the SM interaction are respectively relevant. Furthermore, within the effective -3 -JHEP11(2015)035 1. Models with a generation distinguishing Z boson. This class includes gauged L e −L µ and gauged L e −L τ [91,92], gauged B −αL e −βL µ −γL τ (with α+β +γ = 3) [89,[93][94][95][96][97], and topcolor assisted technicolor [87,98].
Systematic studies of how NSI's can arise in these, and other BSM theories can be found, for instance, in refs. [38,85,86,135]. Thus, the ability to detect NSI's in neutrino experiments would complement the direct searches for new particles at the LHC for a variety of BSM models. Next, we focus our attention to see the role of lepton-flavor-conserving NSI's when neutrinos travel through the matter.

Lepton-flavor-conserving NSI's
The NSI's of eq. (1.3) modify the effective Hamiltonian for neutrino propagation in matter in the flavor basis to where U is the vacuum Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [136][137][138], E is the neutrino energy, the matter-effect parameter a is given by For Earth matter, we can assume N n ≈ N p = N e , in which case N u ≈ N d ≈ 3N e . Therefore, Since we restrict our attention to lepton-flavor-conserving NSI's in this paper, the effective Hamiltonian in the flavor basis takes the form

JHEP11(2015)035
In the absence of off-diagonal terms, we can rewrite the matter effect matrix as follows: where we have assumed that the ε's are small compared to 1. Let Note that the NSI parameter η introduced here is related to the parameter ξ which was used in ref. [82] by Thus, the effective Hamiltonian can be taken to be (2.9) Note that the η = 0 case simply replaces a withâ = a(1 + ζ) in the SM Hamiltonian. Thus, the effective mass-squared differences, mixing angles, and CP-violating phase have the same functional dependence onâ as they had on a in the SM case. In other words, a non-zero ζ simply rescales the value of a by a constant factor. The presence of a non-zero η, on the other hand, requires us to diagonalize H with a different unitary matrix from the η = 0 case, and this will introduce corrections to the effective oscillation parameters beyond a simple rescaling of a. Next, we discuss the constraints that we have at present on these NSI parameters.

Existing bounds on the NSI parameters
Before looking at the matter effect due to the neutrino NSI's, let us first look at what is currently known about the sizes of the ε f C αβ 's and the combinations η = (ε µµ − ε τ τ )/2 and ζ = ε ee − (ε µµ + ε τ τ )/2.

Theoretical expectation
Theoretically, the sizes of the NSI's are generically expected to be small since they are putatively due to BSM physics at a much higher scale than the electroweak scale, or to loop effects. If they arise from the tree level exchange of new particles of mass Λ, which would be described by dimension six operators, we can expect the ε f C αβ 's to be of order . Loop effects involving a heavy particle of mass Λ would be further suppressed by a factor of O(1/4π) or more. Processes that lead to dimension eight operators would be of order O(M 4 W /Λ 4 ). Thus, if we assume Λ = O(1 TeV), we expect the ε f C αβ 's, and consequently η and ζ, to be O(10 −3 ) or smaller.

JHEP11(2015)035
A tighter bound exists for η which has been obtained directly using solar and atmospheric neutrino data in ref. [174], and more recently in ref. [193] using atmospheric and MINOS data. In ref. [174], only NSI's with the d-quarks were considered and the following 3σ (99.7% C.L.) bounds were obtained: [174] is actually constraining ε αβ /3, so this result can be interpreted as Rescaling to 1.64σ (90% C.L.), we find Ref. [193] gives slightly different 90% C.L. bounds of 3 Effective mixing angles and effective mass-squared differences -Neutrino case

Setup of the problem
As we have seen, in the presence of non-zero η and ζ, the effective Hamiltonian (times 2E) for neutrino propagation in Earth matter is given by Table 1. Second column shows the best-fit values and 1σ uncertainties on the oscillation parameters taken from ref. [194]. We use the values listed in the third column as benchmark values for which we calculate our oscillation probabilities in this work. whereâ = a(1+ζ). The problem is to diagonalize H η = H 0 +âηM η and find the eigenvalues λ i (i = 1, 2, 3) and the diagonalization matrix ∼ U as functions ofâ and η. To this end, we utilize the method used in refs. [68,70] where approximate expressions for the λ i 's and ∼ U were derived for H 0 , the η = 0 case, using the Jacobi method [72]. The Jacobi method entails diagonalizing 2 × 2 submatrices of a matrix in the order which requires the the largest rotation angles until the off-diagonal elements are negligibly small. In the case of H 0 , it was discovered in refs. [68,70] that two 2 × 2 rotations were sufficient to render it approximately diagonal, and that these two rotation angles could be absorbed into 'running' values of θ 12 and θ 13 . The procedure that we use in the following for H η is to add on theâηM η term to H 0 after it is approximately diagonalized, and then proceed with a third 2 × 2 rotation to rotate away the additional off-diagonal terms.
As the order parameter to evaluate the sizes of the off-diagonal elements, we use and consider H 0 and H η to be approximately diagonalized when the rotation angles required for further diagonalization are of order 3 = 0.005 or smaller. Note that we are using a slightly different epsilon ( ) here to distinguish this quantity from the NSI's (ε αβ ).
Unless otherwise stated, we use the benchmark values of the various oscillation parameters as given in the third column of table 1 to draw our plots. These values are taken from ref. [194] and correspond to the case in which reactor fluxes have been left free in the fit and short-baseline reactor data with L ≤ 100 m are included. For sin 2 θ 23 , we consider the benchmark value which lies in the lower octant and CP-violating phase δ is assumed to be zero. These choices of the oscillation parameters are well within their 3σ allowed ranges which are obtained in recent global fit analyses [16][17][18]. We also present results considering other allowed values of sin 2 θ 23 and δ which we discuss in section 3.5. In the evaluation of the sizes of the elements of the effective Hamiltonian, we will assume θ 13 ≈ 0.15 = O( ), cos(2θ 12 )/2 ≈ 0.2 = O( ), and | cos(2θ 23 )| ≈ 0.18 = O( ). We also assume that the NSI parameter η is of order 2 = 0.03 (or smaller), since the current 90% C.L. (1.64σ) upper bound on |η| was ∼ 0.05, cf. eqs. (2.16) and (2.18), though we allow it to be as large as 0.1 in our plots to enhance and make visible the effect of a non-zero η.

Change to the mass eigenbasis in vacuum
Introducing the matrix we begin by partially diagonalizing the Hamiltonian H η as and M η (θ 12 , θ 13 , θ 23  Using cos(2θ 23 ) = O( ) and θ 13 = O( ), we estimate the sizes of the elements of M a to be: and those of M η to be Given that we have assumed η = O( 2 ) or smaller, the off-diagonal elements of ηM η are suppressed compared to those of M a , and only become important forâ |δm 2 31 |, or equivalently, β 0.

Diagonalization of a 2 × 2 hermitian matrix
The Jacobi method entails diagonalizing 2 × 2 submatrices repeatedly. For this, it is convenient to note that given a 2 × 2 hermitian matrix, the unitary matrix which diagonalizes this is given by where c ω = cos ω , s ω = sin ω ,  That is, the double signs corresponding to the two possible quadrants for 2ω that satisfy eq. (3.12). In applying the above formula to our problem, care needs to be taken to choose the correct quadrant and sign combination so that the resulting effective mixing angles and masssquared eigenvalues run smoothly from their vacuum values.

η = 0 case, first and second rotations
As mentioned above, we will first approximately diagonalize H 0 , and add on theâηM η term afterwards. Here, we reproduce how the Jacobi method was applied to H 0 in refs. [68,70]. There, a (1, 2) rotation was applied to H 0 , followed by a (2, 3) rotation, which was sufficient to approximately diagonalize H 0 .
Using W , we obtain  where the upper signs are for the δm 2 31 > 0 (normal hierarchy) case and the lower signs are for the δm 2 31 < 0 (inverted hierarchy) case, with As β is increased beyond 0, the λ ± asymptote to for both mass hierarchies. Note that λ − < 0 for the δm 2 31 < 0 case. The β-dependences of λ ± are shown in figure 4. Order-of-magnitude-wise, we have (3.29) In particular, in the range β 0 we have Thus, looking at the sizes of the elements of H 0 in that range we find: where the elements with two entries denote the two different mass hierarchies, O(NH/IH), and we can see that further diagonalization only require angles of order O( 3 ). Therefore, H 0 can be considered approximately diagonal.
The β-dependence of χ is shown in figure 6 for several values of η, both positive (figure 6(a)) and negative (figure 6(b)). Using X, we find The β-dependence of λ X± are shown in figures 6(c) to 6(f).

Effective mixing angles for neutrinos
We have discovered that the unitary matrix which approximately diagonalizes H η is the PMNS matrix U in vacuum can be parametrized as In the following, we rewrite the mixing matrix in matter ∼ U into the analogous form where Q 3 was defined in eq. (3.4).
1. δm 2 31 > 0 case. Using eq. (3.54), it is straightforward to show that where in the last and penultimate lines we have combined the two 12-rotations into one. We now commute R 23 (φ, 0)R 12 (χ, δ) through the other mixing matrices to the left as follows: In the range β 0, the angle θ 12 is approximately equal to π/2 so we can approximate for any φ. On the other hand, in the range β −1 the angle φ is negligibly small so we can approximate for any θ 12 . Therefore, for all β we have In the range β 0, the angle θ 12 is approximately equal to π/2 as we have noted above and we have the approximation given in eq. (3.56). Note that for any χ. On the other hand, in the range β 0, the angle χ is negligibly small so we can approximate for any θ 12 . Therefore, for all β we see that Step 3: commutation of R 12 (χ, δ) through R 13 (θ 13 , 0).
When δm 2 31 > 0 we have θ 13 ≈ π 2 in the range β 1 so we can approximate for any χ. On the other hand, in the range β 0 the angle χ was negligibly small so that we had eq. (3.63). Note that for any θ 13 . Therefore, for all β we see that and using eq. (3.54) we obtain where in the last and penultimate lines we have combined the two 23-rotations into one. The matrix Q 3 on the far right can be absorbed into the redefinitions of Majorana phases and can be dropped.
Thus, we find that the effective mixing matrix ∼ U in the case δm 2 31 > 0 can be expressed as eq. (3.53) with the effective mixing angles and effective CP-violating phase given approximately by (3.72) 2. δm 2 31 < 0 case. Using eq. (3.54), we obtain δ. The first step is the same as the δm 2 31 > 0 case, the only difference being the β-dependence of θ 13 , which is also shown in figure 3 In the range β 0 the angle θ 12 is approximately equal to π/2 as we have noted previously, and we have the approximation given in eq. (3.56). Note that for any ψ. On the other hand, in the range β 0 the angle ψ is negligibly small so that Therefore, for all β we see that and When δm 2 31 < 0 we have θ 13 ≈ 0 in the range β 1 so that for all ψ. On the other hand, in the range β 0 the angle ψ was negligibly small so that we had the approximation eq. (3.75). Note that for all θ 13 . Therefore, for all β we see that and using eq. (3.54) we obtain where in the last and penultimate lines we have combined the two 23-rotations into one. The matrix Q 3 on the far right can be absorbed into redefinitions of the Majorana phases and can be dropped. Thus, we find that the effective mixing matrix ∼ U in the case δm 2 31 < 0 can be expressed as eq. (3.53) with the effective mixing angles and effective CP-violating phase given approximately by (3.84)

Summary of neutrino case
To summarize what we have learned, inclusion of theâηM η term in the effective Hamiltonian shifts θ 23 to θ 23 = θ 23 + χ for the δm 2 31 > 0 case, and to θ 23 = θ 23 + ψ for the δm 2 31 < 0 case. For both cases, θ 23 can be calculated directly without calculating χ or ψ first via the expression Note that as β is increased, θ 23 runs toward π 2 if δm 2 31 η > 0, while it runs toward 0 if δm 2 31 η < 0. The β-dependence of θ 23 is shown in figure 8. The CP-violating phase δ is unaltered and maintains its vacuum value.
The running of the effective mass-squared differences are also modified in the range β 0. For the δm 2 31 > 0 case, λ 1 and λ 2 show extra running, while for the δm 2 31 < 0 case, it is λ 1 and λ 3 that show extra running, cf. figures 6 and 7.

Discussion at the probability level
So far we have focused our attention on how the flavor-diagonal NSI parameter η modifies the running of the effective mass-squared differences, mixing angles, and CP-violating phase as functions ofâ = a(1 + ζ) for the neutrinos. We derived simple analytical expressions -26 -

JHEP11(2015)035
for these effective parameters using the Jacobi method. We have discussed the running of these effective neutrino oscillation parameters for both the normal (δm 2 31 > 0) and inverted (δm 2 31 < 0) neutrino mass hierarchies. The modifications induced by η and ζ in the running of effective oscillation parameters in the case of anti-neutrino are discussed in detail in appendix A. At this point, we look at how these lepton-flavor-conserving NSI parameters alter the neutrino oscillation probabilities for various appearance and disappearance channels.
In the three-flavor scenario, the neutrino oscillation probabilities in vacuum 5 for the disappearance channel (initial and final flavors are same) take the form and for the appearance channel (initial and final flavors are different) we have In the above equations, we define The transition probabilities in eq. (3.86) and eq. (3.87) contain several elements of the PMNS matrix U which are expressed in terms of three mixing angles θ 12 , θ 23 , θ 13 , and a CP-violating phase δ as shown in eq. (3.52). The oscillation probabilities for the antineutrinos are obtained by replacing U αi with its complex conjugate. The neutrino oscillation probabilities in the presence of matter are obtained by replacing the vacuum expressions of the elements of the mixing matrix U and the mass-square differences ∆ ij with their effective 'running' values in matter [21][22][23]: and for the anti-neutrinos To demonstrate the accuracy (or lack thereof in special cases) of our approximate analytical results, we compare the oscillation probabilities calculated with our approximate effective 5 We follow the conventions, notations, and the expressions for various neutrino oscillation probabilities in vacuum as given in the appendices A.1 and A.2 of ref. [70]. running mixing angles and mass-squared differences with those calculated numerically for the same baseline and line-averaged constant matter density along it. For the mixing angles and mass-squared differences in vacuum, we use the benchmark values from ref. [194] as listed in table 1. In some plots, we take different values of sin 2 θ 23 and δ which we mention explicitly in the figure legends and captions. In this paper, all the plots (except in appendix B) are generated considering the line-averaged constant matter density for a given baseline which has been estimated using the Preliminary Reference Earth Model (PREM) [195]. In appendix B, we compare the exact numerical probabilities with lineaveraged constant Earth density and varying Earth density profile for 8770 km and 10000 km baselines.
In figure 9, we present our approximate ν µ → ν e oscillation probabilities (blue curves) as a function of the neutrino energy against the exact numerical results (red curves) considering 6 η = 0.1, ζ = 0 (left panels) and η = −0.1, ζ = 0 (right panels). The upper panels are 6 We take ζ = 0 in our plots since we expect it to be hidden in the uncertainties in the matter density and neutrino energy. drawn for the baseline of L = 2300 km, which corresponds to the distance between CERN and Pyhäsalmi [196][197][198] with the line-averaged constant matter density of ρ = 3.54 g/cm 3 . In the lower panels, we give the probabilities for the baseline of L = 8770 km, which is the distance from CERN to Kamioka [199] assuming ρ = 4.33 g/cm 3 . Here, in all the panels, we assume sin 2 θ 23 = 0.41 (θ 23 = 40 • ), δ = 0 • , and normal mass hierarchy (δm 2 31 > 0). To see the differences in the oscillation probability caused by the NSI parameters, we also give the exact numerical SM three-flavor oscillation probabilities in matter in the absence of NSI's which are depicted by the solid black curves with figure legend 'SM, Exact.' It has been already shown in ref. [70] that our approximate expressions for the η = 0 case match extremely well with the exact numerical results for all these baselines and energies. We also compare our results with the approximate expressions of Asano and Minakata 7 [69] (dashed green curves). The correspondence between our η and ζ and the NSI parameters ε ee , ε µµ , and ε τ τ used in refs. [69,201] can be obtained via eqs. (2.4), (2.7), and (2.9) which suggest the changes: a →â ≡ a(1 + ζ), ε ee → 0, ε µµ → η, and ε τ τ → −η. We can see from figure 9 that for the 2300 km baseline, the Asano-Minakata expressions give better matches compared to our results, while for the 8770 km baseline, our expressions agree better with the exact numerical results.
The accuracy of our analytical approximations as compared to the exact numerical results for different vacuum values of θ 23 is demonstrated in figure 10. Here we consider the minimum (35 • ) and maximum (55 • ) values of θ 23 which are allowed in the 3σ range [194]. We also present the results for the maximal mixing choice. All the plots in figure 10 have been generated assuming δ = 0 • and normal mass hierarchy (δm 2 31 > 0). We consider the same choices of η and ζ as in figure 9 and results are given for 2300 km (upper panels) and 8770 km (lower panels) baselines. As is evident, our approximation provides satisfactory match with exact numerical results for different values of θ 23 . Figure 11 compares our approximate probability expressions (solid curves) against the exact numerical results (dashed curves) assuming four different values of the CP-violating phase δ (0, 90, 180, and 270 degrees) at 2300 km. Here, we consider θ 23 = 40 • and δm 2 31 > 0. These plots clearly show that our approximate expressions work quite well even for nonzero δ and can predict almost accurate L/E patterns of the oscillation probability for finite δ and η. It also suggests that one can explain qualitatively the possible correlations and degeneracies between δ and η using these analytical expressions which cannot be tackled with numerical studies.
In figure 12, we plot the ν µ → ν µ survival probability in the presence of NSI for two different vacuum values of θ 23 (40 • and 50 • ) at 8770 km. We show the matching between the analytical and numerical results assuming η = 0.1, ζ = 0 (left panel) and η = −0.1, ζ = 0 (right panel). We also give the exact numerical standard three-flavor oscillation probabilities in matter without NSI's so that one can compare them with the finite η case. In both the panels, we assume δ = 0 • and δm 2 31 > 0. Figure 12 shows that our approximate expressions match quite nicely with the numerical results. Note that at  higher energies, the impact of NSI's are quite large in the ν µ → ν µ survival channel and there is a substantial difference in the standard and NSI probabilities for both the choices of θ 23 which can be probed in future long-baseline [202,203] and atmospheric [204][205][206] neutrino oscillation experiments.  Another important point to be noted that in the absence of NSI, the standard probability curves for both the values of θ 23 almost overlap with each other at higher energies, whereas with NSI, there is a large separation between them. It immediately suggests that the corrections in the probability expressions due to the NSI terms depend significantly on whether the vacuum value of θ 23 lies below or above 45 • [31]. Figure 12 also indicates that there are degeneracies between the octant of θ 23 and the sign of NSI parameter η for a given choice of hierarchy. For an example, P µµ (θ 23 = 50 • , η = 0.1) in the left panel is almost same with P µµ (θ 23 = 40 • , η = −0.1) in the right panel for the energies above 12 GeV or so. Again, P µµ (θ 23 = 40 • , η = 0.1) in the left panel matches quite well with P µµ (θ 23 = 50 • , η = −0.1) in the right panel. These kinds of degeneracies can be well explained qualitatively with the help of our analytical expressions. We discuss this issue in detail in the next section which is one of the highlights of this work.

Possible applications of analytical expressions
In this section, we discuss the utility of our analytical probability expressions to determine the conditions for which the impact of the NSI parameter η becomes significant. We also give simple and compact analytical expressions to show the possible correlations and degeneracies between θ 23 and η under such situations. We begin our discussion with electron neutrinos.

ν e → ν α oscillation channels
Let us first consider the ν e → ν α (α = e, µ, τ ) oscillation channels in matter in the presence of the NSI parameter η. Since we expect the effect of η to become important in the range β 0, we set s 12 ≈ 1, c 12 ≈ 0 (which is valid in the range β  leads to the following simple expressions: Recall that when δm 2 31 > 0, the effective mixing angle θ 13 increases monotonically toward π/2, going through π/4 around β ∼ 0, while in the δm 2 31 < 0 case, it decreases monotonically toward 0, cf. Figure 3(b). This will cause sin 2 (2θ 13 ) to peak prominently around β ∼ 0 for the normal hierarchy case, but not for the inverted hierarchy case as shown in figure 13(a). As discussed in ref. [70], demanding that sin 2 ( ∼ ∆ 32 /2) also peaks at the same energy leads to the requirements of L ∼ 10000 km and E ∼ 7 GeV. Thus, measuring ∼ P (ν e → ν e ) survival probability at this baseline and energy will allow us to discriminate between the normal and inverted mass hierarchies irrespective of the value of θ 23 [207].
Also, around β ∼ 0, ∼ ∆ 32 is not affected by η provided η 0.1 (see middle and bottom panels of figure 6), allowing this channel to determine the mass hierarchy free from any NSI effect.
To see the effect of η, we need to observe the running of θ 23 . This could be visible in the ν e → ν µ and ν e → ν τ appearance channels for the normal hierarchy case as a change in the heights of the oscillation peaks around E ∼ 7 GeV provided θ 23 deviates sufficiently from the vacuum value of θ 23 at that energy. The running of s 2 23 for various values of η is depicted in figure 13(b). Atâ ∼ δm 2 31 , eq. (3.85) tells us that sible correlations and degeneracies between θ 23 and η. Such correlations and degeneracies could also be found numerically, but the reason for those features will not be so transparent. Note that eq. (4.4) suggests that a positive value of η would enhance ∼ P (ν e → ν µ ) while suppressing ∼ P (ν e → ν τ ), and a negative η would do the opposite. Figure 14 confirms this feature where we plot the ν e → ν µ (ν e → ν τ ) transition probability in the left (right) panel. In both panels, the standard oscillation probabilities in matter without NSI's for four different values of CP phase 8 δ are given by the solid red lines. Blue lines (black dashed lines) depict the approximate analytical (exact numerical) probabilities with η = 0.1 and four different values of δ. Figure 14 infers that the effect of η is visible in these channels provided the vacuum value of θ 23 is sufficiently well known and η is also large enough.

ν µ → ν α oscillation channels
Next, let us consider ν µ → ν α (α = µ, τ ) oscillation channels. In addition to assuming β −2, which allows us to set s 12 ≈ 1, c 12 ≈ 0, we further restrict our consideration to the range β 1, which allows us to set s 13 ≈ 1, c 13 ≈ 0 or s 13 ≈ 0, c 13 ≈ 1 depending on whether δm 2 31 > 0 or δm 2 31 < 0 (see figure 3(b)). With these conditions, we obtain the following simple expressions:  Figure 14 shows that the impact of the CP phase δ is quite weak around β ∼ 0. In this region, θ 12 approaches to π/2 so that s 12 ≈ 1 and c 12 ≈ 0. Therefore, the Jarlskog Invariant [208] in matter J = s 12 c 12 s 13 c 2 13 s 23 c 23 sin δ almost approaches to zero diminishing the effect of δ. This argument works even in the presence of the NSI parameter η. Our simple and compact approximate probability expressions given by eq. (4.2) and eq. (4.3) also validate this point as there are no δ-dependent terms in these expressions. .

(4.8)
In the absence of η, θ 23 does not run and sin 2 (2θ 23 ) will maintain its vacuum value close to one. In the presence of a non-zero η, however, θ 23 will run towards either π 2 or 0 depending on the sign of δm 2 31 η, as was shown in figure 8, and sin 2 (2θ 23 ) will run toward zero in both cases. This is depicted in figure 15.
Let us see how the ∆ factors in eq. (4.6) and eq. (4.8) behave in the range β 1. For the normal hierarchy case, we have while for the inverted hierarchy case, they take the form Since the sign of δλ ij does not affect the value of sin 2 ( ∼ ∆ ij /2), both mass hierarchies lead to the same asymptotic oscillation probabilities. We will therefore only consider the normal Therefore, for fixed baseline L and matter density ρ, as the neutrino energy E is increased, ∼ ∆ 21 /2 damps to zero when η = 0, but asymptotes to a constant value proportional to |η| when η = 0. This is demonstrated in figure 16 for the baseline L = 10000 km with average matter density ρ = 4.53 g/cm 3 . Consequently, when η = 0, the factor sin 2 (2θ 23 ) stays constant while sin 2 ( ∼ ∆ 21 /2) damps to zero as we increase E, while in the η = 0 case, the factor sin 2 (2θ 23 ) damps to zero while sin 2 ( ∼ ∆ 21 /2) asymptotes to a constant value as E is increased. In either case, the ν µ → ν τ oscillation probability is suppressed at high energy.
The difference between the η = 0 and η = 0 cases could manifest itself at the ν µ → ν τ oscillation peak which happens at η + · · · (4.15) which at β ≈ 0.6 yields θ 23 ≈ θ 23 ± 3η . (4.16) In the above equation, upper (lower) sign is applicable for the normal (inverted) hierarchy case. It suggests that there is a degeneracy between the choices of sign of δm 2 31 and the sign of η which give rise to same amount of corrections in θ 23 . To observe the shift in the oscillation probability, we keep terms up to order η 2 since the linear term is suppressed due to the fact that θ 23 is close to π/4. Therefore, we have sin(2θ 23 ) = sin(2θ 23 ) + 2 cos(2θ 23 )δθ 23 − 2 sin(2θ 23 )(δθ 23 ) 2 + · · · ≈ sin(2θ 23 ) ± 6 cos(2θ 23 )η − 18 sin(2θ 23 )η 2 . (4.17) For the benchmark value of s 2 23 = 0.4 (0.6) in the lower (higher) octant, above equation can be written as where upper (lower) signs are for the normal (inverted) hierarchy. The above equations clearly reveal that for a given choice of hierarchy, there are degenerate solutions 9 of the octant of θ 23 and the sign of NSI parameter η (s 2 23 = 0.4, η = ±0.1 and s 2 23 = 0.6, η = ∓0.1), giving rise to same value of sin(2θ 23 ). Figure 17 shows the variation in θ 23 and sin 2 (2θ 23 ) 9 We have already seen this degeneracy in the P (νµ → νµ) oscillation channel in figure 12.  figure 17 where we consider s 2 23 = 0.6. All these observations in figure 17 and figure 18 suggest that our approximate calculations are valid.

Summary and conclusions
Analytical studies of the neutrino oscillation probabilities are inevitable to understand how neutrino interactions with matter modify the mixing angles and mass-squared differences in a complicated manner in a three-flavor framework. In previous papers [68,70], we showed that the neutrino oscillation probabilities in matter can be well understood if we allow the mixing angles and mass-squared differences in the standard parametrization to 'run' with the matter effect parameter a = 2 √ 2G F N e E, where N e is the electron density in matter and E is the neutrino energy. We managed to derive simple and compact analytical approximations to these running parameters using the Jacobi method. We found that for large θ 13 , the entire matter effect could be absorbed into the running of the effective mass-squared differences and the effective mixing angles θ 12 and θ 13 , while neglecting the running of the mixing angle θ 23 and the CP-violating phase δ.

JHEP11(2015)035
In this paper, we extended our analysis to study how the running of the neutrino oscillation parameters in matter would be altered in the presence of NSI's of neutrinos with the matter fermions. Such NSI's are predicted in most of the new physics models that attempt to explain the non-zero neutrino masses, as well as in a wealth of various other BSM models. There, the NSI's are simply the effective four-fermion interactions at the energy scales relevant for neutrino oscillation experiments that remain when the heavy mediator fields of the full theory are integrated out. These NSI's give rise to new neutralcurrent type interactions of neutrinos (both flavor-conserving and flavor-violating) during their propagation through matter on top of the SM interactions, causing the change in the effective mass matrix for the neutrinos which ultimately affect the running of the oscillation parameters and hence change the oscillation probabilities between different neutrino flavors. These sub-leading new physics effects in the probability due to NSI's can be probed in upcoming long-baseline and atmospheric neutrino oscillation experiments.
In this work, we restricted our attention to the matter effect of flavor-conserving, nonuniversal NSI's of the neutrino, relegating the discussion of the flavor-violating NSI case to a separate paper [73]. The relevant linear combinations of the flavor-diagonal NSI's were η = (ε µµ − ε τ τ )/2 and ζ = ε ee − (ε µµ + ε τ τ )/2, where a non-zero ζ led to a rescaling of the matter-effect parameter a →â = a(1 + ζ), while a non-zero η led to non-trivial modifications on how the running oscillation parameters depend onâ.
Utilizing the Jacobi method, as in refs. [68,70], we obtained approximate analytical expressions for the effective neutrino oscillation parameters to study how they 'run' with the rescaled matter-effect parameterâ, and to explore the role of non-zero η in neutrino oscillation. We found that in addition to the two rotations, which were required for the SM matter interaction and were absorbed into effective values of θ 12 and θ 13 , a third rotation was needed to capture the effects of η, which could be absorbed into the effective value of θ 23 . Thus, within the neutrino mixing matrix, the effect of η appears as a shift in the effective mixing angle θ 23 , while the SM matter effects show up as shifts in θ 12 and θ 13 . The CP-violating phase δ remains unaffected and maintains its vacuum value. The running of all the effective neutrino oscillation parameters were presented for both the normal and inverted neutrino mass hierarchies. The changes caused by η in the running of the effective oscillation parameters for the anti-neutrino case are discussed in detail in appendix A.
We have also studied the impact of the lepton-flavor-conserving NSI parameters on the neutrino oscillation probabilities for various appearance and the disappearance channels. To demonstrate the accuracy (or lack thereof in special cases) of our approximate analytical expressions, we compared the oscillation probabilities estimated with our approximate effective 'running' mixing angles and mass-squared differences with those calculated numerically for the same choices of benchmark oscillation parameters, energy, baseline, and line-averaged constant matter density along it. We found that our approximation provided satisfactory matches with exact numerical results in light of large θ 13 for different values of θ 23 , CP-violating phase δ, and for positive and negative values of the NSI parameter η. A comparison of our results with the approximate expressions of Asano and Minakata [69] for the ν µ → ν e appearance channel has also been presented.
Finally, we examined the merit of our analytical probability expressions to identify the situations at which the impact of the NSI's become compelling. It was found that at higher -38 -

JHEP11(2015)035
baselines and energies, the impact of η can be quite significant in the ν µ → ν µ survival channel if |η| is of the order of its current experimental upper bound. A considerable difference between the SM and NSI probabilities can be seen irrespective of the vacuum value of θ 23 , and the sign of the NSI parameter η. We note that this feature may be explorable with the upcoming 50 kiloton magnetized iron calorimeter detector at the Indiabased Neutrino Observatory (INO), which aims to detect atmospheric neutrinos and antineutrinos separately over a wide range of energies and path lengths [209]. Using our analytical approach, we showed in a very simple and compact fashion that the corrections in θ 23 due to the η depend significantly on whether the vacuum value of θ 23 lies below or above 45 • , suggesting a possible degeneracy between the octant of θ 23 and the sign of η for a given choice of mass hierarchy.

Acknowledgments
We would like to thank Minako Honda and Naotoshi Okamura for their contributions to ref. [82], the predecessor to this work. A Effective mixing angles and effective mass-squared differences -Antineutrino case In this appendix, we study the matter effect due to the anti-neutrino NSI's. We again utilize the Jacobi method to estimate how the NSI parameter η alters the 'running' of the effective mixing angles, effective mass-squared differences, and the effective CP-violating phase δ in matter for the anti-neutrinos. Like the neutrino case, we also present here a comparison between our approximate analytical probability expressions and exact numerical calculations towards the end of this appendix.

A.1 Differences from the neutrino case
For the anti-neutrinos, the effective Hamiltonian is given by

JHEP11(2015)035
The differences from the neutrino case are the reversal of signs of the CP-violating phase δ (and thus the complex conjugation of the PMNS matrix U ), and the matter interaction parameterâ = a(1 + ζ). We denote the matter effect corrected diagonalization matrix as U (note the mirror image tilde on top) to distinguish it from that for the neutrinos.

A.2.1 Change to the mass eigenbasis in vacuum
Using the matrix Q 3 from eq. (3.4), we begin by partially diagonalize the effective Hamiltonian H η as where M a (θ 12 , θ 13 , θ 23  and As in the neutrino case, we will first approximately diagonalize H 0 and then add on thê aηM η term later. The Jacobi method applied to H 0 is as follows.

13
, 0 <θ 12 ≤ θ 12 . (A.10) The dependences ofθ 12 and λ ± on β are plotted in figure 19. Note that in contrast to the neutrino case,θ 12 decreases monotonically from θ 12 to zero as β is increased. The β-dependences of s 12 = sin θ 12 and c 12 = cos θ 12 are shown in figure 20(a). As β is increased beyond β = −2, c 12 grows quickly to one while s 12 damps quickly to zero. The productâs 12 stops increasing around β = −2 and plateau's to the asymptotic value of δm 2 21 s 12 c 12 /c 2 13 = |δm 2 31 |O( 2 ) as shown in figure 20(b). That is: Note also that the scales of λ ± are given simply by since no level crossing occurs in this case. 2. Second rotation. Sinceâc 12 continues to increase with β whileâs 12 does not, we perform a (1, 3) rotation on H 0 next.
1. δm 2 31 > 0 case. In the δm 2 31 > 0 caseâs 13 → O( )|δm 2 31 | as β is increased beyond 0. Therefore, we can approximate  Note that The β-dependence of χ is shown in figure 24 for several values of η, both positive (figure 24(a)) and negative ( figure 24(b)). Using X, we find Note that λ X± have the same asymptotics as λ X± for the neutrino case except with the η > 0 and η < 0 cases reversed. The β-dependence of λ X± are shown in figures 24(c) to (f).
2. δm 2 31 < 0 case. In the δm 2 31 < 0 caseâc 13 → O( )|δm 2 31 | as β is increased beyond 0. Therefore, we can approximate   Note that The β-dependence of ψ is shown in figure 25 for several values of η, both positive (figure 25(a)) and negative ( figure 25(b)). Using Y , we find Note that λ Y ± have the same asymptotics as λ Y ± for the neutrino case except with the η > 0 and η < 0 cases reversed. The β-dependence of ψ and λ Y ± are shown in figure 25.

A.3 Effective mixing angles for anti-neutrinos
We have discovered that the unitary matrix which approximately diagonalizes H is U * = U * Q * 3 V W X when δm 2 31 > 0, and U * = U * Q * 3 V W Y when δm 2 31 < 0. Taking the complex conjugate, these are respectively U = U Q 3 V W X * when δm 2 31 > 0, and U = U Q 3 V W Y * when δm 2 31 < 0. In the following, we rewrite the mixing matrix in matter U into the form absorbing the extra mixing angles and CP phase into appropriate definitions of the 'running' parameters θ 12 , θ 13 , θ 23 , and δ. As in the neutrino case, frequent use is made of eq. (3.54).
In the range β 0, the angle θ 12 is approximately equal to zero as we have noted above and we have the approximation given in eq. (A.42). Note that for any χ. On the other hand, in the range β 0, the angle χ is negligibly small so we can approximate for any θ 12 . Therefore, for all β we see that When δm 2 31 > 0 we have θ 13 ≈ 0 in the range β 1 so we can approximate for any χ. On the other hand, in the range β 0 the angle χ was negligibly small so that we had eq. (A.49). Note that for any θ 13 . Therefore, for all β we see that and using eq. (3.54) we obtain where in the last and penultimate lines we have combined the two 23-rotations into one. The matrix Q 3 on the far right can be absorbed into the redefinitions of Majorana phases and can be dropped.

JHEP11(2015)035
We now commute R 13 (φ, 0)R 12 (ψ, δ) through the other mixing matrices to the left and re-express U as in eq. (A. 40), absorbing the extra mixing angles and CP phase into θ 12 , θ 13 , θ 23 , and δ. The first step is the same as the δm 2 31 > 0 case, the only difference being the β-dependence of θ 13 , which is also shown in figure 21(b).
When δm 2 31 < 0 we have θ 13 ≈ π 2 in the range β 1 so that for all ψ. On the other hand, in the range β 0 the angle ψ was negligibly small so that we had the approximation eq. (3.75). Note that for all θ 13 . Therefore, for all β we see that where in the last and penultimate lines we have combined the two 23-rotations into one. The matrix Q 3 on the far right can be absorbed into redefinitions of the Majorana phases and can be dropped.
The running of the effective mass-squared differences are also modified in the range β 0. For the δm 2 31 > 0 case, λ 2 and λ 3 show extra running, while for the δm 2 31 < 0 case, it is λ 1 and λ 2 that show extra running, cf. figures 24 and 25.

A.5 Discussion at the probability level
Now, we demonstrate how these lepton-flavor-conserving NSI parameters affect the antineutrino oscillation probability for the various appearance and the disappearance channels. In figure 27, we compare our approximateν µ →ν e oscillation probabilities (blue curves) as a function of the neutrino energy against the exact numerical results (red curves) assuming η = 0.1, ζ = 0 (left panels) and η = −0.1, ζ = 0 (right panels). The upper (lower) panels are drawn for 2300 km (8770 km) baseline. Here, in all the panels, we assume θ 23   To generate these plots, we consider θ 23 = 40 • and δm 2 31 < 0. δ = 0 • , and inverted mass hierarchy (δm 2 31 < 0). We also compare our results with the approximate expressions of Asano and Minakata [69] (dashed green curves).
The accuracy of our analytical approximations as compared to the exact numerical results for different values of θ 23 is shown in figure 28. All the plots in figure 28   generated assuming δ = 0 • and inverted mass hierarchy (δm 2 31 < 0). We take the same choices of η and ζ like in figure 27 and results are given for 2300 km (upper panels) and 8770 km (lower panels) baselines. We find that our approximation provides satisfactory match with exact numerical results for different values of θ 23 . Figure 29 presents a comparison of our approximate probability expressions (solid curves) against the exact numerical results (dashed curves) assuming four different values of the CP-violating phase δ at 2300 km. Here, we take θ 23 = 40 • and δm 2 31 < 0. In figure 30, we depict theν µ →ν µ survival probability in the presence of NSI for two different values of θ 23 (40 • and 50 • ) at 8770 km baseline. We present the matching between the analytical and numerical results assuming η = 0.1, ζ = 0 (left panel) and η = −0.1, ζ = 0 (right panel). In both the panels, we consider δ = 0 • and δm 2 31 < 0. Figure 30 portrays that our approximate expressions match quite nicely with the numerical results.

B Comparing probabilities with constant & varying Earth density profile
So far, we considered the line-averaged constant Earth matter density for a given baseline which has been estimated using the PREM profile to present our results. Now, it would be quite interesting to study how the exact numerical probabilities would be affected if we consider the more realistic varying Earth density profile instead of the line-averaged constant matter density for the baselines as large as 8770 km and 10000 km. In figure 31, we show the exact numerical probabilities considering the constant and varying Earth density profiles for 8700 km (upper panels) and 10000 km (lower panels) baselines. In the figure legends, the line-averaged constant matter density cases are denoted by 'ρ avg ' and the varying Earth density cases are labelled by 'PREM'. We perform these comparisons for both the SM and NSI scenarios. For the NSI's, we take η = 0.1, ζ = 0 (left panels), and . ν µ → ν e transition probability as a function of neutrino energy E in GeV for 8770 km (10000 km) baseline in upper (lower) panels. We compare the exact numerical probabilities considering the constant (denoted by 'ρ avg '), and varying (labelled by 'PREM') Earth density profiles for both the SM and NSI cases. We also plot our approximate analytical expressions in the presence of NSI. For the NSI's, we take η = 0.1, ζ = 0 (left panels), and η = −0.1, ζ = 0 (right panels). In all the panels, we consider θ 23 = 40 • , δ = 0 • , and normal mass hierarchy. η = −0.1, ζ = 0 (right panels). For the sake of completeness, we also plot the approximate probability expressions in the presence of NSI that we derived in this paper assuming the line-averaged constant Earth matter density based on the PREM profile. Figure 31 clearly shows that though the line-averaged constant matter density probability does not completely overlap with the PREM-based varying matter density profile probability, it is still fairly accurate. Moreover, the effect due to the inclusion of flavor-diagonal NSI's is correctly captured in our approximate analytical expressions.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.