Biprobability approach to CP phase degeneracy from non-standard neutrino interactions

Non-standard interactions (NSI) between neutrinos and matter at long-baseline experiments could make determination of the CP-violating phase $\delta_{13}$ ambiguous due to interference with additional complex phases. Such degeneracies are often studied in the context of specific experiments and a few parameter choices, leaving it unclear how to extract a general understanding of when two sets of parameters may be degenerate or how different types of experiments in principle combine to lift such a degeneracy. This work complements detailed simulations of individual experiments by showing how underlying parameters relate to degeneracies as represented on a biprobability plot. We show how a range of energies near the oscillation maximum $\Delta_{31} = \pi/2$ separates some degenerate probabilities along the CP-conserving direction of biprobability space according to $\delta_{+} \equiv \delta_{13} + \delta_{e\tau}$, while near $\Delta_{31} = 3\pi/2$ degenerate probabilities are separated along the CP-violating direction according to $\delta_{e\tau}$. We apply this to the experimental hints that suggest $\delta_{13} \sim -\pi/2$ to see that this could also be consistent with $\delta_{13}, \delta_{e\tau} = 0$ or $\pi$. The baseline and energy range characteristic of DUNE provides some resolution, but a further improvement comes from beams a few degrees off-axis at $\gtrsim 1000$ km baselines, including some proposed sites for T2HKK.


