Evolution of Neutrino Mass-Mixing Parameters in Matter with Non-Standard Interactions

We explore the role of matter effect in the evolution of neutrino oscillation parameters in the presence of lepton-flavor-conserving and lepton-flavor-violating neutral-current non-standard interactions (NSI) of the neutrino. We derive simple approximate analytical expressions showing the evolution/running of mass-mixing parameters in matter with energy in the presence of standard interactions (SI) and SI+NSI (considering both positive and negative values of real NSI parameters). We observe that only the NSI parameters in the (2,3) block, namely $\varepsilon_{\mu\tau}$ and $(\gamma - \beta) \equiv (\varepsilon_{\tau\tau} - \varepsilon_{\mu\mu})$ affect the running of $\theta_{23}$. Though all the NSI parameters influence the evolution of $\theta_{13}$, $\varepsilon_{e\mu}$ and $\varepsilon_{e\tau}$ show a stronger impact at the energies relevant for DUNE. $\theta_{12}$ quickly approaches to $\sim$ $90^{\circ}$ with increasing energy in both SI and SI+NSI cases. The change in $\Delta m^2_{21,m}$ is quite significant as compared to $\Delta m^2_{31,m}$ both in SI and SI+NSI frameworks. Flipping the signs of the NSI parameters alters the way in which mass-mixing parameters run with energy. We demonstrate the utility of our approach in addressing several important features related to neutrino oscillation such as: a) unraveling interesting degeneracies between $\theta_{23}$ and NSI parameters, b) estimating the resonance energy in presence of NSI when $\theta_{13}$ in matter becomes maximal, c) figuring out the required baselines and energies to have maximal matter effect in $\nu_{\mu}$ $\rightarrow$ $\nu_{e}$ transition in the presence of different NSI parameters, and d) studying the impact of NSI parameters $\varepsilon_{\mu\tau}$ and $(\gamma - \beta)$ on the $\nu_{\mu} \to \nu_{\mu}$ survival probability.

