Correlations and degeneracies among the NSI parameters with tunable beams at DUNE

The Deep Underground Neutrino Experiment (DUNE) is a leading experiment in neutrino physics which is presently under construction. DUNE aims to measure the yet unknown parameters in the three flavour oscillation scenario which includes discovery of leptonic CP violation, determination of the mass hierarchy and determination of the octant of $\theta_{23}$. Additionally, the ancillary goals of DUNE include probing the sub-dominant effects induced by new physics. A widely studied new physics scenario is that of nonstandard neutrino interactions (NSI) in propagation which impacts the oscillations of neutrinos. We consider some of the essential NSI parameters impacting the oscillation signals at DUNE and explore the space of NSI parameters as well as study their correlations among themselves and with the yet unknown CP violating phase, $\delta$ appearing in the standard paradigm. The experiment utilizes a wide band beam and provides us with a unique opportunity to utilize different beam tunes at DUNE. We demonstrate that combining information from different beam tunes (low energy, LE and medium energy, ME) available at DUNE impacts the ability to probe some of these parameters and leads to altering the allowed regions in two-dimensional space of parameters considered.

The Deep Underground Neutrino Experiment (DUNE) is a leading experiment in neutrino physics which is presently under construction. DUNE aims to measure the yet unknown parameters in the three flavor oscillation scenario which includes discovery of leptonic CP violation, determination of the mass hierarchy and determination of the octant of θ23. Additionally, the ancillary goals of DUNE include probing the subdominant effects induced by new physics. A widely studied new physics scenario is that of nonstandard neutrino interactions (NSI) in propagation which impacts the oscillations of neutrinos. We consider some of the essential NSI parameters impacting the oscillation signals at DUNE and explore the space of NSI parameters as well as study their correlations among themselves and with the yet unknown CP violating phase, δ appearing in the standard paradigm. The experiment utilizes a wide band beam and provides us with a unique opportunity to utilize different beam tunes at DUNE. We demonstrate that combining information from different beam tunes (low energy and medium energy) available at DUNE impacts the ability to probe some of these parameters and leads to altering the allowed regions in two-dimensional space of parameters considered.

