On the Determination of Leptonic CP Violation and Neutrino Mass Ordering in Presence of Non-Standard Interactions: Present Status

We perform a global analysis of neutrino data in the framework of three massive neutrinos with non-standard neutrino interactions which affect their evolution in the matter background. We focus on the effect of NSI in the present observables sensitive to leptonic CP violation and to the mass ordering. We consider complex neutral current neutrino interactions with quarks whose lepton-flavor structure is independent of the quark type. We quantify the status of the"hints"for CP violation, the mass-ordering and non-maximality of $\theta_{23}$ in these scenarios. We also present a parametrization-invariant formalism for leptonic CP violation in presence of a generalized matter potential induced by NSI.


Introduction
Experiments measuring the flavor composition of neutrinos produced in the Sun, in the Earth's atmosphere, in nuclear reactors and in particle accelerators have established that lepton flavor is not conserved in neutrino propagation, but it oscillates with a wavelength which depends on distance and energy. This demonstrates beyond doubt that neutrinos are massive and that the mass states are non-trivial admixtures of flavor states [1,2], see Ref. [3] for an overview.
The most recent global analysis of the neutrino oscillation data in Ref. [4] (see also [5,6]) in the context of 3ν-mixing provides us with the determination of the three leptonic mixing angles and the two mass differences with a 3σ/3 precision of 3% for θ 12 , 3% for θ 13 , 9% for θ 23 , 5% for ∆m 2 21 and 2.5% for |∆m 2 31 |. Questions still open in the analysis are the maximality and octant of θ 23 , the ordering of the mass eigenstates, and the status of the leptonic CP violating phase δ CP . Some hints in this respect are emerging -with a special role played by the long-baseline (LBL) accelerator experiments T2K [7,8], and NOvA [9, 10] -but without a consolidated statistical significance yet. The latest results show a preference of the normal ordering (NO) at the 2-3σ level, and a best fit for the complex phase at δ CP = 215 • close to maximal CP violation. The clarification of these unknowns is the main focus of the running LBL experiments and its precise determination is at the center of the physics program of the upcoming LBL facilities, in particular the Deep Underground Neutrino Experiment (DUNE) [11] and the Tokay to HyperKamiokande (T2HK) experiment [12].
Under the assumption that the Standard Model (SM) is the low energy effective model of a complete high energy theory, neutrino masses emerge naturally as the first observable consequence in the form of the Weinberg operator [13], the only dimension five operator that can be built within the SM particle content. In this framework the next operators with observable consequences at low energies appear at dimension six. They include four-fermion terms leading to Non-Standard Interactions (NSI) [14][15][16] between neutrinos and matter (for a recent review see [17]).
Neutral Current NSI can modify the forward-coherent scattering (i.e., at zero momentum transfer) of neutrinos as they propagate through matter via the so-called Mikheev-Smirnov-Wolfenstein (MSW) mechanism [14,18]. Consequently their effect can be significantly enhanced in oscillation experiments where neutrinos travel large regions of matter. Indeed, the global analysis of data from oscillation experiments in the framework of mass induced oscillations in presence of NSI currently provides some of the strongest constraints on the size of the NSI affecting neutrino propagation [19][20][21].
In the presence of such NSI, however, the task of exploring leptonic CP violation in LBL experiments becomes enriched (to the point of confusion) by the possible coexistence of new sources of CP violation [22]. Furthermore, the determination of the mass ordering is jeopardized by the presence of an intrinsic degeneracy in the relevant oscillation probabilities due to a new symmetry of the Hamiltonian describing the neutrino evolution in the modified matter potential [19,20,23,24]. This has resulted in an intense phenomenological activity to quantify these issues and to devise strategies to clarify them, first in proposed facilities like the Neutrino Factory [25][26][27] and most recently in the context of the upcoming experiments .
Very interestingly, it has been argued that NSIs can already play a role in the significance of the "hints" mentioned above [51,52]. In particular in Ref. [51] it was pointed out the discomforting possibility of confusing CP conserving NSI with a non-zero value of δ CP in the analysis of ν e andν e appearance results at T2K and NOvA. Clearly such confusion could lead to an incorrect claim of the observation of leptonic CP-violation in a theory which is CP conserving.
We recently performed a global analysis of oscillation data in the presence of NSI relevant to neutrino propagation in matter in Ref. [21]. For simplicity the analysis in Ref. [21] only constrained the CP conserving part of the Hamiltonian and for consistency the observables most sensitive to CP violating effects, i.e., ν e andν e appearance at LBL experiments, were not included in the fit. Consequently the issue of the possible confusion between real NSI and leptonic CP-violation could not be addressed. Furthermore, under the simplifying assumptions employed, the results of Ref. [21] could not yield any conclusion on the status of the determination of the ordering of the states in the presence of NSIs either.
In this paper we extend the analysis in Ref. [21] to account for the effect of NSI in the observables sensitive to leptonic CP violation and to the mass ordering. Our goal is to quantify the robustness of the present "hints" for these effects in the presence of NSI which are consistent with the bounds imposed by the CP-conserving observables. We start by briefly summarizing the formalism and notation in Sec. 2. In doing so we take the opportunity to present a parametrization-invariant formalism for leptonic CP violation in the presence of a generalized matter potential induced by NSI. In Sec. 3 we describe the strategy employed in the study. Finally in Sec. 4 we present the results, and in Sec. 5 we summarize our conclusions. We present some detail of the construction of the basis invariants for CP violation in Appendix A.