To obtain a better understanding of the neutrino oscillation probabilities as functions of baseline L and/or neutrino energy E in the presence of SI or SI+NSI, it is quite important in the first place to have a clear knowledge on how various mixing angles and mass-squared differences get modified in matter with energy for a given baseline. Simple approximate analytical expressions showing the evolution/running of mass-mixing parameters in matter with energy in the presence of SI and SI+NSI allow us to address several important features that show up in neutrino oscillation in a more general and transparent fashion. This simple and more intuitive way to understand the neutrino oscillation phenomena will likely pave a way to disentangle the various non-trivial correlations/degeneracies that may be present among the various oscillation and NSI parameters. This paper addresses several pressing issues along this direction.
There exist several studies in the literature investigating how the presence of SI and NSI affect the evolution of effective neutrino oscillation parameters (the mixing angles, mass-squared differences, and CP-violating phase) in matter with energy, and eventually how they modify the oscillation probabilities [55,57,[59][60][61][68][69][70][71][72][73][74][75]. In Refs. [71,72,74,75], the authors diagonalize analytically the three-flavor propagation Hamiltonian in constantdensity matter to obtain the exact expressions for the modified mass-mixing parameters in the presence of SI. The authors in Ref. [73] make use of the Cayley-Hamilton approach with a plane wave approximation to derive the expressions for the modified mass-mixing parameters without performing the actual diagonalization of the Hamiltonian. They also briefly discuss how these oscillation parameters get modified with the strength of SI. In Ref. [55], the author diagonalizes the neutrino propagation Hamiltonian in the presence of SI by applying successive rotations and obtain the expressions for the modified mass-1 They appear into the picture due to the Standard Model (SM) W -exchange interactions between the ambient matter electrons and the propagating electron neutrinos, which is popularly known as the 'MSW effect' [30,47,48]. 2 In Ref. [49], the authors performed a detailed comparative study between different expansions for neutrino oscillation probabilities in the presence of SI in matter. They also studied the accuracy and computational efficiency of several exact and approximate expressions for neutrino oscillation probabilities in the context of long-baseline (LBL) experiments. mixing parameters. In Ref. [75], the authors make use of the relations between the Jarlskog invariants in vacuum and matter (Naumov-Harrison-Scott identities [76][77][78]) to derive the expressions for modified mass-mixing parameters in the presence of SI in constant-density matter. In Ref. [69], the authors adopt a perturbative approach towards the SI and NSI effects and discuss the possible modifications of the mass-mixing parameters. In most of these studies, the authors extract the expressions for modified mass-mixing parameters in order to obtain approximate analytical expressions for the neutrino oscillation probabilities.
Using the Jacobi method [79], the authors in Ref. [59] show that the matter effect on neutrino oscillation due to SI could be assimilated into the evolution of the effective mixing angles θ 12 and θ 13 , and the effective mass-squared differences in matter as functions of the Wolfenstein matter term 2 √ 2G F N e E, while the effective values of θ 23 and δ CP remain unaltered. Here, G F is the Fermi muon decay constant, N e is the ambient electron number density, and E is the energy of the neutrino. They obtain the approximate neutrino oscillation probabilities by simply replacing the mass-mixing parameters in the expressions for the probabilities in vacuum with their running in-matter counterparts. Similar approach is adopted by the authors in Ref. [70] to show the evolution of mass-mixing parameters in the presence of lepton-flavor-conserving, non-universal NSI of the neutrino.
In the present work, we perform successive rotations to almost diagonalize the propagation Hamiltonian in the presence of SI and SI+NSI and derive simple approximate analytical expression for the effective mass-mixing parameters in constant-density matter. While deriving our expressions, we retain the terms of all orders in sin θ 13 and α (the ratio of solar and atmospheric mass-squared differences, ∆m 2 21 /∆m 2 31 ) which are quite important in light of the large value of θ 13 . In our study, we also entertain all possible allowed values of θ 23 in vacuum. As far as NSI are concerned, we consider all possible lepton-flavorconserving and lepton-flavor-violating neutral-current (NC) NSI at-a-time in our analysis which affect the propagation of neutrino in matter. We discuss many salient features of the evolution of oscillation parameters with energy for some benchmark choices of baseline and study in detail how these mass-mixing parameters get affected by various combinations of NSI parameters. Our simple analytical expressions enables us to explore the possible degeneracies between θ 23 (which still has large uncertainty) and NSI parameters for a given choice of neutrino mass ordering in a simple manner. For the first time, we show how the famous MSW-resonance condition (θ 13 in matter becomes 45 • ) [30,47,48,80] gets altered in the presence of NC-NSI. We demonstrate how the simple approximate analytical expressions for the running of oscillation parameters in matter help us to estimate the baselines and energies for which we have the maximal matter effect in ν µ → ν e oscillation channel in the presence of various NSI parameters. For simplicity, we perform our calculations in a CP-conserving scenario where the standard Dirac CP phase δ CP and the phases associated with the lepton-flavor-violating NSI parameters are assumed to be zero. We consider both positive and negative values of real NSI parameters in our analysis.
We plan this paper in the following fashion. We start Sec. 2 with a brief discussion on the theoretical formalism of NSI. This is followed by a short summary of the existing bounds on the NC-NSI. In Sec. 3, we describe our method of approximately diagonalizing the effective neutrino propagation Hamiltonian in the presence of all possible NC-NSI in constant-density matter. Subsequently, we derive the expressions for the modified massmixing parameters. In Sec. 4, we study the evolution of θ 23 , θ 13 , and θ 12 in matter with energy in detail for some benchmark choices of baseline and analyze the role of various NSI parameters on their running. We illustrate the impact of SI and various NSI parameters on the running of two modified mass-squared differences in Sec. 5. In Sec. 6, using the expressions for modified mass-mixing parameters, we estimate for the first time a simple and compact expression for the θ 13 -resonance energy in the presence of all possible NC-NSI parameters and identify the NSI parameters that significantly affect the θ 13 -resonance energy. We devote Sec. 7 to exhibit the utility of our approach in determining the baselines and energies for which we can achieve the maximal matter effect in ν µ → ν e transition in the presence of various NSI parameters. Section 8 describes how the NSI parameters in the (2,3) block affect ν µ → ν µ disappearance channel. Finally, we summarize and draw our conclusions in Sec. 9.

Theoretical Formalism of NSI
NSI which arise naturally in most of the neutrino mass models can be of charged-current (CC) or neutral-current (NC) in nature. Both of them can be described with a dimensionsix operator in the four-fermion effective Lagrangian [30,34,40], where, P C indicates the chiral projection operators P L or P R . The dimensionless coefficients ε f C αβ in Eq. 2.1 denote the strength of NC-NSI between the leptons of flavors α and β (α, β = e, µ, τ ), and the first generation fermions f ∈ {e, u, d}. In Eq. 2.2, the dimensionless coefficients ε f f C αβ indicate the strength of CC-NSI between the leptons of α and β flavors (α, β = e, µ, τ ), and the first generation fermions f = f ∈ {u, d}. The hermiticity of these interactions imposes the following conditions: The CC-NSI modify the production and detection of neutrinos and may also lead to charged-lepton flavor violation. The NC-NSI, on the other hand, affect the propagation of neutrinos. Since the coupling strength ε f C αβ enters into the Lagrangian only through vector coupling, we can write ε f αβ = ε f L αβ + ε f R αβ . It is worthwhile to mention here that models employing scalar mediators [81] or other spin structures [82] are also available in the literature. Beyond a simplified model approach, many UV complete models for NSI have also been explored (see, for instance, [83][84][85][86][87]). For a recent comprehensive review of the NSI, see [43]. Using Eqs. 2.1 and 2.2 and the well-known relation G F / √ 2 g 2 W /8m 2 W , it can be shown that the effective NSI parameters (ε) are proportional to m 2 W /m 2 X [62,65,88], where g W is the coupling constant of the weak interaction, m W is the W boson mass ( 80 GeV ∼ 0.1 TeV), and m X is the mass scale where NSI are generated. Thus it can easily be observed that for m X ∼ 1 TeV, the NSI parameters are of the order of 10 −2 .
In the present work, we concentrate on the NC-NSI which appear during neutrino propagation through matter. Here, the effective NSI parameter can be written in the following fashion Here, N f is the first generation (e, u, d) fermion number density in the ambient medium. The effective Hamiltonian for neutrinos propagating in matter in presence of all the lepton-flavor-conserving and lepton-flavor-violating NC-NSI can be written as where, ∆m 2 21 (≡ m 2 2 − m 2 1 ) and ∆m 2 31 (≡ m 2 3 − m 2 1 ) are the solar and atmospheric masssquared differences, respectively. U is the 3 × 3 unitary Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix in vacuum [89][90][91], which can be parametrized using the three mixing angles: θ 12 , θ 23 , θ 13 , and one Dirac-type CP phase δ CP (ignoring Majorana phases) in the following fashion In Eqn 2.5, V CC is the standard W -exchange interaction potential in matter which can be expressed as where Y e = N e /(N p + N n ) is the relative electron number density of the medium and ρ avg is the line-averaged constant matter density. For the Earth matter which is the focus of our paper, it is safe to assume neutral and isoscalar matter, i.e. N n ≈ N p = N e . Under these assumptions, the relative electron number density inside the Earth turns out to be Y e ≈ 0.5. The (1,1) element of the effective Hamiltonian H f (see Eq. 2.5) contains the term ε ee V CC which gets simply added to the standard matter effect term. Since it can mimic the role of standard interaction, it is a wise choice to subtract a common physical phase I(≡ ε ee V CC ) from the right-hand side (R.H.S.) of Eq. 2.5. Then, the effective Hamiltonian takes the form where, ∆ 31 ≡ ∆m 2 31 /2E, α ≡ ∆m 2 21 /∆m 2 31 ,Â ≡ 2EV CC /∆m 2 31 . We define the effective lepton-flavor-conserving diagonal NC-NSI parameters as β ≡ ε µµ − ε ee and γ ≡ ε τ τ − ε ee . We now briefly discuss the present constraints on the effective NC-NSI parameters obtained from the global fit of neutrino oscillation data [92]. Using Eq. 2.4, we can write, where, Y n is the average neutron/proton ration inside the Earth. According to Ref. [92], Y n = 1.051. Here, we have taken into account the fact that N u = 2N p + N n and N d = N p + 2N n , which in turn imply that ε p αβ = 2ε u αβ + ε d αβ and ε n αβ = ε u αβ + 2ε d αβ . Note that the contribution from ε e αβ is not considered in the global 3ν analysis in the presence of NC-NSI parameters [92]. Now, we use the bounds (2σ) on ε u αβ and ε d αβ from the global fit analysis [92] and list the subsequent 2σ bounds on the effective NSI parameters ε αβ in Table 1.

Diagonalization of the Effective Hamiltonian in the presence of NSI
Here, we derive the approximate analytical expressions for the fundamental oscillation parameters in matter considering all possible lepton-flavor-conserving and lepton-flavorviolating NC-NSI 3 which are real i.e., all the phases associated with the non-diagonal elements of the NSI matrix are assumed to be zero.
In order to simplify the subsequent calculations, we perform our analysis in the CPconserving scenario i.e., we take the standard Dirac CP phase δ CP to be zero. The elements 3 The authors in Ref. [93] derived similar expressions in the context of a particular beyond the Standard Model (BSM) scenario where they considered the presence of long-range flavor-diagonal NSI appearing due to abelian Le-Lµ symmetry. In the present work, we adopt a model independent approach and introduce all possible NSI parameters at-a-time in the framework. It allows us to study the evolution of massmixing parameters in a more generalized scheme considering all possible NSI parameters which has rich phenomenological implications in neutrino oscillation. In the above expressions, we use the abbreviations: cos θ ij → c ij , sin θ ij → s ij , and retain the terms of all orders in sin θ 13 and α which are quite essential in light of the large value of θ 13 . To find the effective mixing angles and mass-squared differences in the presence of Earth matter potential (V CC ) and all possible NC-NSI parameters, we need to diagonalize the effective Hamiltonian H f in Eq. 2.8. We approximately diagonalize H f by applying three successive rotations R 23 (θ m 23 ), R 13 (θ m 13 ), and R 12 (θ m 12 ), where R ij (θ m ij ) is the rotation matrix for the (i, j) block with the rotation angle θ m ij . The product of these rotation matrices construct a 3 × 3 unitary matrix where, the off-diagonal terms after the final rotation are quite small (∼ 10 −8 ) and can be safely neglected. Below, we give the expressions for the mixing angles in matter that we derive by equating the small off-diagonal elements to zero after each rotation during the diagonalization process:   where, ∆θ 23 ≡ θ 23 − θ m 23 is the deviation of the modified mixing angle θ 23 from its vacuum value. In the above equations, λ 1 , λ 2 , and λ 3 take the following forms: (3.12) (3.13) (3.14) Note that throughout the entire paper, we consider the propagation of neutrinos inside the Earth and assume normal mass ordering 4 (NMO). In case of antineutrino propagation, one has to reverse the sign of V CC in the above equations which in turn reverses the sign ofÂ. Similarly, to get the corresponding expressions for the inverted mass ordering (IMO), one has to flip the sign of α as well as the sign ofÂ in Eqs. 3.9 to 3.14.

Evolution of Mixing Angles in the presence of NSI
In the present section, we study in detail how the effective mixing angles in matter θ m 12 , θ m 13 , and θ m 23 (we derive their expressions in Sec. 3) get modified as functions of energy and baseline in the presence of all possible NC-NSI. For this study, we consider the three-flavor vacuum oscillation parameters as given in Table 2. To show our results, we consider two benchmark values of the NSI parameters: 0.2 and -0.2.

Running of θ m 23
Approximate analytical expression describing the evolution of the effective mixing angle θ m 23 is given in Eq. 3.9. We can further simplify this expression by neglecting the small 4 There are two possible patterns of neutrino masses: a) m3 > m2 > m1, called normal mass ordering (NMO) where ∆m 2 31 > 0 and b) m2 > m1 > m3, called inverted mass ordering (IMO) where ∆m 2 31 < 0.  Table 2. terms which are proportional to αs 13 ∼ 10 −3 in Eq. 3.9, which enable us to extract the useful physics insights related to the running of θ m 23 in a more concise fashion. With this approximation, the expression showing the evolution of θ 23 in matter in the presence of NSI takes the form where, γ − β = ε τ τ − ε µµ . Two important features emerge from this simplified expression.
• In the limiting case of all NSI parameters equal to zero (which in this case removes the standard matter effectÂ also), one would get back the vacuum mixing angle (i.e., θ m 23 = θ 23 ) irrespective of energy, baseline, and the octant of θ 23 . In other words, it implies that θ m 23 does not run in the presence of standard matter effect.
Note that, in the exact expression of θ m 23 in Eq. 3.9, due to the presence of the tiny terms proportional to αs 13 , θ m 23 slightly deviates from its vacuum value even in the presence of SI.
In Fig. 1, we show the running of θ m 23 (using Eq. 3.9) with energy in presence of NSI parameters (γ − β), ε µτ , taken one-at-a-time for a baseline corresponding to the DUNE experiment i.e. 1300 km. The left column shows the effect of NSI parameter (γ − β) while the right column corresponds to the effect of ε µτ . The black curves in each column depict the SI case for three possible values of θ 23 in vacuum, namely higher octant (θ 23 = 50 • ), maximal mixing (θ 23 = 45 • ), and lower octant (θ 23 = 40 • ). As discussed above, value of θ m 23 in SI case remains almost equal to the value of θ 23 in vacuum. Only very small deviations from the vacuum value of θ 23 can be observed due to the presence of terms proportional to αs 13 in Eq. 3.9, which are neglected in Eq. 4.1. The solid (dashed) red curves in the left column of Fig. 1 illustrate the presence of (γ − β) with a benchmark value of 0.2 (-0.2). We observe that for all the three values of θ 23 mentioned above, θ m 23 monotonically decreases (increases) with energy when (γ − β) is present with a positive (negative) value. In the right column, the solid (dashed) blue curves depict the case when only ε µτ is present with a benchmark value of 0.   non-zero. The four colored curves in each panel illustrate the effect of the four possible sign combinations of (γ − β) and ε µτ while the black curve shows the SI (with standard matter effect and no NSI) case, as shown in the legend. As before, three scenarios of the vacuum mixing angle θ 23 are considered: higher octant (θ 23 = 50 • ), maximal mixing (θ 23 = 45 • ), and lower octant (θ 23 = 40 • ). We note from Fig. 2 that in the presence of (γ − β) with a negative (positive) sign, θ m 23 monotonically increases (decreases) with energy irrespective of the sign of ε µτ and the octant of θ 23 . We also observe that for lower (higher) octant, the decrease (increase) is the steepest when (γ − β) is positive (negative) with negative value of ε µτ . For maximal mixing, the running of θ m 23 appears symmetric around the SI case since the term with cos 2θ 23 in the denominator of Eq. 3.9 vanishes.    Table 2 and we assume NMO.
To show a correlation between NSI strength and the value of θ 23 in vacuum, we have shown in Fig. 3, the evolution of θ m 23 in the plane of [θ 23 −ε µτ ] (top panels) and [θ 23 −(γ −β)] (bottom panels). We demonstrate the effect of baseline by choosing three different baseline lengths as 1300 km, 5000 km, and 8000 km in the three columns, respectively. For each baseline, the energy is chosen as the corresponding value near the first oscillation maximum (E max ) for ν µ → ν e oscillation. For the baseline of 1300 km, we see that θ m 23 decreases (increases) from the vacuum value (θ 23 ) for a positive (negative) ε µτ at higher octant.
However, an opposite trend can be observed at lower octant. For maximal mixing, θ m 23 does not change in presence of ε µτ only. These features are more pronounced for higher baselines since the NSI effect (proportional to matter density) gets enhanced. In the bottom row, in the presence of positive (negative) value of (γ − β), θ m 23 decreases (increases) from the vacuum value, irrespective of the octant or maximal mixing. Larger baselines manifest it more clearly as evident from the steeper slant of the boundaries between different colors.
As mentioned earlier, we have assumed normal mass ordering (NMO) for our analysis. In case of inverted mass ordering (IMO) with neutrino (ν, IMO), the effect of each NSI parameters in θ m 23 running is reversed (i.e., if θ m 23 increases with energy in presence of a particular NSI parameter with normal ordering of mass, in case of inverted mass ordering θ m 23 will decrease with energy). This happens since the termÂ associated with each NSI parameter changes its sign in case of IMO. Also, in case of antineutrino propagation with inverted mass ordering (ν, IMO), running of θ m 23 is almost the same as that of neutrino propagation with NMO (ν, NMO). This is because of the fact that in both cases, sign of A is the same.

