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 $\theta_{12}$ and $\theta_{13}$, while the flavor-diagonal NSI's only affect $\theta_{23}$. The CP-violating phase $\delta$ 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 numu to numu 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 $\theta_{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 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].

non-constant mass-density, if the need arises.
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 neutrino mixing matrix, the SM interaction only affects the running of θ 12 and θ 13 , while the flavor-diagonal NSI's only affect the running of θ 23 . The CP-violating phase δ remains unaffected and maintains its vacuum value.
We note that similar studies have been performed in the past by many authors e.g. in Refs. [69,[74][75][76][77][78][79][80][81]. This work differs from these existing works in the use of the Jacobi method [72] to derive compact analytical formulae for the running effective mass-squared differences and effective mixing angles, which provide a clear and simple picture of how neutrino NSI's affect neutrino oscillation. We also note that we have addressed the same problem with a similar approach previously in Refs. [82][83][84][85]. The current paper updates these works by allowing for a non-maximal value of θ 23 , a generic value of the CP-violating phase δ, a larger range of a, and refinements on how the matter effect is absorbed into the running parameters.
This paper is organized as follows. We start section 2 with a discussion on neutrino NSI's: how and where they arise, how they affect the propagation of the neutrinos in matter, and show that the linear combinations relevant for neutrino oscillation are η = (ε µµ −ε τ τ )/2 and ζ = ε ee − (ε µµ + ε τ τ )/2. This is followed by a brief discussion on the theoretical expectation for the sizes of these parameters, and their current experimental bounds. In section 3, we use the Jacobi method to calculate how the NSI parameter η affects the running of the effective mass-squared differences, effective mixing angles, and the effective CP-violating phase as functions ofâ = a(1 + ζ) for the neutrinos. To check the accuracy of our method, we also present a comparison between our approximate analytical probability expressions and exact numerical calculations (for constant matter density) towards the end of this section. Section 4 describes the advantages of our analytical probability expressions to figure out the suitable testbeds to probe these NSI's of the neutrino. In this section, we also present simple and compact analytical expressions exposing the possible correlations and degeneracies between θ 23 and the NSI parameter η under such situations. Finally, we summarize and draw our conclusions in section 5. The derivation of the running oscillation parameters for the anti-neutrino case is relegated to appendix A where we also compare our analytical results with exact numerical probabilities. In appendix B, we examine the differences in the exact numerical probabilities with line-averaged constant Earth density and varying Earth density profile for 8770 km and 10000 km baselines.

A Brief Tour of Non-Standard Interactions of the Neutrino
In this section, we first briefly discuss the various categories of models that give rise to NSI's of neutrino.
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 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.

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 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 η. 4 The notation used in Ref. [174] is ε = ε d µτ and ε = ε d τ τ − ε d µµ .  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.

Parameter
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

6) 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, where 14) 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 Fig. 4. Order-of-magnitude-wise, we have (3.29) In particular, in the range β 0 we have For the off-diagonal terms, sinceâc 12 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.
Using X, we find The β-dependence of λ X± are shown in Figs. 6(c) to 6(f).
Comparing Eq. (3.39) and Eq. (3.46), we can infer that ψ(η) ≈ χ(−η), the small difference due to the δm 2 21 term in the denominator of the expressions for tan 2χ and tan 2ψ. This can be seen in Fig. 7 where the β-dependence of ψ is shown for several values of η, both positive ( Fig. 7(a)) and negative ( Fig. 7(b)).
Using Y , we find (3.49) Thus, H η− is approximately diagonal. The asymptotic forms of λ Y ± at β 0 are The β-dependence of λ Y ± are shown in Fig. 7(c) to Fig. 7(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).
• δm 2 31 > 0 Case: 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 Note that for any θ 12 . Therefore, for all β we have and 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 Note that for any θ 12 . Therefore, for all β we see that and 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.
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) • δm 2 31 < 0 Case: δ. 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 Fig. 3(b).
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 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.
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)
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. Figs. 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 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 CPviolating phase δ as shown in Eq. (3.52). The oscillation probabilities for the anti-neutrinos 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 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 Fig. 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 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] 6 We take ζ = 0 in our plots since we expect it to be hidden in the uncertainties in the matter density and neutrino energy. 7 For comparison, we take Eq. (36) of Ref. [69] where the authors adopted the perturbation method [80,200] to obtain the analytical expressions for oscillation probability in presence of NSI's for large θ13.
The accuracy of our analytical approximations as compared to the exact numerical results for different vacuum values of θ 23 is demonstrated in Fig. 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 Fig. 10 have been generated assuming δ = 0 • and normal mass hierarchy (δm 2 31 > 0). We consider the same choices of η and ζ as in Fig. 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 .   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 Fig. 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. Fig. 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].

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 β −2, see Fig. 2(a)), which 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. Fig. 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 Fig. 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 Fig. 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 Fig. 13(b). Atâ ∼ δm 2 31 , Eq. (3.85) tells us that showing the possible shift in θ 23 due to η in a simple and compact fashion which clearly establishes the merit of our analytical approximation. Eq. (4.4) also shows compactly possible 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. Fig. 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 δ. Fig. 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.