I. INTRODUCTION
An important goal of the current and upcoming generation of long-baseline neutrino oscillation experiments is to measure the phase δ 13 in the mixing matrix, which violates CP symmetry if its value is not 0 or π. Determining the existence and nature of leptonic CP violation is an important step in understanding the mixing properties in the standard 3-neutrino picture, and could have further implications for early universe leptogenesis [1,2]. However, neutrinos could experience beyond-Standard Model interactions with the matter in between source and detector. This possibility is often studied in a model-independent way as an addition to the matter potential parametrized by a set of non-standard interaction (NSI) parameters [3][4][5][6][7][8]. These NSI parameters include new CP-violating phases that contribute to oscillation probabilities along with δ 13 , leading to potential ambiguity in determining the true underlying parameters.
It is possible, for instance, that for a given choice of underlying parameters and experimental setup (neutrino energy and baseline length), there could be significant CP violation present in the model but not apparent in the data. In other words, nature chooses parameters which then imply oscillation probabilities P (ν µ → ν e ) and P (ν µ → ν e ), while we measure events to obtain these probabilities which we hope will tell us about the parameters. But different underlying parameters can give the same or similar oscillation probabilities, leading to a parameter degeneracy.
Computing the spectrum of events as a function of neutrino energy is well-defined, but working backward to interpret the result in terms of underlying parameters can be unintuitive. Furthermore, results are often presented in terms of sensitivity plots for particular experimental setups, making it more difficult to discern what is a feature of the basic oscillation parameters and what is related to specifics of the experiment. While analytic expressions are available, it can be difficult to tell at a glance what may be the effect of changing a given parameter. Therefore, it's desirable to understand how to concisely organize the information in a way that promotes easy correspondence between oscillation probabilities and favored parameters or models. One method that has been useful for representing other degeneracies, and the situations where they are lifted, is to plot oscillation probability for antineutrinos, P ≡ Prob(ν µ → ν e ), versus oscillation probability for neutrinos, P ≡ Prob(ν µ → ν e ), for values of δ 13 ranging from 0 → 2π [11]. This gives a curve that traces out an ellipse in biprobability space, and is a useful way to visualize degeneracies for a range of parameter values. An example of such for neutrino oscillations with standard interactions is shown by the solid curve in Fig. 1a. In Fig. 1a it is evident that a given P , P could correspond to widely separated values of δ eµ , and even the no-NSI case can be degenerate with them. Fig. 1b shows the range of possibilities by plotting dotted curves for values of δ eµ going from 0 to 2π in steps of π/4. Both plots are for L = 1300 km and correspond to ∆ 31 = π/2, and Fig. 1a shows only the normal mass hierarchy. Parameters and numerical methods are discussed in Sec. II.
The biprobability plot also provides a visualization of parameter degeneracies and their complexity. This is apparent from the dashed curves in Fig. 1a, which show the result when | eµ | = 0.05, for a few values of δ eµ . 1 There are evidently some values of the neutrino and antineutrino probabilities that do not lead back to a unique choice of the underlying parameters. Fig. 1b shows that this is a problem in general: the solid curves are for the normal and inverted mass orderings in the Standard Oscillation case (no NSI), and the dotted curves show biprobability ellipses for | eµ | = 0.05 and δ eµ varying from 0 to 2π in steps of π/4. The point emphasized by Fig. 1 but not unique to its parameter choices is that measuring P and P is not generally enough to uniquely distinguish the underlying parameters, in particular to distinguish δ 13 from δ eµ or δ eτ .
Looking at different baseline lengths and/or neutrino energies helps resolve such degeneracies, and some previous work has incorporated biprobability plots into studies of NSI. The ability of a combination of 3000 km and 7000 km baseline experiments to distinguish CP violation due to NSI was considered in [35]. This work also showed how biprobability plots are affected by NSI, but did not quantify degeneracies in the manner that the present paper does. Biprobability plots have also been used to represent degeneracies in the presence of eτ at NOνA, with examples of degenerate probabilities at NOνA that are separated with the broader energy spectrum at DUNE [15], and to help represent an analytic result, showing that any "apparent" δ 13 and hierarchy in the standard case could correspond to δ 13 = 0 for some values of eµ and δ eµ (and similar for eτ ) [23]. Implications of nonzero eτ , the effect on biprobability plots, and the ability of DUNE to distinguish some values of δ 13 and δ eτ was also considered in [27,36]. However, as mentioned above, it is still desirable to have a general picture of how oscillation parameters influence these degeneracies, separate from detailed experimental scenarios. Some work has used biprobability ellipses to represent eτ degeneracies relevant to an apparent δ 13 ∼ −π/2 [37]. This paper shows how δ 13 and NSI parameters influence the biprobability space representation of ν µ → ν e and ν µ → ν e oscillations. We apply this insight to study degeneracies between δ 13 and NSI phases δ eµ , δ eτ . In the absence of NSI, there is an approximate degeneracy between values of δ 13 with the same sin δ 13 . In the biprobability context, looking at energies away from the oscillation maximum (at ∆ 31 in terminology to be introduced in Sec. II) helps to break degeneracies, largely by increasing the separation between two points along the "CP-conserving" (or P + P ) direction on the biprobability plot. This carries over to the case of nonzero NSI, where two points of the same sin δ 13 and the same δ + ≡ δ 13 + δ eµ or δ 13 + δ eτ are approximately degenerate near the ∆ 31 = π/2 maximum. As the energy varies from this maximum, δ + (rather than δ 13 or δ eτ ) plays the dominant role in separating points. On the other hand, near ∆ 31 = 3π/2 it is δ eτ or δ eµ that is dominant, and this tends to separate degenerate points along the CP-violating (or P − ) direction. We examine the approximate degeneracy between δ 13 = −π/2 (i.e. δ 13 = −π/2 with no NSI, apparently maximal CP violation) and the CP-conserving situation δ 13 = 0, π, δeτ = 0, π for | eτ | = 0.02. Using the above results, we see that examining the energy spectrum around ∆ 31 = π/2 (relevant to DUNE) provides some separation between these points, while ∆ 31 ≈ 3π/2 (relevant to T2HKK) further improves upon this. The purpose of this paper is not to make detailed predictions about whether specific experiments will be able to resolve some parameters at a given statistical significance. In contrast, we seek to understand at a general level how the parameters relevant to the next decade of experiments work together to determine degeneracies and their breaking.
Before we examine degeneracies, Sec. II will describe the numerical and analytic methods used in this paper. Then Sec. III will apply these to understand degeneracies and how they are represented on the biprobability plot and Sec. IV will show the role of the energy dependence in lifting the degeneracy in biprobability space. Sec. V applies these results to the interesting question of which degeneracies could be present if δ 13 has apparently been measured as 3π/2, and how this degeneracy breaking can be represented in biprobability space. Sec. VI summarizes the results and outlines future work that will build on the results presented here.