Running of θ m 13
Eq. 3.10 shows the running of θ m 13 in matter with NSI. We note that all five NSI parameters as well as the standard matter effect (Â) have impact on the running 5 of θ m 13 . It is observed that the value of θ 23 in vacuum (when it is between 40 • and 50 • ) has a very small effect on the running of θ m 13 . So, we simplify the expression for our understanding by assuming that the mixing angle θ 23 in vacuum is maximal i.e.,  such that at lower energy (E 6 GeV), it is almost constant when only one of them is present. It can be explained by the fact that both numerator and denominator of R.H.S. in Eq. 4.2 decreases with energy when ε eµ and/or ε eτ are negative, such that the overall value of θ m 13 remains almost constant at that energy range. However, at higher energy (∼ 10 GeV) value of the denominator is so small that the overall effect led to the rapid increase in the magnitude of θ m 13 with energy. As we can see from Eq. 4.2, in presence of both ε eµ and ε eτ with a negative sign, the numerator decreases faster with energy compared to the previous case due to the additive effect of two NSI parameters. Consequently, value of θ m 13 decreases with energy from its vacuum value, and becomes negative (at E 6 GeV) when the numerator becomes negative.
In the case of IMO, the behavior of θ m 13 in SI as well as in SI+NSI case is significantly different from the NMO case for neutrino. It can be understood from the (λ 3 −Â) term in the denominator of Eq. 4.2. SinceÂ changes its sign, the denominator increases with energy, consequently the value of the θ m 13 decreases from its vacuum value. However, in case of antineutrino (ν) propagation and inverted mass ordering (ν, IMO), sinceÂ does not change its sign, running of θ m 13 is almost similar to neutrino (ν, NMO) case. In the presence of NSI parameters in (2,3) block, λ 1 , λ 2 and c m 13 undergo mild change,retaining almost the same features as that of SI. The presence of ε eµ (ε eτ ) however, adds up to the numerator of Eq. 4.4 7 and the value of θ m 12 at which it saturates, shifts down (up). When both ε eµ and ε eτ are present, they cancel their effect due to the relative sign between them and the running of θ m 12 almost coincides with SI scenario. In the bottom row, we have shown the running of θ m 12 for the NSI with negative strength. Since running of θ m 12 very mildly depend on NSI parameters from the (2,3) sector (bottom left panel), the sign of these NSI parameters do not have any significant effect. In the bottom right panel, we see that role of ε eµ and ε eτ is reversed when the sign of the NSI parameter is changed. Interestingly, at energies around 10 GeV, sudden decrease (increase) of θ m 12 can be observed in the presence of NSI parameter ε eµ (ε eτ ) with negative strength. It happens due to the presence of the term cos θ reduces rapidly to a very small value around that energy (see Fig. 4 and related discussion in Subsec. 4.2). Unlike θ m 13 , the θ m 12 running shows similar behavior in SI as well as in SI+NSI cases for neutrino propagation with IMO (ν, IMO). Also, it shows completely different behavior in case of antineutrino propagation with inverted mass ordering (ν, IMO). It can be understood from the fact that in case of IMO, sign of first and third terms in the numerator of Eq. 4.4 gets flipped, and in the denominator, the sign of λ 1 gets changed. Since the effect from other remaining terms are very small, both numerator and denominator change their sign, and as a result, θ m 12 remains the same as in the case of (ν, NMO). In case of (ν, IMO), only first term in the numerator changes its sign, λ 1 in the denominator remains the same as in case of (ν, NMO). As a result, we see a completely different behavior of θ m 12 .