(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 Fig. 8, and sin 2 (2θ 23 ) will run toward zero in both cases. This is depicted in Fig. 15.
8 Fig. 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.
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 hierarchy case in the following. Recalling (see Eq. (2.2)) that we find .

(4.13)
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 Fig. 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.

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 δ.
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 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.

A Effective Mixing Angles and Effective Mass-Squared Differences -Anti-Neutrino 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 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:
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 Fig. 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 Fig. 20(b). That is:â Note also that the scales of λ ± are given simply by since no level crossing occurs in this case.

Second Rotation
Sinceâc 12 continues to increase with β whileâs 12 does not, we perform a (1, 3) rotation on H 0 next.
Define the matrix W as  . (A.14) The angle φ is in the fourth quadrant when δm 2 31 > 0 (normal hierarchy), and in the first quadrant when δm 2 31 < 0 (inverted hierarchy). The β-dependence of φ is shown in Fig. 21(a) for both mass hierarchies.
Using W , we obtain for both mass hierarchies. Note that λ − < 0 < λ + for the δm 2 31 > 0 case, while both λ ± < 0 for the δm 2 31 < 0 case. The β-dependences of λ ± are shown in Fig. 22. Order-of-magnitude-wise, we have (A. 18) In particular, in the range β 0, we have For the off-diagonal terms, we haveâs 12 = |δm 2 31 |O( 2 ), c 13 = O(1), s 13 = O( ), and Thus, looking at 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 see that further diagonalization require angle of order O( 3 ). Therefore, H 0 is approximately diagonal.
Using X, we find where we have dropped off-diagonal terms of order O( 3 )|δm 2 31 | or smaller. Define the matrix Y as  .
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 Fig. 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 , 0)R 13 (θ 13 , δ)R 12 (θ 12 , 0)Q 3 = R 23 (θ 23 , 0)R 13 (θ 13 , δ)R 12 (θ 12 , 0)Q 3 , 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. (A.40) with the effective mixing angles and effective CP-violating phase given approximately by θ 12 ≈ θ 12 = θ 12 + ϕ , (A.58) • δm 2 31 < 0 Case: Using Eq. (3.54), we obtain 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 Fig. 21(b).
When δm 2 31 < 0 we have θ 13 ≈ π 2 in the range β 1 so that Note 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. Thus, we find that the effective mixing matrix U for anti-neutrinos in the case δm 2 31 < 0 can be expressed as Eq. (A.40) with the effective mixing angles and effective CPviolating phase given approximately by θ 12 ≈ θ 12 = θ 12 + ϕ , θ 13 ≈ θ 13 = θ 13 + φ ,
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. Figs. 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 Fig. 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 = 40 • , δ = 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 Fig. 28. All the plots in Fig. 28 have been generated assuming δ = 0 • and inverted mass hierarchy (δm 2 31 < 0). We take the same choices of η and ζ like in Fig. 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 . Fig. 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.

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 Fig. 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 η = −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. Fig. 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. . ν µ → ν 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.