Formalism
In this work we will consider NSI affecting neutral-current processes relevant to neutrino propagation in matter. The coefficients accompanying the relevant operators are usually parametrized in the form: where G F is the Fermi constant, α, β are flavor indices, and f is a SM charged fermion. In this notation, ε f αβ parametrizes the strength of the vector part of the new interactions (which are the ones entering the matter potential) with respect to the Fermi constant, ε f αβ ∼ O(G X /G F ). In the second line we make explicit that, as in Ref. [21], we assume that the neutrino flavor structure of the interactions is independent of the charged fermion type, so that one can factorize ε f αβ ≡ ε αβ ξ f . Ordinary matter is composed of electrons (e), up quarks (u) and down quarks (d). As in Ref. [21] we restrict ourselves to non-standard interactions with quarks, so that only ξ u and ξ d are relevant for neutrino propagation. A global rescaling of both ξ u and ξ d by a common factor can be reabsorbed into a rescaling of ε αβ , therefore only the direction in the (ξ u , ξ d ) plane is phenomenologically non-trivial.
In this framework, the evolution of the neutrino and antineutrino flavor state during propagation is governed by the Hamiltonian: where H vac is the vacuum part which in the flavor basis (ν e , ν µ , ν τ ) reads Here U vac denotes the three-lepton mixing matrix in vacuum [1,53,54] which we parametrize as [24], Explicitly: where c ij ≡ cos θ ij and s ij ≡ sin θ ij . If all possible operators in Eq. (2.1) are added to the SM Lagrangian the matter part H mat is a function of the number densities for the fermions in the matter N f (x) along the trajectory: where the "+1" term in the ee entry accounts for the standard contribution, and with describes the non-standard part, which we have written in terms of the effective couplings of protons (p) and neutrons (n) after taking into account that matter neutrality implies N e (x) = N p (x). So the phenomenological framework adopted here is characterized by nine matter parameters: eight related to the matrix ε αβ (after neglecting the trace which is irrelevant for neutrino oscillations) plus the direction in the (ξ p , ξ n ) plane.