Evolution of Mass-Squared Differences in the presence of NSI
where, we assume θ 23 = 45 • and we are already familiar with the expressions of θ m ij and λ i . Using the above equations, we can obtain the approximate analytical expressions for the modified mass-squared differences ∆m  Table 2. We assume θ23 = 45 • and NMO.
in the presence of positive (negative) NSI with strength 0.2. SI case shows steady increase with energy, -reaching a value of an order as high as 10 −3 eV 2 from its vacuum order of magnitude 10 −5 eV 2 . In other words, the running of ∆m 2 21,m can make itself comparable in magnitude with that of ∆m 2 31,m . In the presence of positive (top row) or negative (bottom row) NSI (except for negative ε µτ or taking negative β, γ, ε µτ together), the behavior of ∆m 2 21,m does not show significant deviation in magnitude from SI case. But interestingly, depending on the sign of NSI parameter, the magnitude of ∆m 2 21,m in presence of NSI is slightly higher or lower than in presence of SI. In presence of negative ε µτ (when present singly or together with negative β and γ), we see a deviation from SI case at higher energy which can be understood from the variation of λ 1 with energy. With negative ε eµ or ε eτ , however, at E 10.5 GeV ∆m 2 21,m becomes almost constant. It happens due to a sudden increase in the value of θ m 13 around that energy which results in saturation of the value of λ 1 .
In the case of IMO, running of ∆m 2 21,m is almost the same as (ν, NMO) case for both neutrino and antineutrino propagation. However, IMO leads to a significant change in the running of ∆m 2 31,m for both neutrino and antineutrino propagation which is obvious because the vacuum value ∆m 2 31 changes its sign. becomes maximal (45 • ). We note that this resonance is independent of the value of ε eµ or ε eτ (as evident from the right panels of Fig. 4)  In the above equation, we neglect the small terms proportional to αs 2 13 , (γ − β) 2Â2 and the cross-term proportional to αÂ(γ − β)s 13 . Finally, we get the following simpler expression for λ 3 , It is noteworthy to mention that for SI case, we get λ 3 c 2 13 . Putting this back in Eq. 6.1 and using OMSD approximation, we easily obtain the well-known expression for resonance in Eq. 6.2. Equating Eqs. 6.1 and 6.4, we obtain the following final expression for the resonance energy, .
(6.5) The term in the square bracket in the R.H.S. of Eq. 6.5 is the correction over Eq. 6.2. The term 1 2 (β+γ+2ε µτ ) is the correction induced by the presence of NSI, while αs 2 12 c 2 13 / cos 2θ 13 is the modification induced by relaxing the OMSD approximation. Thus it is now also clear analytically that θ m 13 -resonance gets affected only by the NSI parameters in the (2,3) block and not by ε eµ or ε eτ .   Fig. 8 shows the line-averaged constant Earth matter density (ρ avg ) for a given baseline L obtained from the PREM profile [94]. In both the panels of Fig. 8, we indicate by three gray shades, the baselines when it touches the three interior layers of the Earth: crust, mantle, and core. Since ρ avg shows an increase in magnitude (thus increasing V CC in Eq. 6.5) with L, the values of E NSI res itself decreases with L, following similar pattern (for both SI and NSI). As it is also clear from Eq. 6.5, a positive (negative) value of the NSI parameters β, γ, or ε µτ shifts the magnitude of E NSI res to a higher (lower) value than the SI case. Eq. 6.5 also tells us, if it turns out that the NSI parameters are present in Nature with such magnitudes that β + γ = −2ε µτ , then the correction due to NSI vanishes. In that case, if we ignore the minor correction induced by αs 2 12 c 2 13 / cos 2θ 13 term we obtain, E NSI res E SI res .
One of the most important oscillation channel that is probed in LBL experiment is ν µ → ν e appearance channel. This channel will play a significant role in determining the value of the CP phase, neutrino mass ordering, octant of θ 23 from various upcoming neutrino oscillation experiments. So, in this section, we are interested in studying the effect of NSI on ν µ → ν e transition probability maxima at various baselines (L) through the Earth-matter with neutrino beam having energy (E) in the GeV range. In order to study this, in Fig. 9, we plot the ν µ → ν e transition probability in (E -L) plane in SI  Table 2 with θ23 = 45 • and NMO.
case and SI+NSI cases considering a benchmark value of 0.2 for the strength of the NSI parameters. We use the vacuum probability expression for ν µ → ν e appearance [59] and replace the vacuum oscillation parameters with their modified counterparts in matter with SI and NSI (Eqs. 3.9-3.11 and Eqs. 5.1-5.3) and use this modified probability expression for plotting Fig. 9. θ 23 in vacuum is considered to be maximal. We check that Fig. 9 shows very good agreement in both SI and SI+NSI cases with the exact three-flavor oscillation probabilities which are calculated numerically using the GLoBES software [95,96]. Top left panel shows the SI case where no NSI are taken into account. Here, it is observed that the region of maximum appearance probability occurs for the baseline almost passing through the core and the mantle boundary. However, in presence of ε eµ (bottom left panel) or ε eτ (top right panel) this region shifts towards lower baselines. In case of ε µτ (bottom right panel), this region remains almost the same as in the SI case. To show the effect of NSI with negative strength, we similarly plot the oscillograms in Fig. 10 in SI case and in the presence of negative NSI with strength 0.2. Huge differences in the oscillation patterns can be observed in case of ε eµ and ε eτ . Unlike Fig. 9, there is no such region of the maximum transition probability in Fig. 10. In Fig. 9, we notice that for positive values of ε eµ and ε eτ , a significant enhancement in the ν µ → ν e transition probability as compared to SI case, for some choices of L and E where we have large matter effects. On the other hand, in Fig. 10, for negative choices of ε eµ and ε eτ , we see a large depletion in ν µ → ν e transition probability for some choices of L and E where matter effect is suppressed. Now, we make an attempt to understand these features with the help of approximate analytical expressions. After replacing the vacuum oscillation parameters with their modified counterparts in the ν µ → ν e transition probability as mentioned above, we simplify it further by using the approximation that θ m 12 almost saturates to π/2 (see Fig. 5 and the related discussion in Subsec. 4.3). As a result, we obtain the following simplified expression that helps us to explain the broad features observed in Figs (7.1) In Fig. 11, we plot Eq. 7.1 with energy and also the contribution from each term T 1 , T 2 , and T 3 , separately. It is clear from the figure that the variation of T 1 with energy is very less compared to the variation of T 2 and T 3 . Thus, the energy at which maximum of T 1 T 2 T 3 occurs is the same as that of T 2 T 3 . In other words, the maximum of P m νµ→νe occurs at an energy determined by T 2 and T 3 , not T 1 . This feature is also valid for any other baselines. So, P m νµ→νe is maximum when both the following two conditions are satisfied simultaneously.