I. INTRODUCTION
In a seminal paper in 1978, Wolfenstein first proposed the possibility that nonstandard neutrino interactions (NSI) could be responsible for conversion of a given neutrino flavour to another even if neutrinos were massless [1]. However, thanks to the wealth of data accumulated by a variety of oscillation experiments covering different energies and baselines, we now have a fairly clear picture that neutrino oscillations occur due to nonzero neutrino masses. The data from most of the oscillation experiments can be nicely explained by invoking three flavors of neutrinos (ν e , ν µ , ν τ ) which are superpositions of the mass states ν 1 , ν 2 , ν 3 with masses m 1 , m 2 , m 3 respectively. The 3 × 3 mixing matrix appearing in the weak charged current interactions is given by, where s ij = sin θ ij , c ij = cos θ ij and δ is the Dirac-type CP phase. The form of U given in Eq. 1 is referred to as the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) parametrization [2]. If neutrinos are Majorana particles, there can be two additional Majorana-type phases in the three flavour case. However, those Majorana phases play no role in neutrino oscillation studies. We have measured the parameters entering the neutrino oscillation framework to a fairly good precision (see the global fit analyses [3][4][5][6][7]). The best-fit values and 3σ range of neutrino mass and mixings deciphered from oscillation data are given in Table I. Yet, there are some unknowns in the standard mass-induced oscillation framework. These include the question of neutrino mass hierarchy (sign of ∆m 2 31 ), the value of the CP violating phase (δ) and determining the correct octant of θ 23 .
On the theoretical side, neutrino oscillations require non-zero masses while neutrinos are massless in the SM. This implies that one needs to go beyond the SM in order to explain the results of oscillation experiments. The minimal way is to have a new physics model which can give rise to nonzero neutrino masses but the interactions are still described by SM. Once we invoke new physics to accommodate neutrino masses, it is only natural to consider the possibility that the neutrino interactions are described by NSI (as was proposed by Wolfenstein [1]). Clearly, a dominant contribution from such interactions is ruled out by the present data [4][5][6][7]. However a subdominant contribution cannot be ruled out given the present accuracy of the neutrino oscillation experiments. Therefore, the idea proposed by Wolfenstein does not hold true in totality in the current times yet his insight remains in the form of subdominant effects due to NSI on neutrino oscillations.
The fact that parameter degeneracies crop up in the presence of standard interactions (SI) has been well recognized since the past two decades or so [8][9][10][11][12][13]. Identification and resolution of parameter degeneracies is crucial for a clean determination of the oscillation parameters. Besides, any new physics sector (such as NSI considered in the present work) introduces a multitude of parameter degeneracies apart from those in the standard case and the structure of parameter degeneracies is far more complex. There has been a vast body of work done on NSI and neutrino oscillations. For a comprehensive recent review on the topic of NSI in the context of neutrino oscillations, we refer the reader to [14]. The idea of subdominant NSI in neutrino propagation affecting the CP violation studies, neutrino mass hierarchy and octant of θ 23 at upcoming long baseline neutrino experiments has received tremendous attention in neutrino physics in the recent years  mainly because our ability to search for subdominant effects has increased substantially due to the precisely designed experiments.
In order to set the stage for the present work, we summarize the most relevant references dealing with the issue of propagation NSI at long baselines and constraining NSI parameters. By carrying out detailed simulations of the DUNE experiment in the presence of new physics, the authors of [16] focus on whether DUNE would be able to distinguish between different kinds of new physics such as propagation NSI and sterile neutrino. In [17], it was shown that DUNE will improve the constraints over some of the propagation NSI parameters by carrying out sensitivity studies and suggested that a combination of DUNE and T2HK would help in resolving degeneracies among standard and NSI parameters. Reference [20] focuses on LBNO and addresses prospects of probing strength of propagation NSI parameters at long baseline experiments as a function of the oscillation channel, baseline length and detector mass. Correlations between propagation NSI and source as well as detector NSI have been studied in [29]. Reference [35] deals with yet another long baseline experiment, T2HKK and discusses how different configurations of T2HKK would be helpful in constraining the propagation NSI. Reference [37] discusses the issue of parameter degeneracies in the presence of propagation NSI and the authors perform a comparison of the potential of DUNE, T2HK and T2HKK in probing some of the NSI parameters. In [41], the author considers a combination of information from atmospheric neutrinos and long baseline experiment T2HK and its impact on constraining the NSI parameters. It should be noted that the studies carried out so far on constraining NSI terms on DUNE has invariably utilized the standard low energy (LE) flux that peaks around the first oscillation maximum for P µe i.e., around 2 − 3 GeV. We advance in this direction by incorporating different beam tunes at DUNE and understand the role of beam tunes in constraining the NSI parameters. In a recent work, high energy beams have been shown to be helpful in distinguishing the NSI scenario from the standard three neutrino scenario [62]. While the new physics context of the present study is that of propagation NSI, our approach is valid for a variety of new physics models.  [4,5] .
The article is organised as follows. In Sec. II, we give the theoretical introduction to neutral current (NC) NSI which is the new physics scenario considered in the present article. We also mention the present constraints on the NSI terms. In Sec. III, we describe the numerical simulation procedure as well as introduce the beam tunes used. In subsection IV A, we first discuss the impact of individual NSI terms on the behaviour of probabilities (P µe and P µµ ) as functions of δ. In subsections IV B and IV C, we analyze the behaviour of the probability difference between NSI and SI as a function of energy as well as δ. In Sec. V, we do a comparative ∆χ 2 analysis to discuss in detail how the higher energy beams in conjunction with the standard low energy beam impact the sensitivities of parameters. Finally, we summarize our conclusions in Sec. VI. In Appendices A and B, we have given the relevant probability expressions that aid in understanding our results. Appendix C contains the SI-NSI event difference plot for some representative choice of parameters.

II. MODEL : NONSTANDARD INTERACTION DURING PROPAGATION
The new physics scenario considered in the present work is that of propagation NSI which impacts the propagation of neutrinos. Such a scenario can be described by a dimension-six operator involving four fermions, where α, β = e, µ, τ indicate the neutrino flavor, f denotes the matter fermions, e, u, d. The new NC interaction terms can impact the neutrino oscillation physics via flavour changing interaction or flavour preserving interaction. From a phenomenological point of view, only the sum (incoherent) of all the individual contributions (from different scatterers such as e, u or d) contributes to the coherent forward scattering of neutrinos on matter. Normalizing to n e , the effective NSI parameter for neutral Earth matter 1 is given by where n f is the density of fermion f in medium crossed by the neutrino and n refers to neutrons. Also, ε f αβ = ε f L αβ +ε f R αβ which encodes the fact that NC type NSI matter effects are sensitive to the vector sum of NSI couplings.
In the presence of NSI, the Hamiltonian in the effective Schrödinger -like equation governing neutrino evolution can be expressed as where ∆m 2 ij are the mass-squared differences. Here a(x) = 2 √ 2EG F n e (x) is the standard charged current (CC) potential due to the coherent forward scattering of neutrinos, n e is the electron number density and ε αβ (≡ |ε αβ | e iϕ αβ ) are complex NSI parameters. U is the PMNS three flavour neutrino mixing matrix (see Eq. 1).
We now mention the constraints on the NC NSI parameters. The combination that enters oscillation physics is given by Eq. 3. Assuming that the errors on individual NSI terms are uncorrelated, model-independent bounds on NC NSI terms ε αβ were given in Ref. [63]. In particular, one obtains the following: for neutral Earth matter. Direct experimental constrains from neutrino experiments on NSI parameters are more restrictive. The SK NSI search in atmospheric neutrinos crossing the Earth found no evidence in favour of NSI and the study led to upper bounds on NSI parameters [64] given by |ε µτ | < 0.033, |ε τ τ − ε µµ | < 0.147 (at 90% CL) in a two flavour hybrid model [65]. The off-diagonal NSI parameter ε µτ is constrained −0.20 < ε µτ < 0.07 (at 90% CL) from MINOS data in the framework of two flavour neutrino oscillations [66,67].
In what follows, we shall adopt a numerical approach to discuss the impact of various NSI parameters. For the sake of simplicity and clarity, we consider one NSI parameter at a time. Wherever analytic description is feasible, we give approximate analytic expressions which are valid in the present context and additional plots which help in understanding the results obtained numerically (for more details, see Appendices A and B).