NSI-Mass-ordering degeneracy
When considering the neutrino evolution in the most general matter potential described above one finds an intrinsic degeneracy as a consequence of CPT [19,20,23,24]. In brief CPT symmetry implies that the neutrino evolution is invariant if the Hamiltonian H ν = H vac + H mat is transformed as H ν → −(H ν ) * . This requires a simultaneous transformation of both the vacuum and the matter terms. The transformation of H vac implies ∆m 2 31 → −∆m 2 31 + ∆m 2 21 = −∆m 2 32 , θ 12 → π/2 − θ 12 , which does not spoil the commonly assumed restrictions on the range of the vacuum parameters (∆m 2 21 > 0 and 0 ≤ θ ij ≤ π/2). Eq. (2.9) involves a change in the octant of θ 12 as well as a change in the neutrino mass ordering (i.e., the sign of ∆m 2 31 ), which is why it has been called "generalized mass ordering degeneracy" in Ref. [24]. In the SM, this degeneracy is broken by the SM matter potential and that is what allows for the determination of the octant of θ 12 in solar neutrino experiments and of the ordering of the states in LBL experiments. But once the NSI are included the mass order degeneracy can be recovered if together with Eq.(2.9) one also exchanges see Refs. [20,23,24]. As seen in Eq. (2.7) the matrix E αβ (x) depends on the chemical composition of the medium, which may vary along the neutrino trajectory, so that in general the condition in Eq. (2.10) is fulfilled only in an approximate way. It becomes exact in several cases, for example if the couplings to neutrons cancel (ξ n = 0) so that the position-dependent term in Eq. (2.7) disappears. More interestingly the mass ordering degeneracy is also recovered when the neutron/proton ratio, Y n (x), is constant along the entire neutrino propagation path. This is the case for long-baseline experiments, either in accelerators or reactors, as in the Earth the neutron/proton ratio Y n (x) which characterize the matter chemical composition can be taken to be constant to very good approximation. It is also a reasonably good approximation for atmospheric neutrino experiments. The PREM model [55] fixes Y n = 1.012 in the mantle and Y n = 1.137 in the core, so that for atmospheric and LBL neutrino experiments one can set it to an average value Y ⊕ n = 1.051 all over the Earth. Hence for those experiments one gets E αβ (x) ≡ ε ⊕ αβ with: In other words, within this approximation the analysis of atmospheric and LBL neutrinos holds for any combination of NSI with up and down quarks (and also with electrons if it had been considered) and it can be performed in terms of the effective NSI couplings ε ⊕ αβ , which play the role of phenomenological parameters. In particular, the best-fit value and allowed ranges of ε ⊕ αβ are independent of the (ξ p , ξ n ) couplings.