(7.4)
Note that under the OMSD approximation, the resonance energy condition in Eq. 6.1 takes a very simple form: (λ 3 −Â − s 2 13 ) = 0 and we make use of this expression in Eq. 7.3 to obtain Eq. 7.4, which exactly matches with the expression derived by the authors in Ref. [98]. Replacing ∆m 2 32,m in Eq. 7.2, we finally have the following condition for the maximal ν µ → ν e transition probability in the presence of NC-NSI.
ρ avg L NSI (2n + 1)π × 5.18 × 10 3 tan 2θ 13   The second factor in the R.H.S. of Eq. 7.8 is the correction introduced by the NSI parameters. As shown in Fig. 11, the modified θ m 23 does not run significantly and is also close to 45 • for maximal mixing of θ 23 . Using the approximation θ m 23 45 • , we simplify Eq. 7.8 to the following. ρ avg L max NSI ρ avg L SI 1 1 − (β + γ + 2ε µτ ) /2 + √ 2(ε eµ + ε eτ )/ tan 2θ 13 km g/cm 3 . (7.10) In Fig. 12, we plot the R.H.S. of Eq. 7.10 in presence of SI and also in presence of NSI parameters (β, ε µτ , ε eτ taken one-at-a-time with a magnitude of 0.2) in the ρ avg and L plane. In the presence of NSI parameter γ (ε eµ ), Eq. 7.10 is same as in the presence of β (ε eτ ). Again, three gray regions correspond to the baseline length passing through the crust, crust-mantle, and crust-mantle-core. In the same figure, we also show the line-averaged constant Earth matter density according to the PREM profile (yellow curve). The points of intersections of the hyperbolic curves (corresponding to Eq. 7.4 for the SI case and Eq. 7.10 for the SI+NSI cases) with the yellow curve give the baseline lengths, required to achieve the maximum of P m νµ→νe . When there are no NSI present in the scenario, the baseline length required for maximum transition probability is around 10600 km, which almost  touches the core of the Earth. This same feature is also observed in the oscillogram plot in Fig. 9 which is plotted using full three-flavor ν µ → ν e oscillation probability expression. When NSI parameters from (2,3) block with positive (negative) strength are present oneat-a-time, the required baseline length increases (decreases) slightly. But interestingly, when the NSI parameter ε eµ or ε eτ is present, the required baseline length for maximum ν µ → ν e transition decreases drastically to around 6700 km (intersection between red and yellow curve in Fig. 12), which passes through only crust and mantle of the Earth 8 . It is evident from Eq. 7.10 that since the role of ε eµ and ε eτ are on the same footing, the presence of ε eµ induces an effect identical to that of ε eτ with the same magnitude. But the oscillograms (Figs. 9 and 10) for ε eµ and ε eτ look quite different. This is because of the fact that θ m 12 saturates to a value higher or lower than 90 • in the presence of ε eµ or ε eτ (see Fig. 5). Since θ m 12 is not exactly 90 • , we have non-zero contributions from some other terms in ν µ → ν e oscillation probability expression, which affect the oscillograms in the presence of ε eµ and ε eτ in a different fashion. In the presence of negative ε eµ and ε eτ , we observe from Fig. 10 that we no longer achieve the maximum transition in ν µ → ν e oscillation channel. It is because of the fact that in this case, the baseline length required for the maximum ν µ → ν e appearance probability turns out to be longer than the Earth's diameter (see Eq. 7.10). Therefore, it is not possible to attain the maximum ν µ → ν e transition inside the Earth for negative values of ε eµ and ε eτ as evident from Fig. 10. We observe from Fig. 9 and Fig. 10 that in the presence of non-zero ε µτ , there are slight changes in L and E as compared to SI case for which we obtain maximum possible ν µ → ν e transition.
8 Impact of NSI in ν µ → ν µ disappearance channel So far, we have focused on ν µ → ν e appearance channel which is one of the most important channel probed in LBL experiments. However, another crucial channel, ν µ → ν µ disappearance channel can be probed in LBL and atmospheric neutrino experiments. This channel can play an important role in precision measurement of the atmospheric oscillation parameters. In this section, we discuss the effect of NSI in ν µ → ν µ survival probability. Since NSI parameters from the (2,3) block have significant impact on this channel [68,100], only these NSI parameters have been considered. To get the broad feature, we simplify the analysis by assuming ∆ 21 0 and θ 13 0. Under these approximations, ν µ → ν µ disappearance probability expression reduces to [101,102]  (1 + ε µτÂ ) , (8.4) where, we use the approximation θ m 12 → π/2 in the expression of m 2 1,m in Eq. 5.3. So, using Eqs. 8.3 and 8.4, ν µ → ν µ disappearance probability in presence of NSI parameters from (2,3) sector can be written as, (1 + 2ε µτÂ ) .
If we only consider the off-diagonal NSI parameter ε µτ , the expression boils down to the simplied expression already derived in [101]. From the approximate expression in Eq. 8.5, some broad features about the impact of NSI on the ν µ → ν µ survival channel can be observed. We see that the parameter (γ − β) always appears in second order in Eq. 8.5, while other NSI parameter ε µτ has a linear dependence. For the same reason, the sign of (γ − β), unlike the sign of ε µτ , does not affect the disappearance probability. Since the strength of NSI parameters are not very large, it is expected that the impact of (γ − β) will be always small compared to ε µτ . To find out whether these features remain intact even if we assume non-zero θ 13 and finite ∆m 2 21 , in Fig. 13, we have plotted the ν µ → ν µ survival probability as a function of baseline (x-axis) and energy (y-axis) commonly known as oscillogram plot. We first consider the full three-flavor vacuum expression of ν µ → ν µ survival probability without any approximation [59] and replace the vacuum parameters with their modified expressions in matter with NSI which we have derived in this work. As before, θ 23 in vacuum is assumed to be 45 • . In the top row, we compare the survival probabilities in presence of the parameter (γ − β) with negative (left column) and positive (right column) values and compare with the SI case (middle column). We see that there is not any notable variation in oscillation pattern except for the magnitude of the survival probability, which decreases at high energies (E) and baselines (L). Some small differences in the pattern appear at some (L, E) region. It appears due to non-zero θ 13 which brings the matter effect into the picture and the finite value of ∆ 21 which gives correction due to solar term. As predicted from Eq. 8.5, the sign of (γ − β) does not affect the ν µ → ν µ disappearance probability.
In the bottom row of Fig. 13, we plot the same but in presence of ε µτ with negative (positive) value in the extreme left (right) panel and show the results for SI case in the middle column. It is observed that the presence of ε µτ can lead to significant differences in the pattern of ν µ → ν µ disappearance probability as compared to SI case. When ε µτ is positive (negative), one can observe a significant shift in the oscillation dip (blue regions emerging from the origin) from the SI case towards higher (lower) energies in Fig. 13. This feature can be explained from the approximate expression in Eq. 8.5. In that expression, the value of P m νµ→νµ is mainly determined by the first term in R.H.S. since second term is suppressed by the NSI parameters appearing quadratically. At the first term in R.H.S. of  Table 2.
Eq. 8.5, minimum occurs at higher (lower) energy compared to SI case for a given baseline L when ε µτ is present with positive (negative) strength 9 . Depending on the sign of ε µτ , the regions representing the oscillation dip tend to bend upward or downward with increase in baseline length compared to SI case.