II. NUMERICAL AND ANALYTIC COMPUTATION OF OSCILLATION PROBABILITIES
In this section, we describe the specific methods and parameter choices used to compute oscillation probabilities in this paper. In particular, we consider NSI affecting propagation (i.e. matter potential) but not affecting interactions in the detector. In this case, NSI appear as parameters ij in the flavor-basis interaction Hamiltonian 2 : where the off-diagonal elements have complex phases defined by ij ≡ | ij | exp (iδ ij ), the parameter encoding the matter potential is and the mixing matrix in terms of s ij ≡ sin θ ij , c ij ≡ cos θ ij is where L is the baseline length, E is the neutrino energy, N e is the number density of electrons, ρ is the density of matter that neutrinos propagate through, and G F is the Fermi constant. Table I gives the values of various neutrino oscillation parameters used for the calculations in this paper. These are based on the 2016 Particle Data Group (PDG) Review of Particle Physics [38]. In particular, one dimensionless combination that will come up often in the context of perturbation expansion is 3 Since the purpose of this paper is to gain a general understanding of the behavior of biprobability plots within reasonable parameter ranges, not rigorous implementation experimental constraints or simulation of specific experiments or detectors, we examine a range of NSI parameters inspired by the different bounds listed in Table I. These are completely model-independent bounds ("NSI Range A") and more restrictive but model-dependent bounds ("NSI Range B") listed in section VI.A of [39]. In this work, we will examine eµ , eτ ∼ O(0.01) to O(0.1). The reason for the lower limit is that for small enough αβ , CP phase degeneracies will not be significant, a point we will return to in Sec. III. In matter, P and P may generally be found by numerical integration: the time-evolution (equivalently, distance-evolution) of amplitudes follows from the Hamiltonian. Oscillation probabilities shown in this paper are the result of such a numerical   [38]. The NSI parameter ranges are based on model-independent ("NSI Range A") and more model-dependent ("NSI Range B") bounds as described in section VI.A of [39].
integration using the Hamiltonian Eq. (1), with matter density ρ = 3.0 g/cm 3 . Realistic variations in ρ will change the results quantitatively but not qualitatively. We will restrict our attention to the electron neutrino (antineutrino) appearance processes ν µ → ν e (ν µ → ν e ) studied at T2K, NOνA, and DUNE. These experiments aim for a combination of baseline length and neutrino energy giving ∆ 31 ≡ ∆m 2 31 L/(4E) ∼ ±0.003 L[km] / E[GeV] = π/2, near the oscillation maximum. Results in this paper use L = 1300 km when necessary in order to give a sense of the size of effect that may be relevant to DUNE, but this should not be interpreted as a precise statement of what will be measured at the DUNE experiment. Relating our results to the capabilities of specific experiments is left to future work.
As an alternative to the numerical approach, it is often convenient to treat the parameters ∆m 2 21 /∆m 2 31 , sin(θ 13 ), and | αβ | as small parameters of the same order [9,40,41]. Then oscillation probabilities may be represented by approximate analytic solutions. While not exact, this approach is still helpful in illuminating trends seen in numerical results, as we will see in the remaining sections. Before we get to the full perturbative expressions, we will point out some important general features.
In vacuum, the Standard Model oscillation probabilities can be expressed in the form which is the parametric form for an ellipse in the P, P plane as δ 13 varies from 0 to 2π [11]. Using the perturbation expansion [40] in matter, the oscillation probabilities still take the form of Eq. (5), but with the coefficients N , C 13 , etc. altered by the matter effects. In the perturbation expansion with NSI present, the form of Eq. (5) remains [41], i.e. varying δ 13 still traces an ellipse in biprobability space. For later use in this paper, we will refer to the standard oscillation results as P 0 , P 0 and the effect of NSI as a ∆P eµ , ∆P eµ or ∆P eτ , ∆P eτ that would be added to P 0 , P 0 : and the same expression for eτ with ∆P eµ → ∆P eτ . It's interesting to note that the probabilities can also be expressed in a form that factors out sines and cosines of the "hidden sector" phase δ eµ rather than δ 13 : and the same can be done for δ eτ . While the form of Eq. (5) is an ellipse along which δ 13 continuously varies from 0 to 2π, Eq. (7) describes a biprobability ellipse along which δ eµ continuously varies for one constant value of δ 13 . Examples of these "hidden sector" ellipses from numerical integration are shown in Fig. 3a, superposed on the standard biprobability ellipse that varies δ 13 when eµ = 0. Previous works have also used this form [27,35,36,42], and in Sec. III we will find it particularly useful for representing degeneracies.
In the standard oscillation case, the perturbative expansion [9,40,41] gives probabilities where + (−) corresponds to P 0 (P 0 ) and J r ≡ c 12 s 12 c 2 13 s 13 c 23 s 23 (with c 12 ≡ cos(θ 12 ) etc.). With eµ = 0 but all other NSI parameters absent [41], the probabilities in Eq. (8) are adjusted by With eτ = 0 but all other NSI parameters absent [41], the probabilities are adjusted by where δ + ≡ δ 13 + δ eµ or δ 13 + δ eτ (11) To keep the notation simple, we use δ + in the case of either nonzero eµ or eτ . In this paper we will consider either one at a time to be nonzero, so from context the definition of δ + will be clear.

III. CP PHASE DEGENERACIES IN THE PRESENCE OF NSI
We will now use the methods reviewed in Sec. II to develop an understanding of degeneracies in the presence of NSI. In Sec. III A we will put the perturbative expressions into a more convenient form for the examination of CP Violation in biprobability space, and in Sec. III B we will use this along with numerical solution of probabilities in order to study degeneracies. Ref. [27] also used this combination of methods to study the effect of NSI as seen on biprobability plots and agrees with our results when there is overlap, but here we have a narrower scope and focus more on quantifying phase degeneracies and their breaking.