New sources of CP violation
Another obvious consequence of the effect of NSI in the matter potential is the appearance of new sources of CP violation. As seen in Eqs. (2.5) and (2.6) the Hamiltonian contains now four possible physical phases which in that parametrization are the Dirac CP-phase in U vac and the three arguments of the off-diagonal elements in H mat . We emphasize that only one of the phases can be assigned to the vacuum part of the Hamiltonian and only one can be assigned to the matter Hamiltonian in a basis independent form. The other two phases are relative phases between both pieces of the Hamiltonian. This is more transparent if one uses the parametrization of Ref. [19] which for the relevant part of the NSI matrix in Eq. (2.1) reads: and A is a real number. The advantage of this formalism is that, when propagating in a medium with constant chemical composition (for example, the Earth with it is possible to choose A in such a way that the matrix H mat mimics the structure of the vacuum term in Eq. (2.3): (2.14) Written in this form, it is explicit that the two phases α 1 and α 2 included in Q rel are not a feature of neutrino-matter interactions, but rather a relative feature between the vacuum and matter terms.
To further illustrate this point we can introduce four basis invariants to characterize leptonic CP violation in the presence of non-standard neutrino interactions as described in the appendix A. They are built as products of hermitian matrices, concretely the charged lepton mass-squared matrix S = M M † , the neutrino mass-squared matrix S ν = M ν M ν † , and the ε matrix. One of them can be chosen as to contain only parameters which appear in the neutrino evolution in vacuum, so it gives a measure of leptonic CP violation in neutrino oscillations in vacuum: where the third equality shows its form in the charged lepton mass basis, with v(m e , m µ , In writing this expression we have kept U vac = R 23 (θ 23 , δ 23 )· R 13 (θ 13 , δ 13 ) · R 12 (θ 12 , δ 12 ) to explicitly show the physical phase appearing in this invariant δ CP ≡ δ 12 −δ 13 +δ 23 and which, for convenience, in our parametrization in Eq. (2.5) we have kept as attached to the 12 rotation by setting δ 13 = δ 23 = 0. Equation (2.15) corresponds to the usual Jarlskog invariant for the leptonic sector.
Unlike for the case of the invariants in Eqs. (2.15) and (2.16), there is not a clear physical set up which could single out the contribution from Eq. (2.17) and Eq. (2.18) (or any combination of those) to a leptonic CP violating observable. From their explicit expressions one also finds that setting the α i phases in Q rel as well as δ NS to zero is not enough to cancel those "matter-vacuum interference" CP invariants. On the contrary, setting those phases to zero one reintroduces a dependence on the phase convention for the phases in the vacuum mixing matrix.
Admittedly the discussion above is only academic for the quantification of the effects induced by the NSI matter potential on neutrino propagation, because the relevant probabilities cannot be expressed in any practical form in terms of these basis invariants and one is forced to work in some specific parametrization. What these basis invariants clearly illustrate is that in order to study the possible effects (in experiments performed in matter) of NSI on the determination of the phase which parametrizes CP violation in vacuum without introducing an artificial basis dependence, one needs to include in the analysis the most general complex NSI matter potential containing all the three additional arbitrary phases.

Analysis Framework
As mentioned in the introduction this work builds upon the results of our recent comprehensive global fit performed in the framework of three-flavor oscillations plus NSI with quarks [21] to which we refer for the detailed description of methodology and data included. In principle, for each choice of the (ξ p , ξ n ) couplings the analysis depends on six oscillations parameters plus eight NSI parameters, of which five are real and three are phases. To keep the fit manageable in Ref. [21] only real NSI were considered and ∆m 2 21 effects were neglected in the fit of atmospheric and long-baseline experiments. This rendered such analysis independent of the CP phase in the leptonic mixing matrix and of the ordering of the states.
In this work we extend our previous analysis to include the effect of the four CP-phases in the Hamiltonian as well as the ∆m 2 21 effects, in particular where they are most relevant which is the fit of the LBL experiments. In order to do so while still maintaining the running time under control, we split the global χ 2 in a part containing the data from the long-baseline experiments MINOS, T2K and NOvA (accounting for both appearance and disappearance data in neutrinos and antineutrino modes), for which both the extra phases and the interference between ∆m 2 21 and ∆m 2 31 are properly included, and a part containing CP-conserving observables where the complex phases can be safely neglected and are therefore implemented as described in Ref. [21]. In what follows we label as "OTH" (short for "others") these non-LBL observables which include the results from solar neutrino experiments [56][57][58][59][60][61][62][63][64][65][66], from the KamLAND reactor experiment [67], from medium-baseline (MBL) reactor experiments [68][69][70], from atmospheric neutrinos collected by IceCube [71] and its sub-detector DeepCore [72], and from our analysis of Super-Kamiokande (SK) atmospheric data [73]. 1 Also, in order to quantify the impact on the fit of the data on coherent neutrino-nucleus scattering collected by the COHERENT experiment [75] we will show the results after including COHERENT as part of OTH as well. It should be noted, however, that COHERENT bounds are only applicable for models where the mediator responsible for the NSI is heavier than about 10 MeV, as discussed in Ref. [21].
Schematically, if we denote by w the five real oscillation parameters (i.e., the two mass differences and the three mixing angles), by η the direction in the (ξ p , ξ n ) plane, by |ε αβ | the five real components of the neutrino part of the NSI parameters (two differences of the three diagonal entries of ε αβ , as well as the modulus of the three non-diagonal entries 2 ), and by φ αβ the three phases of the non-diagonal entries of ε αβ , we split the global χ 2 for 1 As in Ref. [21] we include here our analysis of SK atmospheric data in the form of the "classical" samples of e-like and µ-like events (70 energy and zenith angle bins). As discussed in Ref. [4] with such analysis in the framework of standard 3-nu oscillations we cannot reproduce the sensitivity to subdominant effects associated with the mass ordering and δCP found by SK in their analysis of more dedicated samples [74]. For that reason we include SK atmospheric data but only as part of the "OTH" group.
2 More precisely, in our analysis of OTH experiments we consider both the modulus and the sign of the non-diagonal ε αβ , i.e., we explicitly account for the all the CP-conserving values of the three phases: φ αβ = 0, π. However, in the construction of χ 2 OTH these signs are marginalized, so that only the modulus |ε αβ | is correlated with χ 2 LBL .
To make the analysis in such large parameter space treatable, we introduce a series of simplifications. First, we notice that in MBL reactor experiments the baseline is short enough to safely neglect the effects of the matter potential, so that we have: Next, we notice that in LBL experiments the sensitivity to θ 12 , ∆m 2 21 and θ 13 is marginal compared to solar and reactor experiments; hence, in χ 2 LBL we can safely fix θ 12 , θ 13 and ∆m 2 21 to their best fit value as determined by the experiments included in χ 2 OTH . However, in doing so we must notice that, within the approximations used in the construction of χ 2 OTH ( w, |ε αβ |, η), there still remains the effect associated to the NSI/mass-ordering degeneracy which leads to the appearance of a new solution in the solar sector with a mixing angle θ 12 in the second octant, the so-called LMA-Dark (LMA-D) [76] solution. Although LMA-D is not totally degenerate with LMA, due to the variation of the matter chemical composition along the path travelled by solar neutrinos, it still provides a good fit to the data for a wide range of quark couplings, as found in Ref. [21]. Concretely, after marginalization over η we get that the parameter region containing the LMA-D solution lies at Therefore, when marginalizing over θ 12 we consider two distinct parts of the parameter space, labelled by the tag "REG": one with θ 12 < 45 • , which we denote as REG = LIGHT, and one with θ 12 > 45 • , which we denote by REG = DARK. Correspondingly, the fixed value of θ 12 used in the construction of χ 2 LBL,REG is the best fit value within either the LMA or the LMA-D region: sin 2θlight 12 = 0.31 or sin 2θdark 12 = 0.69, respectively. The best fit values for the other two oscillation parameters fixed in χ 2 LBL,REG are the same for LMA and LMA-D: ∆m 2 21 = 7.4 × 10 −5 eV 2 and sin 2θ 13 = 0.0225. Further simplification arises from the fact that for LBL experiments the dependence on the NSI neutrino and quark couplings enters only via the effective Earth-matter NSI combinations ε ⊕ αβ defined in Eq. (2.11). It is therefore convenient to project also χ 2 OTH over these combinations, and to marginalize it with respect to η. In addition we neglect the small correlations introduced by the common dependence of the atmospheric experiments in χ 2 OTH and the LBL experiments on ∆m 2 31 and θ 23 , and we also marginalize the atmospheric part of χ 2 OTH over these two parameters. This means that in our results we do not account for the information on ∆m 2 31 and θ 23 arising from atmospheric experiments, however we keep the information on ∆m 2 31 from MBL reactor experiments. With all this, we can define a function χ 2 OTH,REG depending on six parameters:

Results
In order to quantify the effect of the matter NSI on the present oscillation parameter determination we have performed a set of 24 different analysis in the eleven-dimensional parameter space. Each analysis corresponds to a different combination of observables. The results of the long-baseline experiment MINOS are always included in all the cases, so for convenience in what follows we define χ 2 OTHM ≡ χ 2 OTH + χ 2 MINOS . To this we add χ 2 LBL with LBL = T2K, NOvA, and T2K+NOvA, and for each of these three combinations we present the results with and without the bounds of COHERENT. In addition, we perform the analysis in four distinctive parts of the parameter space: the solar octant "REG" being LIGHT or DARK, and the mass ordering being normal (NO) or inverted (IO). For illustration we show in Figs. 1 and Fig. 2 all the possible one-dimensional and two-dimensional projections of the eleven-dimensional parameter space accounting for the new CP violating phases, parametrized as φ αβ ≡ arg(ε αβ ) = arg(ε ⊕ αβ ). In both figures we show the regions for the GLOBAL analysis including both T2K and NOvA results and also accounting for the COHERENT bounds. In Fig. 1 we present the results for the LIGHT sector and Normal Ordering, while in Fig. 2 we give the regions corresponding to the DARK sector and Inverted Ordering; in both cases the allowed regions are defined with respect to the local minimum of each solution. From these figures we can see that, with the exception of the required large value of ε ⊕ ee − ε ⊕ µµ in the DARK solution, there is no statistically significant feature for the ε ⊕ αβ parameters other than their bounded absolute values, nor there is any meaningful information on the φ αβ phases. The most prominent non-trivial feature is the preference for a non-zero value of ε ⊕ eµ at a ∆χ 2 ∼ 2 level, associated with a φ eµ phase centered at the CP-conserving values π (0) for the LIGHT (DARK) solution. More on this below.
In order to quantify the effect of the matter NSI on the present determination of δ CP and the mass ordering we plot in Fig. 3 the one-dimensional χ 2 (δ CP ) function obtained from the above χ 2 GLOB,REG after marginalizing over the ten undisplayed parameters. In the left, central and right panels we focus on the GLOBAL analysis including T2K, NOvA, and T2K+NOvA respectively. The upper (lower) panels corresponds to the results without (with) the inclusion of the bounds from COHERENT. In each panel we plot the curves obtained marginalizing separately in NO (red curves) and IO (blue curves) and within the REG = LIGHT (full lines) and REG = DARK (dashed) regions. For the sake of comparison we also plot the corresponding χ 2 (δ CP ) from the 3ν oscillation analysis with the SM matter potential (labeled "NuFIT" in the figure).
For what concerns the analysis which includes T2K but not NOvA, i.e., the left panels in Fig. 3, we find that: The panels show the two-dimensional projections of the allowed parameter space after marginalization with respect to the undisplayed parameters. The different contours correspond to the allowed regions at 1σ, 2σ and 3σ for 2 degrees of freedom. Note that as atmospheric mass-squared splitting we use ∆m 2 3 = ∆m 2 31 for NO. Also shown are the one-dimensional projections as a function of each parameter. For comparison we show as dotted lines the corresponding one-dimensional dependence for the same analysis assuming only standard 3ν oscillation (i.e., setting all the NSI parameters to zero).
• The statistical significance of the hint for a non-zero δ CP in T2K is robust under the inclusion of the NSI-induced matter potential for the most favored solution (i.e., LIGHT-NO), as well as for LIGHT-IO. This statement holds irrespective of the inclusion of the bounds from COHERENT.
• The ∆χ 2 for DARK solutions exhibits the expected inversion of the ordering as well as the δ CP → π − δ CP transformation when compared with the LIGHT ones. This is a consequence of the NSI-mass-ordering degeneracy discussed in Eqs. (2.9) and (2.10).
• We notice that ∆χ 2 min,OTHM+T2K,DARK,IO = ∆χ 2 min,OTHM+T2K,LIGHT,NO because of the breaking of the NSI-mass-ordering degeneracy in the analysis of solar experiments as a consequence of the sizeable variation of the chemical composition of the matter crossed by solar neutrinos along their path. However as we see in the upper left panel This suggests that the DARK-IO solution can provide a perfect fit to T2K data, i.e., there is an almost total loss of sensitivity to the ordering in T2K. Indeed what the inequality above shows is that within the allowed DARK parameter space it is possible to find areas where the fit to T2K-only data for IO are slightly better than the fit for NO in the LIGHT sector (and than NO oscillations without NSI). For the same reason we also find that because there are DARK solutions with NO which give a "less bad" fit than the degenerate of the LIGHT-IO minimum. For example, in DARK-NO we find that the best fit value for ∆m 2 31 can be slightly larger than the best-fit |∆m 2 32 | in LIGHT-IO, which leads to a slightly better agreement with the results on ∆m 2 31 from MBL reactors. These solutions, however, involve large NSI parameters, in particular ε eµ and ε eτ , which are disfavored by COHERENT. Consequently in the lower left panel we read an slightly higher ∆χ 2 min,OTHM+T2K,DARK,IO 3. and a very similar difference in the quality of the fits between the orderings in LIGHT and DARK regions (with different sign of course), so χ 2 min,OTHM+T2K,DARK,NO − χ 2 min,OTHM+T2K,DARK,IO

5.1.
• For the same reason, without including COHERENT, the statistical significance of the hint of CP violation in T2K is reduced for the DARK solutions with respect to the LIGHT ones. We find that CP conservation (CPC), that is, a fit with all phases either zero or s π, lies at • However we still find that even without COHERENT  So the CL for CPC as naively read from the curves of δ CP still holds, or what is the same, there is no leptonic CP violation "hidden" when there is no CP violation from δ CP .
For the global combination including NOvA without T2K (central panels in Fig. 3) we notice that: • The sensitivity to δ CP diminishes with respect to that of the oscillation only analysis both in the LIGHT and DARK sectors. The constraints from COHERENT have a marginal impact on this.
• Within the DARK sector, IO is the best solution as expected from the NSI/massordering degeneracy, but it is still disfavored at ∆χ 2 min,OTHM+NOvA,DARK,IO ∼ 3 because of SOLAR+KamLAND, Eq. (3.3). By chance, this happens to be of the same order of the difavoring of IO in the pure oscillation analysis, ∆χ 2 min,OTHM+NOvA,OSC ∼ 3.5 (although the physical effect responsible for this is totally different).
In the global analysis including both T2K and NOvA (right panels in Fig. 3) we find qualitatively similar conclusions than for the analysis without NOvA, albeit with a slight washout of the statistical significance for both δ CP and the disfavoring of IO due to the tensions between T2K and NOvA. Such washout is already present in the oscillation-only analysis and within the LIGHT sector it is only mildly affected by the inclusion of NSI. However one also observes that in the favored solution, LIGHT-NO, maximal δ CP = 3π/2 is more allowed than without NSI. This happens because, as mentioned above, in NOvA the presence of NSI induces a loss of sensitivity on δ CP , so in the global analysis with both T2K and NOvA the behavior observed in T2K dominates.
In the global analysis there remains, still, the DARK-IO solution at without (with) including COHERENT. The status of the non-maximality and octant determination for θ 23 is displayed in Fig. 4 where we show the one-dimensional χ 2 (sin 2 θ 23 ) obtained from χ 2 GLOB,REG including both T2K and NOvA and after marginalizing over all the undisplayed parameters (so these are the corresponding projections to the left panels in Fig. 3). In the left (right) panel of Fig. 4 we show the results without (with) COHERENT. As seen in the figure the global analysis  including NSI for both orderings, for both LIGHT and DARK sectors, and irrespective of the inclusion of COHERENT, still disfavors the maximal θ 23 = π/4 at a CL ∼ 2σ ÷ 2.5σ. The main effect of the generalized NSI matter potential on θ 23 is the relative improvement of the CL for the first octant, which, compared to the best fit in the second octant is now at less than 2σ. For example for the two most favored solutions it is now at for the analysis without (with) COHERENT. This is so because the disfavoring of the first octant in the global oscillation-only analysis is driven by the excess of appearance events in NOvA. These events can now be fitted better with θ 23 in the first octant when including a non-zero ε eµ to enhance the ν µ → ν e flavor transition probability. For this reason, as seen in the left panel in Fig. 4 in the case of DARK-NO without COHERENT, the first octant becomes favored. But, as mentioned above, these solutions required large NSI which are disfavored by COHERENT and therefore once this is included in the analysis, the second octant becomes the best fit in all regions of parameter space.

Summary
In this work we have extended the analysis in Ref. [21] to account for the effect of NSI affecting neutrino propagation in matter on the observables sensitive to leptonic CP violation and to the mass ordering. We have quantified the robustness of the present hints for these effects in the presence of NSI as large as allowed by the global oscillation analysis itself. We conclude that the CL for the preference for a CKM-like CP phase close to 3π/2 in T2K, which is the one that drives the preference in the global analysis, remains valid even when including all other phases in the extended scenario. On the contrary the preference for NO in LBL experiments is totally lost when including NSI as large as allowed by the global analysis because of the intrinsic NSI/mass-ordering degeneracy in the Hamiltonian which implies the existence of an equally good fit to LBL results with IO and reversed octant of θ 12 and δ CP → π − δ CP (so in this solution the favored δ CP is also close to 3π/2). In the global analysis the only relevant breaking of this degeneracy comes from the composition dependence of the matter potential in the Sun which disfavors the associated LMA-D with CL below 2σ. Finally, we have also studied the effect of NSI in the status of the non-maximality and octant determination for θ 23 and find that for both orderings, for both LIGHT and DARK sectors, and irrespective of the inclusion of COHERENT, maximal θ 23 = π/4 is still disfavoured in the global fit at a CL ∼ 2σ ÷ 2.5σ. Including NSI, however, results into a decrease of the preference for the second octant with respect to the first octant to less than 2 σ.

A Invariants for leptonic CP violation with NSI
In this appendix we derive a set of basis and rephasing invariants that characterize leptonic CP violation in the presence of non-standard neutrino interactions. We follow the methodology introduced in Refs. [77,78] for generalizing the construction of such invariants in quark sectors first introduced for three generations in [79,80]. Notice that we are working with Dirac neutrinos which is all it is needed when interested in CP violation in neutrino oscillations.