III. SIMULATION PROCEDURE AND BEAM TUNES
The proposed Deep Underground Neutrino Experiment (DUNE) has a baseline of 1300 km and a 40 kt liquid argon far detector (FD) is placed at an on-axis location. The primary scientific goals of DUNE include the measurement of leptonic CP violation, the determination of the neutrino mass ordering and the precision measurement of the neutrino mixing parameters [59,60,68,69].
In order to simulate DUNE, we use the GLoBES package [70,71] with the most recent DUNE configuration file provided by the collaboration [72] and implement the density profile given by Preliminary Reference Earth Model (PREM) [73]. We assume a total runtime of 7 years with 3.5 years in the neutrino mode and another 3.5 years in the antineutrino mode.  [72,74]. The target is a thin Be cylinder 2 interaction lengths long. The target location is given with respect to the upstream face of Horn 1. The LBNF neutrino beamline decay pipe length has been chosen to be 194 m. Decay pipe lengths of up to 250 m could be accommodated on the Fermilab site and were an option in previous designs of the beamline.
We consider two beam tunes obtained from a G4LBNF simulation [75,76] of the LBNF beam line using NuMI-style focusing.
LE beam: The standard ν µ beam which peaks around a relatively lower energy of ∼ 2.5 GeV (corresponding to the first oscillation maximum for the ν µ → ν e appearance channel) is referred to as an LE beam in our analyses. It is generated by an 80 GeV proton beam delivered at 1.07 MW with protons on target (pot) of 1.47 × 10 21 .

ME beam:
The second beam is has the characteristic that it is larger at higher energies ( 4 GeV onwards) and we refer to this beam as medium energy (ME) beam. The ME beam is generated by a 120 GeV proton beam delivered at 1.2 MW with a pot of 1.1 × 10 21 . beam as used in [72]. ME beam refers to the flux peaking at a higher energy. See Table II for more details.

ME LE
Both the LE and ME fluxes are shown in Fig. 1. The LE flux peaks around 1.5 GeV to 3.5 GeV but after that it falls off rapidly. In contrast, the ME flux is almost flat from 2 -6 GeV and after that it falls off but at a much slower rate compared to the LE flux and it remains substantially higher than the LE flux even beyond 6 GeV. At ∼ 2.5 GeV, the ME flux is ∼ 25 − 35% smaller than the LE flux. Hence, in our analyses of probing the NSI parameters, we use a combination of LE and ME flux together, so as to extract information on new physics from both the lower energy (1 − 3 GeV) and the higher energy ( 4 GeV) regime as much as possible. We compare the results with those obtained using the LE beam only for the same total runtime of the experiment. The beamline parameters assumed for the different design fluxes used in our sensitivity calculations are given in Table II.
Our analysis includes both appearance (ν µ → ν e ) and disappearance (ν µ → ν µ ) channels, simulating both signal as well as background. The simulated background includes contamination of antineutrinos (neutrinos) in the neutrino (antineutrino) mode, and also misinterpretation of flavors, as discussed in detail in [72]. To analyze the NSI scenario, we utilise the GLoBES extension called snu.c which is described in [77,78].
To calculate the sensitivity with which the NSI parameters can be probed, one can define the (statistical) ∆χ 2 as follows 2 : Here, the SI case is treated as true while the NSI parameters are allowed to vary in the test dataset. The sum over the number of channels runs over the ν µ → ν e and ν µ → ν µ channels and the corresponding antineutrino channels, ν µ →ν e andν µ →ν µ . The index j indicates the sum over all the energy bins ranging from E = 0 − 20 GeV. We have a total of 71 bins of non-uniform widths (64 bins with uniform bin width of 125 MeV in the energy range E = 0 − 8 GeV and 7 bins with variable width beyond 8 GeV) [72]. The detector configuration, efficiencies, resolutions and systematic uncertainties for DUNE are listed in Table. III.
Detector details Normalisation error Energy calibration error Signal Background Signal Background DUNE Runtime (yr) = 3.5 ν + 3.5ν νe : 5% νe : 10% νe : 2% νe : 10% 40 kton, LArTPC νµ : 5% νµ : 10% νµ : 5% νµ : 10% We have used the standard oscillation parameters in Table I, taken from Ref. [4,5]. For the neutrino mass hierarchy, we assume a spectrum corresponding to normal hierarchy in the true dataset. Since DUNE has no sensitivity to the solar parameters and since θ 13 is rather well measured by current reactor and long baseline experiments, we keep these values fixed to their current best-fit values, while marginalizing over θ 23 (in the present 3σ range) and δ ([−π, π]), if not plotting them. In addition, we marginalize over the atmospheric mass-squared splitting, ∆m 2 31 , allowing for the two possible mass hierarchies. When studying a non-diagonal NSI parameter, ε αβ , we also marginalize over its corresponding phase, ϕ αβ in the range [−π, π]. Therefore, if we study two non-diagonal complex parameters simultaneously, we marginalize over a total of five parameters.
In our analysis, we consider two diagonal NSI parameters and three off-diagonal NSI parameters with both their moduli and phases. If we also include the yet unknown CP phase, δ, we have a total of nine parameters. We depict ∆χ 2 correlations among these nine parameters (δ, ε ee , |ε eµ |, ϕ eµ , |ε eτ |, ϕ eτ , |ε µτ |, ϕ µτ , ε τ τ ) considering them pairwise at a time and the number of such combinations is 36.