Summary and Concluding Remarks
In this work, we derive the expressions for the evolution of the fundamental mass-mixing parameters in the presence of SI and SI+NSI considering all possible lepton-flavor-conserving and lepton-favor-violating NC-NSI. In order to derive these expressions, we use a method of approximate diagonalization of the effective Hamiltonian by performing successive rotations in (2,3), (1,3), and (1,2) blocks. In our study, we present the results for the benchmark value of the DUNE baseline of 1300 km and also discuss the results for few other baselines. 9 At oscillation dip, the argument of the cosine term in Eq. 8.5 should be approximately equal to (2n + 1)π/2 where n = 0,1,2.... This roughly implies that (1 + 2εµτÂ) ∆m 2 31 L 4E π/2. We consider both positive and negative values of real NSI parameters with benchmark values of ± 0.2.
In the presence of SI only, the 2-3 mixing angle in matter (θ m 23 ) receives a tiny correction which is independent of energy and the strength of the matter potential. It is observed that only the NSI parameters in the (2,3) block, namely ε µτ and (γ − β) ≡ (ε τ τ − ε µµ ) influence the evolution of θ m 23 . In the presence of negative (positive) value of (γ − β), θ m 23 increases (decreases) with energy. For the maximal value of θ 23 in vacuum, the change in θ m 23 is negligible in the presence of ε µτ . If θ 23 belongs to the upper octant then θ m 23 increases (decreases) for negative (positive) choices of ε µτ . We notice a completely opposite behavior if θ 23 lies in the lower octant. We also study the modification in θ m 23 as a function of energy when both the NSI parameters ε µτ and (γ − β) are present in the scenario with their all possible sign combinations. We unravel interesting degeneracies in [θ 23 -(γ − β)] and [θ 23 -ε µτ ] planes for three different combination of L and E and discuss how our simple approximate analytical expression showing the evolution of θ m 23 plays an important role to understand these complicated degeneracy patterns.
In contrast to θ m 23 , θ m 13 is more sensitive in matter in the presence of SI and SI+NSI. Therefore, an accurate understanding of the running of θ 13 in matter is crucial to correctly assess the outcome of the oscillation experiments in the presence of NC-NSI. θ m 13 goes through an appreciable change even in SI case depending on the choice of mass ordering and whether we are dealing with neutrinos or antineutrinos. Compared to SI case, the relative change in θ m 13 for (ν, NMO) is somewhat suppressed (enhanced) in the presence of positive (negative) NSI parameters in the (2,3) block, namely γ ≡ (ε τ τ − ε ee ), β ≡ (ε µµ − ε ee ), and ε µτ . For positive ε eµ and/or ε eτ , θ m 13 for (ν, NMO) approaches the resonance (θ m 13 = 45 • ) faster than SI case, but after crossing the resonance energy, SI takes over. For negative ε eµ or ε eτ , running of θ m 13 is suppressed almost up to the resonance energy and then, it increases very steeply compared to SI case.
As far as the solar mixing angle is concerned, θ m 12 approaches to 90 • (sin θ m 12 → 1, cos θ m 12 → 0) very quickly as we increase the neutrino energy in SI case. While the NSI parameters in the (2,3) block (γ, β, ε µτ ) have minimal impact on the running of θ m 12 , ε eµ and ε eτ affect the evolution of θ m 12 substantially. In the presence of positive (negative) ε eµ , the change in θ m 12 with energy qualitatively remains the same with the saturation value turns out to be around 80 • (100 • ). In the presence of positive (negative) ε eτ , the saturation happens around 100 • (80 • ).
Out of the two mass-squared differences, the evolution of the solar ∆m 2 21,m is quite dramatic as compared to that of the atmospheric ∆m 2 31,m in matter. Both in SI and SI+NSI cases, as we increase the energy and go beyond 10 GeV, the value of ∆m 2 21,m increases to almost 20 times as compared to its vacuum value for 1300 km baseline. At the same time, the value of ∆m 2 31,m does not change much compared to its vacuum value. We demonstrate the utility of our approach in addressing some interesting features that we observe in neutrino oscillation in presence of matter. It is well known that θ m 13 can attain the value of 45 • (MSW-resonance condition) for some choices of L and E in the presence of standard matter effect. Now, in this work, for the first time, we show how the θ 13 -resonance energy gets modified in the presence of NC-NSI with the help of simple, compact, approximate analytical expressions. We observe that only the NSI parameters in the (2, 3) block affects the θ 13 -resonance energy.
We study in detail how the NC-NSI parameters affect ν µ → ν e oscillation probability which plays an important role to address the remaining unknown issues, namely CP violation, mass ordering, and the precision measurement of oscillation parameters. In this paper, for the first time, we derive a simple approximate analytical expression for E, L and its corresponding ρ avg to have the maximal matter effect which in turn gives rise to maximum ν µ → ν e transition in the presence of all possible NC-NSI parameters. This analytical expression reveals that in SI case, maximum ν µ → ν e transition occurs around the baselines almost passing through the core of the Earth ( 10600 km) under the assumption that the θ 13 -resonance energy coincides with the energy that corresponds to the first oscillation maximum. However, in the presence of positive ε eµ or ε eτ in matter, it is observed that the baseline length for maximum ν µ → ν e transition probability gets reduced to a much lower value (L ≈ 6700 km) for the benchmark choices of ε eµ or ε eτ equal to 0.2. On the other hand, if we consider negative ε eµ or ε eτ , the required baseline for maximum transition probability increases to a value that exceeds the diameter of the Earth.
We also study in detail how the NSI parameters in the (2,3) block affect ν µ → ν µ disappearance channel which plays an important role in atmospheric neutrino experiments. We observe that the off-diagonal NSI parameter ε µτ has the dominant effect as compared to the diagonal NSI parameter (γ − β). It happens because (γ − β) appears in the secondorder in the approximate ν µ → ν µ oscillation probability expression, whereas ε µτ shows a liner dependence. Also, the sign of ε µτ has a significant impact on ν µ → ν µ disappearance channel. However, this oscillation channel is not sensitive to the sign of the NSI parameter (γ − β). We hope that the analysis performed in this paper will take our understanding of the evolution of the oscillation parameters in the presence of all possible NC-NSI a step forward.