A. Applying the Perturbative Expressions
It will be convenient to note that P takes the schematic form of P = (CP even terms) + (CP odd terms) , so P = (CP even terms) − (CP odd terms) .
While previous works (e.g. [21]) have considered such a decomposition, here we will take this idea further and see that it leads to a useful result. In biprobability space, Eq. (12) suggests that the rotated coordinates are natural to look at in the context of CP violation, since We can see this explicitly using the perturbative expressions: Eq. (8) gives Eq. (9) gives and Eq. (10) gives Using the known values of some parameters as outlined in Table I Eq. (16) becomes and Eq. (17) becomes We focus on long-baseline experiments with ∆ 31 ∼ π/2 or 3π/2, so the ∆ 31 and 1/∆ 31 factors above will be ∼ O(1) while sin(2∆ 31 ) will be equal or close to zero. Sec. IV will take a closer look at the effect of the ∆ 31 -dependent terms as the energy varies. We are now in a better position to quantitatively study degeneracies using the results of Sec. III A and the "hidden sector" biprobability ellipses that vary δ eµ or δ eτ for some fixed value of δ 13 . In this section, we will see that these are useful because to a good approximation, at ∆ 31 ≈ π/2 the "center" is determined by sin(δ 13 ) and the dimensions (i.e. major and minor axes) are determined by eµ or eτ . In other words, hidden sector ellipses with the same magnitude NSI parameter and same δ 13 or π − δ 13 are approximately degenerate, and a given point on the ellipse only determines δ + . Previous works (see e.g. [15]) have used an approximation that helps give a sense of when points are degenerate: points separated by 0.01 on the biprobability plot won't be resolved. While this does not fully represent a particular experiment's ability to resolve parameters, our aim here is to study at a general level when degeneracies may exist or be broken, so this convenient rule of thumb is sufficient for our purposes. In Sec. IV, we will see that this degeneracy no longer holds for an energy-baseline combination away from ∆ 31 ≈ π/2. From Eq. (18) it is apparent that, in the absence of NSI, the probabilities in the "CP-violating direction" P − 0 are controlled by both sin δ 13 and matter effects, while probabilities in the "non-CP-violating direction" are controlled by cos δ 13 . Writing the probabilities in this form, we can see how matter effects help break the degeneracy between normal and inverted mass ordering. Under NH ↔ IH, ∆ 31 changes sign while the magnitude is only slightly affected. Since both terms multiplying aL/E change sign under ∆ 31 → −∆ 31 (∆ 31 and sin(2∆ 31 ) change sign, while sin 2 (∆ 31 ) does not), the effect at first order is just to reflect the ellipse across P − = 0. This role of the matter effect was noted in [11] and is evident in the solid curves (standard oscillation case) in Fig. 1b of this paper.
Another insight from the case of standard interactions relates to δ 13 . Taking as an example L = 1300 km and ρ = 3 g/cm 3 , Eq. (18) gives There is evidently a degeneracy between a given δ 13 and π − δ 13 , since the sine of each is the same. Fig. 2a shows this approximate degeneracy (note that the approximate result giving Eq. (18) does not exactly predict P and P , as they are not exactly degenerate, but it does accurately describe the general behavior). Now considering nonzero eµ or eτ , in Eq. (19) the coefficient of the δ + term is appreciably greater than the coefficient of the δ eµ term, for both ∆P + eµ and ∆P − eµ . The same is true in Eq. (20). This suggests that hidden sector ellipses (constant δ 13 , δ eµ varying from 0 to 2π) with δ 13 and π − δ 13 should approximately overlap. Example of biprobability plots where solid curves correspond to the standard ellipse that varies δ 13 , with | eµ | = 0 and dashed curves correspond to the "hidden sector ellipse" that varies δ eµ with fixed δ 13 , as discussed after Eq. (7). Fig. 3a shows (for NH) hidden sector ellipses for four values of δ 13 ; these can be thought of as "how a given point on the biprobability curve changes due to NSI." In particular, there is an approximate degeneracy between the δ 13 = 0 and δ 13 = π curves, which Fig. 3b focuses on, also giving specific values of δ eµ as discussed in the text.
It is useful to consider properties of the "hidden sector ellipses" of fixed δ 13 and varying δ eµ or δ eτ , such as their width or the location where they are centered. These may be estimated in a straightforward way based on Eq. (18), Eq. (19), and Eq. (20). The lowest-order terms have ∆P + eµ, eτ ∝ sin(δ + ), ∆P − eµ, eτ ∝ cos(δ + ), i.e. the form of an ellipse whose major (minor) axis is aligned with the P + (P − ) axis. As δ eτ varies with all other parameters fixed, points on the hidden sector ellipse vary between: for nonzero eµ and for nonzero eτ , where the numerical values are obtained for L = 1300 km, ρ = 3 g/cm 3 and ∆ 31 = π/2. The approximate center of the ellipse (with ∆P ± eµ, eτ = 0) is determined by Eq. (21). Therefore, at the leading order the width of the hidden sector ellipses is determined by eτ and the center is determined by sin(δ 13 ). So one can say at this level of approximation that hidden sector ellipses of the same sin δ 13 are centered at the same point, and with the same eµ or eτ they have approximately the same major and minor axis, so they overlap. It's interesting that the effect of NSI on probabilities depends so strongly on δ + , and we can gain some insight by noticing that the δ eµ or δ eτ terms in Eq. (16) or Eq. (17) are suppressed by a power of ∆m 2 21 /∆m 2 31 . In the ∆m 2 21 /∆m 2 31 → 0 limit the number of independent phases that control oscillation probabilities is reduced (see e.g. the phase reduction theorem in Section IV.C of [41]). In the absence of NSI, we can see from Eq. (15) that taking ∆m 2 21 /∆m 2 31 → 0 eliminates the effects of δ 13 , so there would be no intrinsic CP violation. With only one nonzero NSI element eµ or eτ , taking the limit ∆m 2 21 /∆m 2 31 → 0 would leave only one phase, δ + = δ 13 + δ eµ or δ 13 + δ eτ . Using the actual nonzero value of ∆m 2 21 , the individual phases δ eµ or δ eτ are not absent but suppressed by small but nonzero ∆m 2 21 /∆m 2 31 . Therefore, the result that hidden sector ellipses of constant δ 13 , π − δ 13 approximately overlap with a δ + degeneracy can be thought of as a consequence of the mass hierarchy. 4 The degeneracy depends on the magnitude | eµ | as well, as shown in Fig. 4. If eµ is "large enough" ( 0.01), there are approximately overlapping hidden sector ellipses. (This is consistent with a similar observation in [27].) For eµ ∼ 0.01, the "low P + " side of one approximately overlaps the "high P + " side of the other. For "small enough" eµ ( 0.01) the hidden sector ellipses are about as well-separated in the P + direction as the original P + 0 values. In this case the value of δ eµ does not significantly interfere with determination of δ 13 (this is expected: as eµ → 0 in Eq. (19), the NSI effects ∆P ± eµ go away). It is evident from Eq. (18), Eq. (19) and Eq. (20) that the qualitative behavior we have just described does not change under ∆ 31 → −∆ 31 , which well approximates the difference between normal and inverted hierarchy. We have confirmed numerically that the same characteristics of degeneracies hold for the inverted hierarchy as suggested by the above reasoning. Furthermore, comparing the form of Eq. (19) with Eq. (20), we can see that the behavior just described for eµ = 0, eτ = 0 should carry over to the case where eµ = 0, eτ = 0. This is also supported by numerical solutions; an example is shown in Fig. 7a.
Thus, if one knows the presence of nonzero eµ or eτ , there is at first order a degeneracy where only δ + is determined. A further complication in the degeneracy picture is not knowing whether eµ or eτ is nonzero. We saw this in Fig. 1 as overlap between biprobability ellipses with and without NSI. More quantitatively, comparison of Fig. 2a with Fig. 6a (showing hidden sector ellipses for δ 13 = π/4 and 3π/4) shows that the point on the biprobability plot corresponding to eµ = 0 and δ 13 = π/2 may also correspond to eµ = 0.05 and either δ 13 = π/4, δ eµ ∼ 5π/4 (black curve) or δ 13 = 3π/4, δ eµ ∼ 7π/4 (blue curve). Other points on these hidden sector curves intersect the eµ = 0 curve near δ 13 = 0 or π -if nonzero eµ isn't suspected, this would simply look like there is little or no intrinsic CP violation in the leptonic sector. Returning to Eq. (19) and Eq. (20), we can see that the sin δ 13 and cos δ + terms compete to determine the position of a hidden sector ellipse in the P − direction, while cos δ 13 and sin δ + compete in the P + direction. As we will see in Sec. IV, the δ + dependent terms help resolve this degeneracy as energy is varied. 4 "Mass hierarchy" used in its original sense to mean difference in scales |∆m 2 21 /∆m 2 31 | 1, not the common practice of referring to mass ordering, i.e. sign of ∆m 2 31 .