IV. A SCAN OF PARAMETER SPACE AT THE LEVEL OF PROBABILITY
In order to obtain insight into the correlations and degeneracies among the various NSI and SI parameters that may impact the signals at DUNE, the first step is, naturally, to look at the relevant oscillation probabilities. We consider the following oscillation channels that are accessible 3 at DUNE : In what follows, we consider the relevant parameters that include two of the diagonal NSI parameters (ε ee , ε τ τ ) and the moduli and phases of the three off-diagonal NSI parameters (ε eµ ,ε eτ ,ε µτ ). A detailed assessment of the role of individual NSI terms on the different oscillation channels has been carried out in [79,80]. Based on the analyses in [79,80], we can conclude that among all NSI parameters, ε eµ and ε eτ mainly impact the appearance channel (ν µ → ν e ) while ε ee has a milder impact. It is clear that ε eµ enters ν µ → ν e channel. The almost maximal mixing in the 2 − 3 sector ensures that ε eτ also impacts this channel with similar strength as ε eµ (see Appendix A and the discussion in Sec. IV of [79]). Similarly, the disappearance channel (ν µ → ν µ ) is more sensitive to the presence of NSI parameter ε µτ (see Appendix B and the discussion in Sec. IV of [79]).
In the following sub-sections, we perform a scan of the parameter space at the probability level. We first discuss the fixed energy and fixed baseline snapshots of probabilities (subection IV A). We then discuss SI-NSI degeneracies in the context of DUNE as a function of energy keeping δ fixed at the best-fit value (subsection IV B). Further, we go on to the discussion of SI-NSI degeneracies as a function of δ (keeping the energy fixed) in subection IV C.
A. Snapshots of Pµe and Pµµ at fixed energy and fixed baseline In Fig. 2, we fix the baseline at 1300 km and show the impact of NSI parameters 4 on snapshots of P µe and P µµ as a function of δ at certain (appropriately chosen) fixed energy values. This aids in identification of parameters that may have the largest impact at the level of probabilities, though at specific energy values. For the ν µ → ν e channel (top row in Fig. 2), we choose the fixed value of energy to be E = 2.5 GeV. This value corresponds to the first oscillation maximum for P µe . On the other hand, P µµ is very close to zero at 2.5 GeV while it is substantial at higher values of SI |ɛ eμ |, φ eμ = 0.1, 0 |ɛ eμ |, φ eμ = 0.1, π/4 δ, φ eμ ∈ [-π:π] Pμ e (2.5 GeV) 0 0.05 energy. Hence we depict curves for P µµ (bottom row in Fig. 2) at 5 GeV. The grey bands show the variation of the probability when the relevant phases (δ, ϕ αβ ) are allowed to vary in the range [−π, π]. As a reference, the SI case is shown as a solid black line in all the plots. As far as P µe at 2.5 GeV (top row) is concerned, we note that the effect of ε eµ or ε eτ is more pronounced when compared to the other NSI terms. The presence of ε eµ or ε eτ modifies the overall amplitude and the location of the peaks/dips of the probabilities while the presence of a nonzero ϕ eµ or ϕ eτ brings in additional phase shifts. We note that ε µτ has a much smaller effect on P µe . ε ee and ε τ τ also have a miniscule effect on the amplitude of P µe . On the other hand, P µµ at 5 GeV (bottom row) gets affected most by the presence of the ε µτ term. ε eτ , ε ee , ε τ τ have practically no impact on P µµ . ε eµ induces some phase dependence on P µµ .
In what follows, we generalise the above discussion and study the energy dependence of the SI-NSI degeneracies for P µe and P µµ and also vary the NSI terms instead of keeping their values fixed.