IV. LIFTING PHASE DEGENERACIES
As mentioned in Sec. I, variation in oscillation probabilities with neutrino energy and baseline length can help break degeneracies, a situation we now study using the above results. For convenience and motivated by experiments like DUNE and T2HK, at a given baseline length L we define E 0 as the neutrino energy for which ∆ 31 = π/2. Then an arbitrary neutrino energy can be parametrized as so the parameter x ≡ (E − E 0 )/E 0 is the fractional difference in energy from ∆ 31 = π/2. To get a sense of the numbers involved, for a situation where ∆ 31 = π/2 at an energy of 2.5 GeV, the parameter x will vary between x = ±0.2 as energy varies from 2.0 GeV to 3.0 GeV. This is motivated by DUNE, which will have a wide band of energies. By comparison, a narrow band beam like NOνA is often well approximated as a single energy, as in [15] for example. With this in mind, we will first expand to lowest order in x about x = 0 (i.e. expand about E = E 0 ) terms in P and P that depend on ∆ 31 . Then we will compare this with the behavior near x = − 2 3 , motivated by the off-axis beams proposed for T2HKK.
Based on the definition Eq. (24), we can make the exact substitutions and expand to first order in x. Then Eq. (18) becomes and Eq. (20) becomes We can see from Eq. (27) that the δ 13 , π − δ 13 degeneracy is broken by increasing the separation in the P + direction, since sin(δ 13 ) = sin(π − δ 13 ) but cos(δ 13 ) = − cos(π − δ 13 ), so the cosine term in P + is responsible for separating these points based on δ 13 . For instance, δ 13 = 0, π gives cos δ 13 = ±1, so that x = 0 moves these in opposite directions along P + . Furthermore, since representative numbers give 5.4 × 10 −5 (1300) (3) (0.05) ≈ 0.011 (and correspondingly less for smaller or L < 1300 km), it is the standard δ 13 terms that dominate as the energy changes. Comparison with Fig. 2b shows the advantage and limitations of the perturbative solutions Eq. (27), Eq. (28), Eq. (29). The qualitative behavior agrees, with the δ 13 degeneracy is broken by stretching the ellipse in the P + direction, but as expected with a first-order perturbative expansion, there are clearly higher-order effects missing that prevent the precise description of this behavior.
It is also evident from the perturbative results that decreasing the energy, so that the fraction x is negative, enhances the effect of sin δ 13 , cos δ + , and the matter effect due to the ∼ (1 − x) terms, while increasing the energy by a fraction will diminish these terms. This means that the energy spectrum below E 0 should be most useful in breaking degeneracies; this is also seen numerically in Fig. 5 and Fig. 6. The influence of the x cos δ 13 term in P + 0 is evident, as x < 0 moves the δ 13 = 0 curve toward lower P + and the δ 13 = π curve to higher P + relative to each other. Fig. 5 shows for the normal hierarchy how NSI ellipses with constant δ 13 = 0 and δ 13 = π at ∆ 31 = π/2 (Fig. 5a) are affected as the energy varies: increased by x = +0.1 in Fig. 5a, and decreased by x = −0.1 in Fig. 5b and x = −0.2 in Fig. 5d). Fig. 6 shows the same result for δ 13 = π/4 and 3π/4. Fig. 7 shows a similar result for hidden sector ellipse degeneracy with eτ = 0.05. Fig. 7a has ∆ 31 = π/2 and the hidden sector ellipses that vary δ eτ for fixed δ 13 = π/4, 3π/4 are approximately degenerate. As in the eµ case, and as suggested by  Fig. 5a shows the "hidden sector" biprobability curves for the normal hierarchy with L = 1300 km and ∆ 31 = π/2. Fig. 5b shows the situation with energy increased from the peak by a fractional amount x = 0.1, while Fig. 5c and Fig. 5d show the energy reduced from the peak, they correspond to fractional changes of x = −0.1 and x = −0.2 respectively.
Eq. (29), variation in the energy works to lift the degeneracy. This is shown in Fig. 7b, where the energy has been shifted by x = −0.1.
Finally, we point out how variations in P ± relate back to the spectrum of probabilities (or number of events at the detector) as a function of neutrino energy. Neutrinos and antineutrinos having identical spectra would mean that for any energy, the location on the biprobability plot is along P + (i.e. on the P − = 0 axis). Increasing along the P + direction means that the height of both neutrino and antineutrino spectra are increased equally, while increasing (decreasing) along the P − direction means that for the given energy, the probabilities in the neutrino spectrum are increased (decreased) relative to the antineutrino spectrum.
The above reasoning implies that for x ∼ −2/3, i.e. ∆ 31 ∼ 3π/2, points would be well-separated in the P − direction based on δ eτ , with only a subdominant effect from δ + . In Sec. V we will put this to use for the study of the degeneracy at apparent δ 13 ∼ −π/2. Fig. 9d shows numerical results that support the analytic development here.
A relevant situation is a long-baseline experiment with detectors off-axis by a few degrees, such as the proposed T2HKK sites [43]. We can evaluate these proposed T2HKK sites using the approximate relation E peak ν ≈ (30 MeV) /θ for the peak This means that the Unjang and Minjuji sites can separate points further in biprobability space according to δ eτ , in contrast with separation based on δ + as we have seen so far. (Fig. 9d referenced above shows results for parameters relevant to the Unjang site.) It is worth noting that the overall neutrino flux decreases with increased off-axis angle. This is an important consideration in evaluating a specific experiment's ability to distinguish two parameter choices, but we leave it for future work as it is beyond the scope outlined in Sec. I.
The results of the previous sections may now be applied to gain insight into a case of current interest: the implication of probabilities that suggest δ 13 ≈ −π/2. While no values of δ 13 have been ruled out experimentally, recent results from the T2K and NOνA experiments favor δ 13 = 0, π. T2K gives δ 13 in the range 1.06π to 1.86π at 90% confidence level [33], while NOνA gives δ 13 in the range 0 to 0.12π or 0.91π to 2π at 68.3% confidence level, with best fit 1.21π [34]. As previously stated, this paper does not aim for precise analysis of any experiment, but we will take this as motivation to specifically point out the implications of our work for the P, P points centered around δ 13 = −π/2. Earlier results consistent with these have been discussed in [42], which shows a few examples of the biprobability curves that we have called "hidden sector" curves, corresponding to eτ = 0.3 and δ 13 = 0 or π. They concluded that this particular choice of parameters is consistent with T2K and NOνA results, i.e. there is a degeneracy of the type we have considered in this paper. This section presents results consistent with [42] for varying eτ and δ 13 , and also considers degeneracy breaking. Other work has considered the implications of apparent δ 13 = −π/2 for NSI and the possibility of sterile neutrino mixing [19].
Recalling our previous discussion of the approximate δ 13 , π − δ 13 degeneracy in the standard oscillation case, the "endpoints" of the narrow biprobability ellipse at δ 13 = ±π/2 are special cases because the approximate δ 13 , δ 13 − π degeneracy goes away. However, the possibility of eτ = 0 means the other degeneracies we've considered may still apply. Specifically, a hidden sector ellipse for some δ 13 = −π/2 may overlap with point where the standard oscillation case predicts δ 13 = −π/2. And as we've seen, this means that the π − δ 13 hidden sector ellipse would also approximately overlap at this point. The larger eτ is, the farther this δ 13 or π − δ 13 would be from the apparent value of δ 13 = −π/2. This is illustrated in Fig. 8a for eτ = 0.02, 0.05. On each hidden sector curve, we've marked one specific point that happens to be close to the δ 13 = −π/2 point of the standard standard oscillation case. In general, it's evident that larger eτ directly corresponds to larger deviation of δ 13 from the apparent value. In this case, all marked points have δ + ∼ 2π → 0, so that none of the δ + dependent terms in the perturbation expansion help us distinguish the points by varying energy. Furthermore, the δ 13 's we consider are all within π/4 of the original point in question (δ 13 = −π/2). (Larger eτ improves this situation.) This means the δ 13 dependent terms, which are the dominant terms in the perturbation expansion, are also restricted in their use.
Varying the energy still does help break the degeneracy, as seen in Fig. 8b where the energy is increased by x = +0.2, and in Fig. 8c where the energy is decreased by x = −0.2. Looking at energies below the E 0 where ∆ 31 = π/2 is still helpful in separating the points, but the close proximity of all of these points in parameter space of CP phases means that the first-order approximation we've used is less helpful. We leave the consideration of higher-order terms and their effect on this picture to future work. An interesting situation is shown in Fig. 9. A measurement that apparently gives maximal CP violation in the standard picture would be consistent with nonzero NSI and absence of CP violation in both the standard and hidden sectors. For apparent δ 13 = −π/2 at the baseline and energy of NOνA, the points for | eτ | = 0.2, δ 13 = π, δ eτ = π and | eτ | = 0.2, δ 13 = 0, δ eτ = 0 are essentially degenerate with the no-NSI δ 13 = −π/2 result. At the baseline and energy of DUNE (with x = 0), these points separate by an amount ∼ 0.01 in the P direction. However, looking at x = −0.2, we see that δ 13 , δ eτ = π is well separated from the other two points (which are not distinguished according to the ∆P 0.01 rule of thumb referenced in Sec. III B).
This can be understood using Eq. (18) and Eq. (20). At ∆ 31 = π/2, the expressions for P ± = P ± 0 + ∆P ± eτ approximately coincide for the above three parameter choices. Looking at the coefficients of x in the expressions will then reveal how varying the energy does or does not break this degeneracy. In this case, it's convenient to return to the un-rotated coordinate P = 1 √ 2 (P + + P − ) (using the definition Eq. (13)). Schematically, this expression looks like P = (· · · ) + 1 √ 2 x (· · · ). For δ 13 = −π/2 and eτ = 0 the term in parentheses multiplying x gives −0.058, while for δ 13 = 0, | eτ | = 0.2, δ eτ = 0 this term has the similar value −0.053. However, for δ 13 = π, | eτ | = 0.2, δ eτ = π this term gives −0.11, approximately twice the value for the other two. This is consistent with Fig. 9, where varying energy so that x < 0 splits the point corresponding to δ 13 , δ eτ = π away from the other two along the P direction. Taken together, the numerical results and the perturbation expansion illustrate the role of a broad energy spectrum in addition to different baseline lengths in breaking specific degeneracies.
far-separated along the P − direction from δ 13 = 0, δ eτ = 0. To reiterate the conclusion of Sec. IV B: when switching from ∆ 31 ∼ π/2 → 3π/2, the terms in P − switch from being dominated by δ + to δ eτ . We also reiterate that off-axis beams have a lower neutrino flux, so this further separation in biprobability space does not automatically mean easier resolution. In this paper we have focused on degeneracies in terms of the question "when do different parameters give similar predictions?" and note that this is a case where very different parameters (no CP violation vs. maximal CP violation) give the same prediction at some experiments, and quite different predictions at another.

VI. SUMMARY AND FUTURE WORK
Here we summarize the main features found in the body of this paper, and discuss future directions that aim to use this framework. In the absence of NSI, a given δ 13 leads to probabilities approximately degenerate with π − δ 13 when the baseline and neutrino energy are chosen to give ∆ 31 = π/2. For fixed baseline, looking at neutrino energies away from the E 0 giving ∆ 31 = π/2 works to break this degeneracy. On a biprobability plot, decreasing the energy (x < 0) from E 0 works best to increase separation in the P + direction, improving the possibility for resolving different parameters.
In the presence of nonzero NSI parameter eµ (or eτ ), the effect on the oscillation probabilities near ∆ 31 = π/2 is dominated by the sum δ + = δ 13 + δ eµ (or δ 13 + δ eτ ), with a smaller contribution from δ eµ (or δ eτ ) alone. In combination with the first result listed above there is an approximate degeneracy among pairs of "hidden sector" ellipses with δ 13 and π − δ 13 , as seen in Fig. 5 and Fig. 6 for eµ , and Fig. 7 for eτ . At leading order, these hidden sector ellipses have a center determined by sin(δ 13 ) and major/minor axis length determined by | eµ | or | eτ |. As in the no-NSI case, looking at neutrino energies below the E 0 corresponding to ∆ 31 = π/2 tends to increase separation between two previously-degenerate points in the P + direction of a biprobability plot. These results are consistent with previous work on degeneracies in the presence of nonzero eτ [23], but provide further insight into how individual parameters affect the degeneracy as presented in biprobability space.
For probabilities consistent with δ 13 ∼ −π/2 in the standard oscillation case, nonzero eτ can produce the same oscillation probabilities for certain values of δ 13 and δ eτ . The larger eτ is, the more δ 13 can vary away from the apparent value of −π/2. In particular, apparent δ 13 = −π/2 (maximal CP violation in the standard case) is also consistent with | eτ | = 0.2 and δ 13 , δ eτ = 0 or δ 13 , δ eτ = π (NSI with no CP violation). For experimental parameters of DUNE, these points are more easily distinguishable than at NOνA. In particular, both numerical results and the perturbative approach in this paper show that the δ 13 , δ eτ = π point becomes well-separated as energy varies, but the δ 13 , δ eτ = 0 point remains more difficult to distinguish from maximal CP violation in the standard scenario. This is consistent with results presented in [15,42], but specifically addresses the degeneracy between maximal, standard CP violation and NSI without CP violation while providing a clear relationship between the underlying parameters and their effect on degeneracies as energy is varied.
Furthermore, we have shown that baseline and energy combinations relevant for long-baseline off-axis beams like T2HKK may dramatically help split degenerate points along the P − direction, including points not distinguishable just from the energy spectrum at an experiment like DUNE. In the earlier terminology where x is the fractional change in energy from the ∆ 31 = π/2 maximum, such an off-axis experiment may provide x ≈ −2/3. In comparison with x ≈ 0, the relative importance of the δ + and the δ eτ terms in the expression for P − is reversed.
The perturbative solutions have provided insight that supports each of these observations, but the first-order solutions presented in [9,40,41] only take us so far. It may also be useful examine the next-order terms in [41], but this is probably not consistent as long the assumption that θ 13 is small remains. Removing this assumption leads to perturbative expressions that can adjust the previous results by ∼ 0.001 to 0.01 for the L and E that interest us, and this is the same order as some effects we've considered here [48]. This paper has considered approximate degeneracies and discussed how looking at varying energies tends to break degeneracies as viewed on a biprobability plot. This is a valid approximation for now because the difference between the black and blue curves in Fig. 7a, for example, is much finer than any current experiment can distinguish. To help give context to the implications of how far separated biprobablility points are, we have roughly thought about a 0.01 separation in biprobability space as being degenerate. However, we have specifically avoided precise statements that characterize exactly when two parameters are degenerate because the answer depends on further experimental details. A thorough examination of this for experiments such as DUNE and T2HKK is an important next step to build on this work.