B. Energy dependence of the SI-NSI degeneracies
To quantify the impact of NSI terms, let us define a quantity, |∆P αβ | = |P N SI αβ −P SI αβ |(α, β = e, µ), which is absolute value of probability difference between the SI and NSI scenarios.
Our results are given in Fig. 3 in the form of heatmaps as functions of energy and the strength of the NSI parameter for |∆P µe | (top row) and |∆P µµ | (bottom row). The NSI phases are taken to be zero and the standard oscillation parameters have been pinned to their best-fit values (see Table I). If we carefully examine the top row of Fig. 3, we note that |∆P µe | is mostly affected by |ε eµ | and |ε eτ |. Note that, the impact of |ε eµ | or |ε eτ | is most prominent around 2 − 3 GeV. One can derive a useful conclusion here regarding difference in impact of |ε eµ | and |ε eτ | on |∆P µe |. As we go beyond ∼ 4 GeV, |ε eτ | gradually makes |∆P µe | smaller (red region), while |ε eµ | makes |∆P µe | stay at a high value (blue region) which is almost independent of energy. This, in turn, suggests that one may be able to probe ε eµ more effectively than ε eτ by use of higher energy beam tune. The other NSI terms ε µτ , ε ee or ε τ τ do not induce much change, keeping |∆P µe | 0.005 for most of the energy range.
From the bottom row of Fig. 3 corresponding to |∆P µµ |, we note that ε µτ plays an important role. |∆P µµ | is large (blue) in most of the energy range as long as |ε µτ | 0.02. This is to be contrasted with other NSI terms, as even a small value of ε µτ can induce a large impact on |∆P µµ |. As can be noted, a higher energy beam may be able to probe ε µτ via this channel effectively. If we look at the impact of |ε eµ | and |ε eτ |, we note that |ε eµ | gradually makes |∆P µµ | larger at E 5 GeV (indicated by blue region on the top right side of the panel) while |ε eτ | does not seem to impact  |∆P µµ |. Thus, for the disappearance channel as well, it appears that the higher energy beam may prove more useful in probing ε eµ than ε eτ . For an analytic understanding of the energy dependence of ∆P αβ in the presence of NSI, see Appendices A and B. We next consider the case of nonzero phases ϕ αβ . In Fig. 4, heatmaps corresponding to |∆P µe | (top row) and |∆P µµ | (bottom row) in the two-dimensional plane of individual NSI phase (ϕ αβ ) and energy are shown. The moduli of NSI terms (|ε eµ |, |ε eτ | or |ε µτ |) were kept fixed at 0.05. From Fig. 4, we note that |∆P µe | (top row) is most affected by ϕ eµ or ϕ eτ while ϕ µτ has almost no effect. Around 2 − 4 GeV, ϕ eµ and ϕ eτ produce similar qualitative features indicating SI-NSI degeneracy (red band) occurring at a pair of values given by ϕ eµ ≈ 0, ±π and ϕ eτ ≈ −0.6π, 0.4π. At energies beyond 4 GeV, |∆P µe | is very small (∼ 0) almost uniformly in the presence of ϕ eτ while in the presence of ϕ eµ , it still exhibits a relatively high value (0.005 − 0.01) in large region of the parameter space. For a quantitative understanding of this feature, we refer the reader to Appendix A. For |∆P µµ |, Fig. 4 (bottom row) shows that it is affected most significantly by ϕ µτ , showing sharp SI-NSI degeneracy around ϕ µτ ≈ ±π/2. This arises because of the fact that |∆P µµ | ∝ cos ϕ µτ (See Eq. B3 in Appendix B). We also note that |∆P µµ | remains close to zero in the presence of ϕ eτ and shows moderate variation for in the presence of ϕ eµ .
Finally, we would like to mention that the qualitative features of Fig. 4 remain unchanged even if the moduli of the relevant off-diagonal NSI terms (|ε eµ |, |ε eτ | or |ε µτ |) are increased.

C. δ-dependence of SI-NSI degeneracies
In Figs. 5 and 6 we depict the heatmaps for |∆P µe | and |∆P µµ | in the δ −ε αβ plane. In these plots, the first (second) row corresponds to a fixed energy of 2.5 (5) GeV. We can derive the following conclusions in connection with P µe (see also Appendix A) :  Table I.
• In the case of ν µ → ν e channel (Fig. 5), the NSI terms ε eµ and ε eτ have relatively larger impact than the other NSI parameters. For ε eµ and ε eτ , the degenerate regions (|∆P µe | 0.05) are narrowly concentrated around a pair of values of δ (see Table IV below) These sharp SI-NSI degenerate regions exist even for 5 GeV but at somewhat different values of δ. For ε eτ , the degenerate region seems to be larger at 5 GeV in contrast to 2.5  GeV. This is not seen in the case of ε eµ (this observation is consistent with Fig. 3). Note that the locations of SI-NSI degenerate regions is roughly independent of the size of |ε eµ | and |ε eτ |. For |ε µτ |, the degenerate region is broader and shows a soft feature of peaking at δ ≈ ±π for 2.5 GeV. For ε ee and ε τ τ , the degenerate regions have similar structure showing no CP phase dependence.
• For the ν µ → ν µ channel (Fig. 6), as mentioned earlier, it is more appropriate to look at 5 GeV (the bottom row). As expected, ε µτ has the largest impact and its effect is independent of the CP phase, δ (see also Appendix B). The impact of ε eµ is also important with two sharp peaks occurring around δ ≈ ±π/2. The other terms such as ε eτ , ε ee , ε τ τ have almost no effect at 5 GeV (here also the results are consistent with Fig. 3).
To complete the discussion, we now discuss the effect of non-zero phases. We keep the moduli of the respective NSI terms fixed at |ε αβ | = 0.05 and plot heatmaps corresponding to |∆P µe | and |∆P µµ | in the ϕ αβ − δ plane in Fig. 7 and Fig. 8 respectively. As before, we show our results for two different values energy, 2.5 GeV (top row) and 5 GeV (bottom row). We make the following observations from these plots: • In the case of the ν µ → ν e channel (Fig. 7), we see degenerate regions in the case of ϕ eµ and ϕ eτ (where |∆P µe | 0.005) slanted at an angle of 135 o . In the case of ϕ µτ , |∆P µe | remains close to zero and stays within ∼ < 0.005 in the entire ϕ µτ − δ space 5 . For 5 GeV, the pattern remains very similar for ϕ eµ , but the extent of degeneracy increases for ϕ eτ , as expected from our previous analyses.
• In case of the ν µ → ν µ channel (Fig. 8), we focus on the bottom row. Here ϕ eµ shows the SI-NSI degenerate regions roughly mimicking straight lines at 135 o slope, whereas ϕ eτ shows no effect. ϕ µτ manifests itself by rendering |∆P µµ | to a much larger value ( 0.02) for most of the parameter space, but there exist two sharp degenerate regions occurring at ϕ µτ ≈ ±π/2 with no δ dependence.

V. PROBING THE NSI PARAMETER SPACE AT THE LEVEL OF χ 2
In the present section, we numerically explore the NSI parameter space at the level of χ 2 using the standard LE as well as ME beam tunes. Our main results are summarized in Fig. 9 where we depict contours at a confidence level (c.f.) of 99%. The solid cyan (black hatched) contours correspond to LE (LE + ME) beams. More specifically, the regions enclosed by these contours depict the regions where there is SI-NSI degeneracy for those pair of parameters. Below, we discuss some noteworthy features as can be observed from Fig. 9  beam tune is used (cyan region) and when a combination of low and medium energy (LE + ME) beam tune is used (black hatched region), keeping the total runtime same (3.5 years of ν + 3.5 years ofν run) for both scenarios. In the latter case, the total runtime is distributed between the LE beam (2 years of ν + 2 years ofν) and the medium energy beam (1.5 years of ν + 1.5 years ofν). The panels with a light yellow (white) background indicate significant improvement (no improvement) by using LE + ME beam over using LE only. The numbers in the light yellow shaded panels correspond to the area lying outside the contour for the two cases (cyan for LE and black for LE+ME) expressed as a percentage of the total parameter space plotted. These numbers quantify the improvement over the LE only option when the ME beam tune is used in conjunction with the LE beam tune in these panels.
1. Let us first consider the panels with ε eµ (either |ε eµ | or ϕ eµ or both) which are shown in light yellow colour. We note that use of different beam tunes (ME in conjunction with the LE beam) offers visible improvement of results (shrinking of contours) in these pairs of parameters. This is one of the key results of the present article. In order to explain the observed pattern, let us recollect from Figs. 3 and 4 that the presence of |ε eµ | or ϕ eµ leads to large difference between SI and NSI scenarios even at larger values of energies i.e., E 4 GeV. Thus, with the LE+ME option we are able to place tighter constraints on the parameter space corresponding to parameters |ε eµ | and ϕ eµ .
From Eq. 7, the ∆χ 2 in the presence of two NSI parameters, say, a and b, can be written as 7 : For e.g., if we focus on the |ε eµ |-|ε eτ | plane, we have where the sum is over all the energy bins (0−20 GeV) and the minimization is performed over δ, θ 23 , ∆m 2 31 , ϕ eµ , ϕ eτ . From the probability level discussion (Fig. 3), we can assess the impact of the NSI terms |ε eµ | and |ε eτ | on P µe and P µµ . In the case of |∆P µe |, at low values of energy, the impact of the two NSI parameters is quite similar. But, at higher energies, the effects due to |ε eµ | tend to be larger than effects due to |ε eτ |. This means that ME beam is expected to alter the degenerate region more in the case of |ε eµ | and less in the case of |ε eτ |. That the smaller contribution from the disappearance channel is in the same direction as the larger contribution from the appearance channel (with |ε eµ | and |ε eτ | acting in opposite directions) can also be seen from the plot.
2. We next consider the remaining panels in which we see that there is very little or no improvement of results after using the ME beam along with the LE beam. If we look at the pair of parameters, |ε eτ | − ε ee , ϕ eτ − ε ee , ε τ τ − ε ee , ϕ eτ − |ε eτ | and ε τ τ − |ε eτ | in particular, we note that the degenerate regions get enlarged slightly. This is because of the fact that the the presence of ε eτ , unlike ε eµ , actually adds to the SI-NSI degeneracy at higher energies (see Figs. 3 and 4 and the discussions in Sec. IV B).
3. For the panels with |ε µτ | and ϕ µτ as one of the parameters, there is very marginal improvement (except when |ε eµ | or ϕ eµ is present) in the degenerate contours using the LE+ME beam. To see how the ∆χ 2 arises in panels showing the parameter space associated with |ε µτ |, let us take for example, the pair of parameters, |ε µτ | and |ε eτ | and express the ∆χ 2 (Eq. 9) as where the sum is over all the energy bins (0−20 GeV) and the minimization is carried over δ, θ 23 , ∆m 2 31 , ϕ µτ , ϕ eτ . Now, from Eq. B3, we know that in leading order, |∆P µµ (ε µτ )| is independent of δ and is directly proportional to cos ϕ µτ . Minimization over ϕ µτ ∈ [−π, π] will always then find the constant, energy-independent value of ϕ µτ ≈ ±π/2 which makes the ∆χ 2 contribution due to P µµ vanishingly small 8 . Thus, even when |ε µτ | is present, the ∆χ 2 receives a dominant contribution from the ν µ → ν e channel. This is more clear from the panels showing the parameter space associated to ϕ µτ (i.e., where ϕ µτ is not marginalised). The magnitude of ∆χ 2 in such panels is dominantly contributed by the ν µ → ν µ channel for all values of ϕ µτ ≈ ±π/2. But around ϕ µτ ≈ ±π/2, the contribution from the ν µ → ν µ becomes very small and the ν µ → ν e channel dominates, as we have also verified numerically. This explains the appearance of degenerate contours at ϕ µτ ≈ ±π/2 as well.
4. All the parameter spaces showing ε ee (entire 2nd column and the top panel of the 1st column) have an additional degeneracy around ε ee ≈ −2, in addition to the true solution at ε ee ≈ 0. This extra solution comes due to the marginalisation over the opposite mass hierarchy. Similar degeneracy has also been observed in previous studies: in [25,36,82] (in the context of NSI) and also in [83] in the context of Lorentz violating parameters.

VI. CONCLUDING REMARKS
In [62], the authors combined the experimentally feasible option of using a combination of beam tunes and demonstrated that it was possible to extricate any two physics scenarios more efficiently using experimental handles. In the present article, we address the question of constraining the parameter space of NSI parameters at DUNE by exploiting wide band nature of the beam. We systematically study correlations among various parameters using two beam tunes (LE and ME) and illustrated that to probe a subset of NSI parameter space more effectively, it is advantageous to use a combination of LE and ME tuned beams as opposed to using only the standard LE beam tune.
We provide a systematic and comprehensive description of the overall impact of the NSI parameters on the relevant probabilities (for ν µ → ν e and for ν µ → ν µ ) as a function of energy as well as the CP phase. In the Appendices, we provide analytic expressions of all the relevant expressions for the SI-NSI probability differences in the presence of individual NSI parameters (taken one at a time). These aid in our understanding of the dependencies of oscillation probabilities. We then quantify the differences in terms of a ∆χ 2 quantity and connected the features obtained to the probability level description. In Fig. 9, we have illustrated the ∆χ 2 correlations among the various parameters in the new parameter space appearing in the presence of NSI at a confidence level of 99%. Our key findings can be summarized as follows. The degenerate contours in the space associated with parameters, |ε eµ | and ϕ eµ (shown as panels shaded in light yellow colour in Fig. 9) shrink significantly when we use the LE+ME beam as opposed to LE beam alone. For a quantitative estimate of the improvement, one can compute the area of the parameter space outside each contour (i.e., above the confidence level of 99%) and express the area as the percentage of the total parameter space plotted. It is evident from the pair of numbers (cyan for LE and black for LE+ME) indicated in the light yellow panels that the LE+ME beam leads to improvement over the LE beam. For the remaining NSI parameters, we see marginal or no improvement in terms of constraining the parameters using LE+ME beam in comparison with LE beam. Our detailed analysis also provides explanation for distinguishing features of the ∆χ 2 contours for different parameters.
A few remarks on connection with the existing work that deal with constraining NSI parameters at DUNE [16,17] are in order. It should be noted that the standard beam (available in 2015) was used for these analyses. It can be observed that the contours in our analyses indicate better resolution capability of DUNE and they roughly resemble the contours of [16,17] in shape. The slight differences in the contours may arise from the difference in the experimental inputs such as beam configuration, detector details, exposure and the best-fit values used.
Although the entire analysis has been carried out in the context of NSI and DUNE, we would like to mention that the approach is fairly general and can be easily translated to other new physics contexts such as non-unitarity, CPT violation, Lorentz violation etc. We point out that if we utilize the full potential of a given long baseline experiment such as DUNE which has wide band beam and allows for tunability of beam, we can reduce the degeneracy (using the same experiment with multiple beam tune options) among some relevant choice of parameters in the parameter space as is suggested by the probability level discussion in Sec. IV.
Appendix A: Analytic understanding of the behaviour of ∆Pµe Here we look at the expressions for probability difference between SI and NSI and make an attempt in understanding how the individual NSI parameters affect the SI-NSI degeneracy. We calculate these expressions by making use of the probability expressions from [79] upto first order in ε αβ 's. Using the expressions for P µe in presence of a single NSI parameter (|ε eµ |, |ε eτ | or ε ee ) we arrive at the following three equations: ≈ 2A∆ sin ∆|ε eτ |s 2(23) 2s 13 s 23 D eτ 1 sin(δ + ϕ eτ + γ eτ 1 ) − αs 2(12) c 23 D eτ 2 sin(γ eτ Here A = a/∆m 2 31 = 2 √ 2EG F n e /∆m 2 31 . By making the substitution A → A(1 + ε ee ) [18] we also have, where ∆ = ∆m 2 31 L 4E . When ∆P µe becomes close to zero, it becomes difficult to separate NSI from SI and we have a SI-NSI degeneracy. We plot the terms in Eq. A1, Eq. A2 and Eq. A3 as functions of δ for an energy of 2.5 GeV and also at a higher energy of 5 GeV in Fig. 10 with fixed values of the NSI amplitude and zero NSI phase as indicated in the figure. For ε eµ or ε eτ , the second term (blue) is very small (scaled down by the additional factor α ∼ 10 −2 compared to the first term) and also independent of the CP phase δ. It is the first term (green) which mainly dictates the behaviour of ∆P µe in presence of ε eµ or ε eτ . We note the locations (see Table V) where the overall value of ∆P µe (red) becomes zero in Fig. 10. These locations are indeed consistent with the locations of the red spikes in Fig. 5 as listed in Table  IV. The origin of these special values of δ can easily be understood as follows.
On a closer inspection of the first term in Eq. A1 and Eq. A2, we observe that it is proportional to D eµ 1 for ε eµ and to D eτ 1 for ε eτ . From Fig. 11 (left panel), we observe that around 2.5 GeV both D eµ 1 and D eτ 1 have similar magnitude 9 . But as the energy increases further D eµ 1 keeps on increasing while D eτ 1 decreases. This indicates that at higher energies, |∆P µe (|ε eµ |)| increases while |∆P µe (|ε eτ |)| becomes smaller. This explains why the degeneracy increases for higher energy in presence of ε eτ compared to ε eµ in the ν e appearance channel, as observed from our simulation earlier(see Fig. 3).
In Fig. 12, we plot the terms of Eq. B1, B2 and B3 for two fixed energies 2.5 GeV and 5 GeV as functions of δ. We have already observed before that for the disappearance channel, it is the higher energy range that contributes more. To understand ∆P µµ , we will thus refer to the more relevant bottom row of Fig. 12. It is clear from the figure (first and second column) that the two terms for ∆P µµ act in the same direction for ε eµ (thereby increasing the overall |∆P µµ |), but show opposite behaviour for ε eτ , leading to an overall very small ∆P µµ through cancellation in the latter case. Looking back at Eqns. B1 and B2, we note that both term I and II are roughly proportional to cos δ. But due to the presence of a relative sign in the coefficient of A|ε eτ | in the second term, this behaves in almost opposite direction of the first. 10 This has an interesting consequence that ∆P µµ (ε eτ ) is significantly small at higher energies unlike ∆P µµ (ε eµ ). This is also manifestly evident from our simulation (Fig. 3: bottom row, first and second columns). Additionally, in Fig. 6 (bottom row, first and second column) we have also observed the appearance of two red peaks around ±π/2 for ε eµ and mostly reddish region (implying very small ∆P µµ ) in presence of ε eτ . 10 s 13 ∼ 0.15, As 23 |εeµ| ∼ Ac 23 |εeτ | ∼ 0.03, D(|εeµ|) ∼ D(|εeτ |) 0.15 αs 13 s 2(12) ∼ 0.004, αs 2(12) ∼ 0.027, θ(|εeµ|) ∼ θ(|εeτ |) 10 o Thus in the first term of the Eq. B1 and Eq. B2, cos δ part is dominating and in the second term, θ(ε αβ ) is very small,-making the overall ∆Pµµ approximately proportional to cos δ for ease of understanding. Finally, we see from Eq. B3 and the corresponding third column of Fig. 12 that ∆P µµ (ε µτ ) is independent of the CP phase δ and its value is quite significant (except around 2.5 GeV) compared to that in presence of ε eµ or ε eτ . This corroborates the observations in Fig. 3 (bottom row, third column) and Fig. 6 (third column). In order to illustrate the SI-NSI degeneracy at the level of event rates, we can define the following quantity where N αβ stands for the number of events for ν α → ν β . The results are shown in Fig. 13. The top row depicts the