New physics from oscillations at the DUNE near detector, and the role of systematic uncertainties

We study the capabilities of the DUNE near detector to probe deviations from unitarity of the leptonic mixing matrix, the 3+1 sterile formalism and Non-Standard Interactions affecting neutrino production and detection. We clarify the relation and possible mappings among the three formalisms at short-baseline experiments, and we add to current analyses in the literature the study of the ${\nu}_{\mu} \rightarrow {\nu}_{\tau}$ appearance channel. We study in detail the impact of spectral uncertainties on the sensitivity to new physics using the DUNE near detector, which has been widely overlooked in the literature. Our analysis shows that this plays an important role on the results and, in particular, that it can lead to a strong reduction in the sensitivity to sterile neutrinos from ${\nu}_{\mu} \rightarrow {\nu}_e$ transitions, by more than two orders of magnitude. This stresses the importance of a joint experimental and theoretical effort to improve our understanding of neutrino nucleus cross sections, as well as hadron production uncertainties and beam focusing effects. Nevertheless, even with our conservative and more realistic implementation of systematic uncertainties, we find that an improvement over current bounds in the new physics frameworks considered is generally expected if spectral uncertainties are below the $5\%$ level.


Introduction
The neutrino oscillation picture is still far from being complete. While T2K [1] and NOvA [2] will keep collecting data in the coming years, new facilities will be needed to test at high confidence level whether CP is violated in the leptonic sector, and to determine the neutrino mass ordering with a statistical significance above the 5σ level. For that purpose, before the end of the decade two major long-baseline neutrino oscillation experiments will go online: the Deep Underground Neutrino Experiment [3] (DUNE) in the US, and the Tokai-to-Hyper-Kamiokande [4,5] (T2HK) experiment in Japan.
With their larger detectors and more powerful beams, the upcoming generation of oscillation experiments will also be able to thoroughly test the robustness of the standard three-family oscillation framework at an unprecedented level of precision. In fact, while there is plenty of evidence indicating that at least two of the Standard Model (SM) neutrinos are massive, the exact mechanism responsible for their generation is still unknown and requires the addition of new degrees of freedom to the SM particle content, which may lead to the observation of new physics effects. The most minimal extension that can generate neutrino masses is the inclusion of Majorana right-handed neutrinos through the celebrated type-I Seesaw mechanism [6][7][8][9] which leads to unitarity deviations for the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) 3×3 active-light sub-block of the full neutrino mixing matrix and, depending on the mass scale, to additional oscillation frequencies in the probabilities [10][11][12][13][14]. Alternatively, following a completely model-independent perspective, one could add higher-dimensional effective operators to the SM Lagrangian stemming from a new physics theory at high energies. Such operators may induce new interactions affecting neutrino production, propagation or detection processes, usually referred to as Non-Standard neutrino Interactions (NSI) [15][16][17][18][19][20].
It is clear that, in order to search for new physics, systematic uncertainties should be reduced as much as possible. At long-baseline neutrino oscillation experiments, these are relatively large and stem from two key limitations: (i) our inability to compute the neutrino flux analytically; and (ii) the lack of a microscopic model describing the neutrino-nucleus interaction cross section at GeV energies, where nuclear effects are very relevant. In order to reduce these, long-baseline experiments typically use a near detector (ND), located at a sufficiently short baseline (typically L ND ∼ O(500) m) so that standard neutrino oscillations have not yet developed. The ND data is used to determine the unoscillated spectrum with high accuracy, and its combination with the far detector (FD) data leads to a partial cancellation of systematic uncertainties (see, e.g., Ref. [21] for a review on the theoretical and experimental challenges that have to be faced for this to take place).
A priori, the ND data alone can be used to constrain several of the new physics scenarios outlined above and, in some cases, can provide complementary information to that available at the FD. This is clear in the case of an eV-scale sterile neutrino: for example, at DUNE the new oscillation frequency would roughly match the ND distance and induce an oscillating pattern [22][23][24], while at the FD this effect would be completely averaged-out. In the case of NU of the mixing matrix (or, equivalently, NSI in production and detection), the ND data can be used to probe parameters which would otherwise essentially be inaccessible at the FD due to its much lower statistics. Remarkably, this type of new physics effects can induce a neutrino flavor transition already at zero distance. Although a zero-distance effect in the context of non-unitarity of the PMNS matrix has been studied in the context of Neutrino Factories [25][26][27][28][29], to the best of our knowledge this issue has not been explored for the DUNE ND, with the notable exception of Ref. [24] where different configurations for the ND complex were considered. Our work differs from the study in Ref. [24] in several aspects. In particular, we make use of the latest DUNE configuration details in our simulations (beam configuration, detector efficiencies and resolutions, etc). Additionally, while in Ref. [24] the zero-distance effect was only considered for the ν µ → ν e channel, our analysis also includes a study of the sensitivity to anomalous ν τ appearance in the ND. While some ν τ would in principle be produced at the target from the decay of heavy mesons and τ leptons, the flux at DUNE would be too small to lead to an observable number of events in the ND. Thus, the observation of a ν µ → ν τ signal at the ND would automatically point to new physics.
However, while the ND can be used on its own, from the reasoning above it also becomes clear that it will be affected by a much larger level of systematic uncertainties than FD data. In fact, while some of these uncertainties will be correlated among different energies (leading to an overall "rescaling" effect for the event rates), in some cases they can lead to modifications of the expected shape of the event spectrum. Therefore, when studying the sensitivity of neutrino NDs to new physics it is important to distinguish between normalization and shape uncertainties, as these will affect the sensitivity to experimental observables differently. It is however remarkable that most studies of the sensitivity to new physics using the ND at long-baseline experiments only include normalization uncertainties. In this work, we study in detail the impact of shape uncertainties on the sensitivity of the DUNE ND to three main new physics scenarios: non-unitarity (NU) of the PMNS mixing matrix, sterile neutrino oscillations, and Non-Standard Interactions (NSI) affecting neutrino production and detection processes. In particular, we will quantify the expected loss in sensitivity when shape uncertainties are included in the analysis, which stresses the importance of keeping these under control. In this sense, our approach is similar to the one adopted in Refs. [24,30]. However, note that in Ref. [30] the focus was put on the far detector analysis at DUNE and the sensitivity to NSI in neutrino propagation, while in Ref. [24] the authors considered a different set of near detector configurations (baseline, energy resolution, etc).
The paper is organized as follows. Section 2 presents the formalism used and clarifies the relation and possible mapping between the NSI, NU and sterile neutrino formalisms for short-baseline experiments. The most important simulation details, event rates and our choice of systematic uncertainties is discussed in Sec. 3, while our results are presented in Sec. 4. We summarize and conclude in Sec. 5. Appendices A, B and C include more technical details.

Theoretical framework and notation
As mentioned above, the simplest extension able to account for the light neutrino masses and leptonic mixing observed in neutrino oscillation experiments, consists in the addition of singlet fermions to the SM field content. The light neutrino masses and mixing can be successfully generated via the seesaw mechanism, or its different low scale realizations as the inverse and direct seesaw models [31][32][33][34][35], for masses of the new fermions in the range between the O(eV) and up to (near) the GUT scale. These new states are typically named in the literature as right-handed neutrinos or heavy neutral leptons (HNLs) when their mass lies around or above the electroweak scale, and sterile neutrinos when this new physics scale is close to the O(eV) region as relevant for the LSND [36], MiniBooNE [37] and reactor anomalies [38,39] (see e.g., Refs. [40][41][42][43][44][45] for recent global fits in the context of light sterile neutrinos).
In all generality, the full mixing matrix connecting the mass and flavor basis (including light and heavy states) can be written as Here N is a 3 × 3 non-unitary matrix corresponding to the PMNS active-light sub-block, and Θ is the 3×n sub-block parameterizing the mixing between active and heavy neutrinos, with n the number of new states. The R (S) sub-block gives the mixing between the sterile and light (heavy) states and does not play any role in neutrino oscillations. The phenomenology strongly depends on the mass scale of the new states. According to their impact on neutrino oscillations, we can essentially distinguish two possible scenarios that we will describe below in more detail: (i) NU generated by new states above the electroweak scale, which can be integrated out from the low energy spectrum so that the low energy new physics effects are directly encoded in the active-light PMNS sub-block N (which is no longer unitary); and (ii) kinematically accessible sterile neutrinos, which can be produced in the neutrino beam and lead to additional oscillation frequencies in the probabilities.
Notice that even in the sterile neutrino case, the N sub-block is not unitary and thus this scenario can be considered as a source of deviations from unitarity from low scales. Indeed, in [46] it was shown that when the new squared-mass differences are large enough as compared with the ratio between the experiment's baseline and energy, the effects in longbaseline neutrino oscillations coincide at leading order (linear on the parameters quantifying the deviations from unitarity) with those generated in the NU scenario stemming from heavy scales. In such a case, an averaged-out limit in which the oscillation probabilities become independent of the new frequencies is achieved and the effects in the far detector are virtually equivalent for both cases. As we will show below, for the appearance channels in short-baseline experiments it is necessary to go beyond the linear order, which results into a factor of two difference in the oscillation probabilities between the two scenarios in this case.

Parameterization
The most popular parameterizations of the non-unitary matrix N are given by where η is an hermitian matrix [47,48], while (I − T ) is a lower triangular matrix [49][50][51][52] given by 1 1 Notice that in [51] the subindices for the α parameters are indicated by numbers while we use flavor indices. In addition, here α ββ are small parameters which directly parameterize the deviations from unitarity while in [51] αjj are close to one, in such a way that αjj=1-α ββ .
and U and U are unitary matrices that resemble the standard PMNS matrix up to small corrections driven by α and η. Both parameterizations are equally general and a mapping between them can be found in [46]. In this work we will use the lower triangular parametrization but our results can thus be trivially mapped to the hermitian parameterization. As we will see, our analysis can be performed in terms of the α parameters for both scenarios, NU sourced from heavy new physics and sterile neutrino oscillations. In the literature, the sterile neutrino analyses are typically done as a function of the mixing between the active and heavy (mostly sterile) states. The connection between both formalisms can be easily done just considering the unitarity of the full mixing matrix UU † = I, which implies the following simple relation between the N matrix and the active-heavy mixing: (2.4) where in the right hand side we have introduced the triangular parameterization for N given by Eqs. (2.2) and (2.3). From this equation we arrive to Even though both scenarios can be studied with the same formalism using the same parameterization, the bounds applying to each case are remarkably different. If the deviations from unitarity are generated from very heavy scales, integrating out the heavy states leads to modifications of the charged-current and neutral-current couplings of the active neutrinos and very stringent constraints are derived from precision electroweak and flavor searches [25,[53][54][55][56][57][58][59][60]. On the other hand, when the new states are light enough to be kinematically produced with negligible masses as compared with the energy scale of these experiments, unitarity is effectively restored in those observables and the just mentioned constraints do not apply. In such a case, the present bounds stem from neutrino oscillation experiments and are less stringent. A summary of the current bounds in both scenarios is shown in Tab. 1, extracted from Ref. [46]. Note that the constraints in the middle column apply for new squared-mass differences associated to the sterile neutrinos in the range ∆m 2 100 eV 2 , and are thus relevant for both near and far detectors of most longbaseline experiments when the sterile neutrino oscillations are in the averaged-out regime. The bounds shown in the right column apply for ∆m 2 ∼ 0.1 − 1 eV 2 , being relevant if the sterile neutrino oscillations are in the averaged-out regime only in the far detector.

Non-Unitarity from new physics above the electroweak scale
In this scenario the heavy states are integrated out from the low energy spectrum and are thus not kinematically accessible in the experiment. We will focus our analysis on the sensitivity provided by the DUNE ND complex and, therefore, we will restrict our study to  Table 1.
Upper bounds on the Non-Unitarity (NU) framework using the α parametrization, extracted from Ref. [46]. The constraints are shown at 2σ and 95% CL (1 d.o.f.) for NU stemming from very heavy scales and low scale physics (averaged-out light sterile neutrinos) respectively. The value quoted in parenthesis for the α µe element corresponds to the bound obtained from µ → eγ. The limits for the off-diagonal parameters without a reference are derived indirectly from bounds on the diagonal parameters via |α αβ | ≤ 2 √ α αα α ββ . See [46] for further details.
neutrino oscillations in vacuum 2 . The associated oscillation probability is given by [25,46] (2.6) where L and E correspond to the baseline and the neutrino energy, respectively. Notice that this theoretical probability should be convoluted with the neutrino flux and cross section in order to obtain the number of events. At the experimental level, in a far detector analysis the oscillation probability is obtained from the ratio between the number of events observed in the FD and an extrapolation of the results from the ND. Taking into account the corrections from NU, this results in the following experimentally inferred oscillation probability [46] (2.7) The normalization factor can play a relevant role in long-baseline neutrino oscillation studies [46], which is often overlooked in the literature. However, at very short distances, as the ones considered in this work, the theoretical and experimentally inferred oscillation probabilities for the appearance channels coincide.
In particular, we get 3 : (2.8) From this equation it is clear that, if deviations from unitarity are induced by new physics above the electroweak scale, there is a non-zero probability of observing flavor transitions already at zero distance. This renders the neutrino NDs as a powerful tool to constrain or potentially discover a new physics signal encoded in α γβ .

Sterile neutrinos & Non-Unitarity from new physics at low scales
If the new states are light enough to be kinematically produced in the neutrino beam, they can directly participate in the oscillations. Therefore, a priori, the oscillation phenomenology is expected to be different from the NU case considered above. However, it has already been shown that in the average-out limit both scenarios can be virtually equivalent at the phenomenological level regarding the far detector physics [46]. Here we will focus on the ND physics paying particular attention to the differences or similarities between the two approaches. For simplicity, for the sterile neutrino framework we will consider the 3 + 1 case in which only one extra sterile neutrino is added. At very short baselines (as is the case for the DUNE ND), the three-family oscillation frequencies are very suppressed and the oscillation probabilities can be written as to a very good approximation. It is common in the literature to parametrize these probabilities in terms of effective mixing angles for the appearance and disappearance channels, as In the averaged-out regime (∆m 2 41 L/E 1) the oscillations are too fast to be distinguished at the detector 4 , and we find Notice that a normalization factor as the one present in Eq. (2.7) would not play any role here since it would only give a subleading contribution to Eq. (2.8). 4 See appendix C for details on the numerical implementation of this regime in the simulations.
where in the right-hand side we have just introduced Eq. (2.5) particularized to the n = 1 case under consideration.
In the absence of an auxiliary detector with an even closer location to the neutrino source, the ND has a sensitivity to the diagonal α ββ parameters through disappearance channels which is ultimately determined by the normalization uncertainty. Due to the limited knowledge of the neutrino flux and cross section, the region of α ββ that can be experimentally probed by the ND is already ruled out by present experiments. Therefore, in this regime, only the appearance channels are relevant to constrain the new physics effects with the ND.
Furthermore, we would like to stress that the difference between Eq. (2.8) and Eqs. (2.11)-(2.12) is only a global factor 2. This means that the effects at the ND produced by NU stemming from heavy new physics and averaged-out sterile neutrinos are practically indistinguishable. It is potentially possible to distinguish if NU effects measured in the ND are generated by heavy or light new physics [66] if complementary observables are added to the analysis. However, when the corrections are sourced by new physics above the electroweak scale, the current constraints from other observables become very stringent (see Tab. 1) and essentially exclude the possibility of observing any signal at DUNE. Therefore, we will focus on the case in which the deviations from unitarity are generated by new physics at low scales (i.e., sterile neutrinos in the average-out limit). For brevity, in the rest of this work we will refer to this scenario simply as NU. In any case, our results can be trivially recasted to the heavy NU case rescaling our sensitivity limits by the corresponding factor of two.

NSI
The Non-Standard neutrino Interaction (NSI) framework is a model-independent phenomenological approach in which the new physics effects in neutrino oscillation experiments are parameterized in terms of four-fermion effective operators (see e.g. Refs. [67,68] for a QFT description of NSI interactions in the context of SMEFT and reactor experiments). While obtaining sizable NSI at low energies from new physics at high energies without entering in conflict with charged lepton processes is challenging [69,70], it has been shown in recent literature that a possible way out of this argument is the inclusion of new degrees of freedom at low scales [71][72][73][74] or in radiative models of neutrino masses [75] (see [76] for a recent review on the viability of NSI models in this context). In this work we will follow a purely phenomenological approach and report the expected bounds on NSI without considering the underlying theory that can lead to the four-fermion effective operators at low energies.
NSI would lead to corrections in neutrino production, neutrino detection, or in neutrino propagation through matter. Since we are considering the ND complex, there are no effects in neutrino propagation and only NSI in production and detection will be discussed here. The appearance oscillation probability (γ = β) in presence of NSI can be parameterized as It is clear that, analogously to the NU case, a flavor transition can already occur at zero distance in presence of NSI. The main difference between these two approaches at the phenomenological level is that production and detection are correlated in the NU framework, while this is not generally the case in the NSI scenario. Indeed, NU can be considered a particular case of NSI in production and detection. In particular, comparing the above equation with Eqs. (2.11)-(2.12) it is straightforward to see that there is a mapping for the appearance channel physics between NSI and the low scale NU frameworks [46]: The NSI effects in disappearance channels show up already at linear order in , P ββ =

Simulation details
This section describes the main details used in our simulations for the DUNE experiment.
In order to simulate the expected event rates at the DUNE ND, we use the same configuration as in the DUNE Technical Design Report (TDR), see Refs. [3,77]. All our calculations are performed using the GLoBES library [78,79], together with the files provided by the collaboration as ancillary material with Ref. [80]. Below we discuss the main relevant aspects of the simulation, while we refer the interested reader to Refs. [3,77,80] for additional details.

Beam configuration and exposure
The main experimental configuration details are summarized in Tab. 2. In particular, for the nominal running mode we consider equal exposure in neutrino and antineutrino running modes (3.5 years each, giving a total of 7 years of data taking), with a 1.2 MW beam. The nominal neutrino flux for DUNE is shown in Fig. 1 as a function of the neutrino energy (dashed red line). For comparison, the solid green line in the same plot shows the ν τ charged-current (CC) cross section which, due to the large τ mass, does not kick in until the neutrino energy is at least 3 GeV. As shown in the figure the nominal neutrino flux peaks at much lower energies and, as a result, the ability to probe ν µ → ν τ conversion at the ND is going to be severely limited by statistics. For this reason, we have also considered adding 3.5 additional years of running time (in neutrino mode) using the High-Energy (HE) beam configuration. This flux peaks at considerably higher energies as shown by the dotted blue line in Fig. 1 and therefore turns into a higher event rate for CC ν τ interactions.

Simulation of the ν τ -like event sample
We use the ν τ cross section provided with the ancillary material of Ref. [81], which corresponds to the ν τ CC cross section on Ar simulated with GENIE v2.8.4 [82,83]. Regarding the detection and reconstruction effects, the collaboration does not provide efficiencies or energy smearing matrices for a ν τ CC signal. We describe below our treatment of signal and backgrounds for the ν τ -like event sample. Once a τ lepton has been produced at the detector, it will decay promptly. Given the neutrino energy is in the few GeV range, the distance traveled by the τ is not sufficiently long to allow for a successful particle identification at DUNE. However, its decay will still leave a visible signature in the detector: these events are characterized by either a hadronic shower or a charged lepton, plus missing energy carried away by the outgoing neutrino (or neutrinos, in the case of a leptonic decay). Given that the τ branching ratio into the hadronic decay channel is much larger than that of the leptonic channels (65%, compared to a 17% for each of the leptonic decay channels [84]), in our analysis we include only those events where the τ decays hadronically. 5 For these events the main background comes from neutrino NC interactions, which would also yield a hadronic shower plus missing energy. Although the number of background events in this case is very significant, several cuts can be applied to increase the signal-to-background ratio, see Refs. [85][86][87]. Some of the most relevant signal/background discriminators include the energy of the most energetic pion, the number of pions produced in the shower, or the total visible energy of the event excluding the leading pion [87]. As a benchmark value we set the signal detection efficiency at 30% based on Ref. [86].
The second relevant aspect of ν τ detection involves energy smearing. As outlined above, part of the incident neutrino energy will be carried away by the ν τ produced in the τ decay. Therefore, the visible energy will be significantly lower than the true incident neutrino energy. We implement this effect following Ref. [86]: for a given E true ν , we assume that the observed energy distribution follows a Gaussian with mean value 0.45E true ν and width 0.25E true ν . As for the background from NC events, we use the same migration matrices provided by the experimental collaboration for the ν e CC channel. However, in this case we replace the background rejection efficiencies (which are provided for the e-like sample) by a constant 0.5%, based on Ref. [86]. As a final comment we stress that the signal and background rejection efficiencies considered here are probably somewhat conservative and could be improved with the use of more sophisticated machine learning techniques, as demonstrated in Ref. [87] for the event sample where the τ decays into electrons.

Event rates
The expected event rates are shown in Fig. 2 for the ν e -like and ν τ -like samples, as indicated by the labels. The total event rates for the different channels are summarized in Tab. 4 in App. A for the nominal and HE beam modes. For illustration purposes, we show in Fig. 2 separately the background and signal contributions, for a value of the NU parameters which saturates the present bounds. As can be seen from the figure, even in the NU scenario the signal and background event rates present a very different dependence with the observed energy. In the ν e -like sample, this is because the leading contribution to the background comes from the intrinsic contamination of the beam, which has a very different shape compared to the leading ν µ component since it stems mostly from kaon instead of pion decays. In the case of the ν τ -like sample, on the other hand, the fraction of energy carried away by the outgoing neutrino in NC scattering events and τ decays will be different, which translates into a different amount of migration towards lower values of the observed energy.
Overall, the expected event rates shown in Fig. 2 illustrate the need to go beyond a naive implementation of systematic errors in terms of normalization uncertainties, and to include shape uncertainties as well. The impact will be even larger in the sterile neutrino case (not shown here for conciseness), where the signal event rates would show an oscillatory pattern.

Systematic uncertainties included in the fit
At conventional long-baseline neutrino oscillation experiments the neutrino flux is produced from meson decays (predominantly pions and kaons) which, after being produced are focused into a decay pipe. Thus, the neutrino flux prediction cannot be computed analytically and relies on Monte Carlo simulations. Its associated systematic errors are largely dominated by hadron-production uncertainties, which can be as large as 10%-15% (see e.g. Ref. [77] in the context of DUNE). Additional flux uncertainties come from beam focusing effects and can modify the expected shape of the neutrino flux, although these are typically kept at the few percent level.
As shown in Fig. 2, the nominal flux at DUNE peaks at approximately 2.5 GeV. At these energies, there is a similar contribution to the total neutrino event rates from quasielastic (QE), resonance production (RES) and deep-inelastic scattering (DIS) processes. While the DIS cross section can be accurately described starting from the neutrino interactions with the partons in the nucleons, the situation is very different for RES and QE processes, for which nuclear effects are relevant. For example it has been shown that the impact of two-particle-two-hole (2p2h) effects and final state interactions can lead to a significant bias in the neutrino energy reconstruction [88][89][90][91][92], which in turn can have a large impact on the measurement of physics observables (see e.g. Ref. [21] for a recent review). Therefore, in what follows special attention has been paid to the role of shape uncertainties in our fits.

Event sample
Contribution Benchmark 1 Benchmark 2 Benchmark 3 σ norm σ shape σ norm σ shape σ norm σ shape Table 3. Assumed prior uncertainties affecting the normalization and shape of the event rates in our simulations, for the three benchmark scenarios considered in this work. We assume the same uncertainties (although completely uncorrelated) for the HE and nominal beam configurations, as well as for the neutrino and antineutrino modes. The background contributions are separated into intrinsic beam contamination (intrinsic cont.), flavor mis-identification (flavor mis-ID) and neutralcurrent (NC) backgrounds. We have numerically checked that for the ν e -like and ν τ -like events, including a shape uncertainty also for the signal has a completely negligible impact in the analysis, as expected. For additional details on the systematics implementation, see App. B.
which includes the same systematics implementation as in Ref. [93]. In order to account for shape uncertainties in our fit, a set of bin-to-bin uncorrelated nuisance parameters has been included for the most significant contributions to the total background. Additionally, for ν µ -like events, the sensitivity is expected to be limited by the systematics on the ν µ → ν µ channel, and therefore shape uncertainties are also included for the signal in this case. An overall nuisance parameter (bin-to-bin correlated) is also included to allow for an overall change in normalization, for each signal and background channel separately, for all channels. A pull term is then added to the χ 2 for each nuisance parameter, and the final χ 2 is obtained after minimization over all nuisance parameters included in the fit. Table 3 summarizes the prior uncertainties included in our simulations, while a more detailed description of our χ 2 implementation and the correlations implemented can be found in App. B. As can be seen from Tab. 3, the priors assumed for the normalization uncertainties range between 5% and 20%, depending on the particular channel. Our reasoning for choosing these values is based on the fact that, although some of the predictions for the fluxes and cross sections rely on the prediction from Monte Carlo simulations (with uncertainties that can be as large as 10% -20%), others may be reduced combining different measurements at the ND. For example, in the ν µ -like sample we assume a 10% normalization uncertainty, mainly driven by the large hadron production uncertainties affecting the flux prediction from simulations. However, a measurement of the ν µ CC event sample at the ND can be used to normalize the flux for the contribution to the signal from ν µ → ν e anomalous events, as well as the backgrounds from flavor mis-identification to this search. Thus, for these contributions we assume a smaller uncertainty, at the 5% level. Similarly, while the neutrino NC cross section in the few GeV range is poorly known [94], a dedicated measurement of NC events at the DUNE ND can be used to reduce the uncertainties on the NC background for the ν e -like and ν τ -like samples. In the case of shape uncertainties, driven by cross section and beam focusing effects, we consider two different benchmark values at the 2% and 5% as outlined in Tab. 3, plus a third case where shape uncertainties are not included in the analysis.

Results
This section summarizes the main results of our work. We will discuss separately the results obtained for the two frameworks under study, namely: NU of the leptonic mixing matrix coming from new physics at low scales, and sterile neutrinos participating in oscillations. For completeness, we will also provide the results for NSI in production and detection, which can be obtained from the mapping between NU and NSI as discussed in Sec. 2.

Non-unitarity
In the non-unitarity case the observable impact on the event rates at the DUNE ND would be through a change in normalization. As a consequence, the sensitivity via disappearance channels to α ee and α µµ , see Eq. (2.13), is expected to be limited by the size of the systematic errors affecting the normalization of the event rates for the ν e -like and ν µ -like samples. For this reason, we do not consider these parameters here and focus instead on the determination of the off-diagonal parameters α µe and α τ µ , which would lead to observable differences in the energy dependence of the total event rates, as shown by the histograms in Fig. 2.
Our results are shown in Fig. 3 for α µe (left panel) and α τ µ (right panel), respectively. The different lines show the results for different choices of systematic uncertainties as indicated by the labels, as a function of running time. Two main features can be seen from both panels in the figure. First, the results depend severely on the choice of systematic errors (as expected). Second, if shape uncertainties are included in the analysis the results are completely dominated by systematics and the sensitivity does not improve significantly by increasing the number of years of data taking in each mode. There is a certain improvement in sensitivity due to the changes between neutrino and antineutrino modes, as well as from the addition of extra data in the HE mode: this is due to new information coming into the fit from a different dependence of the signal-to-background ratio with the energy in the different running modes. However, after an initial (sharp) increase in sensitivity the experiment enters again a systematics-dominated regime and the results stop improving with new data.
Our results show an expected improvement with respect to current bounds for α µe , as shown in the right panel. In the case of α τ µ , however, the outcome will depend on the final level of systematic uncertainties: even in the case in which a 2% shape uncertainty is assumed the improvement over current bounds on this parameter is marginal. This is so because of the much larger background levels expected for the ν τ samples with respect to the signal (as shown in Fig. 2).

Sterile neutrinos
Before discussing our results, let us point out that the determination of the sensitivity to sterile neutrinos from oscillations is non-trivial from the statistical point of view. In particular, the results obtained under the assumption that Wilks' theorem [95] applies may differ from the ones obtained via Monte Carlo simulation [43,[96][97][98]. The sensitivity contours in this work have however been computed under the assumption that Wilks' theorem applies, in order to ease the comparison with previous literature.
In order to discuss the sensitivity to sterile neutrino oscillations it is important to keep in mind that at the DUNE ND the event rates would be sensitive to several oscillation probabilities: P ee , P µµ , P µe and P µτ . Each of them will be sensitive to a different set of mixing matrix elements as described in Sec. 2, see Eq. (2.10). In the 3 + 1 framework it is also common to parameterize the full mixing matrix as U = R 34 S 24 S 14 R 23 S 13 R 12 , where R ij is a rotation matrix in the ij sector with mixing angle θ ij , while S 14 are complex rotation matrices which also include a CP phase δ ij . The new rotation angles driving the mixing between the mostly sterile mass eigenstate and the active neutrinos can be trivially mapped to the triangular parameterization, and also be written in terms of mixing matrix elements as: .
In Fig. 4 we present the results obtained for the sensitivity to a sterile neutrino inducing non-standard P ee , P µµ and P µτ probabilities at the DUNE ND. In each case, the sensitivity comes mainly from the measurement of the ν e -like, ν µ -like and ν τ -like samples separately. Therefore, no combination between different samples is performed for any of the results shown in this figure. For comparison, we compare our results to the present dominant constraints from Daya Bay/Bugey-3 [101], NOMAD [65], CHORUS [99], MI-NOS/MINOS+ [101], and Icecube [100]. Note that Eq. (2.9) implies that the oscillation probability is invariant under the change |U α4 | 2 → (1 − |U α4 | 2 ). Therefore, a null result in the disappearance channels can be accommodated not only for small values of U α4 (as expected) but for large values as well. Very large values of the mixing matrix elements in the fourth column of the PMNS are however excluded by present neutrino oscillation data [102]. Thus, in order to be consistent with current bounds we only consider values of the mixing matrix elements such that |U µ(e)4 | 2 < 0.5.
As in the case for the NU scenario, in the sterile neutrino case we see again a large dependence of the results with the implementation of systematic uncertainties. The effect is dramatic on the sensitivity to P µµ oscillations, as shown in the middle panel of Fig. 4, where shape uncertainties at the level of 2% lead to a decrease in sensitivity of over an order of magnitude with respect to the scenario where only normalization uncertainties are considered. In this case, only a mild improvement over present constraints from MINOS could be expected at DUNE. The reason is that, in this channel, the sensitivity is almost entirely driven by the observation of a small oscillatory pattern on top of the expected (large) number of ν µ CC events (see Tab. 4 in App. A), which can be easily hidden by the inclusion of shape uncertainties. The effect is similar (although a bit milder) for the P ee case shown in the top panel of Fig. 4; however, current constraints in the ∆m 2 41 ∼ 10 eV 2 region are not as strong as in the P µµ case and, therefore, there is still room for a considerable improvement (up to one order of magnitude) if shape uncertainties are at the 2% level. Finally, we also present our results for the P µτ channel in the lower panel of Fig. 4. In this case, the sensitivity is severely limited by the reduced statistics; however, for shape uncertainties at the level of 2%, there is room for a factor ∼ 5 improvement over current bounds from CHORUS [99] and NOMAD [65] for mass-squared splittings ∆m 2 41 50 eV 2 , specially if HE data is included in the fit (indicated by the shaded blue band). Figure 4 also shows how the average-out regime (in which the sensitivity becomes independent of ∆m 2 41 ) is achieved for ∆m 2 41 O(100) eV 2 . Notice that, for the P µµ and P ee channels, the sensitivity in this region is basically limited by the normalization Figure 4. Expected sensitivity to the sterile neutrino scenario, for oscillations in the P ee , P µµ and P µτ channels for the top, middle and lower panels, respectively. The shaded regions show current constraints from other experiments [65,[99][100][101], while the colored lines indicate the sensitivity for the DUNE ND, for different choices of systematic uncertainties as indicated by the legend. In each case, the region to the right of the colored lines would be disfavored at 90% CL (2 d.o.f.). Finally, the colored blue bands indicate the increase in sensitivity due to the addition of 3.5 years of data taken in the HE mode.
uncertainty, see Eq. (2.13). Since we are using similar normalization uncertainties for both channels (see Tab. 3), a similar limit is obtained for |U e4 | 2 and |U µ4 | 2 in this regime. For the P µτ channel, we also observe that in the averaged-out region the low-scale NU limit shown in Fig. 3 is recovered, as expected from Eq. (2.12). Figure 5, on the other hand, shows our results for oscillations in the P µe channel. Again, for comparison the shaded pink, purple and gray areas show current dominant constraints from previous experiments, namely, NOMAD [64], KARMEN [103] and the combination of MINOS/MINOS+ and Daya Bay/Bugey-3 data [101]. Note that while NOMAD and KARMEN directly probed the P µe channel, the combination of MINOS/MINOS+ and Daya Bay/Bugey-3 provides an indirect constraint from P ee and P µµ disappearance channels. For reference, the shaded yellow/orange regions indicate the fraction of parameter space that is favored at 99% CL by the LSND [36] and MiniBooNE [37] anomalies.
In principle the sensitivity to the P µe channel comes from the measurement of the ν elike sample, which is directly sensitive to transitions between muon and electron neutrinos. However, at DUNE the combination with data from the ν µ -like sample would greatly enhance the sensitivity since it would allow for a simultaneous constraint on | U e4 | 2 , | U µ4 | 2 , and | U e4 | 2 | U µ4 | 2 . This can be seen from the comparison between the top and lower panels in Fig. 5: while in the top panel the analysis is performed using only information from the ν e -like sample (which is sensitive to P µe and P ee ), in the lower panel we also include the information on P µµ from the ν µ -like sample.
In fact, the combination of data from different channels for a sterile neutrino search was already performed for the DUNE experiment in Ref. [104], leading to an impressive sensitivity to sterile neutrinos which reached values of the effective mixing angle of sin 2 2θ µe 10 −6 in the region where ∆m 2 41 ∼ O(1 − 10) eV 2 . However, in Ref. [104] only normalization uncertainties were included in the analysis. Interestingly, even though our implementation of normalization uncertainties is very different than in Ref. [104] we recover a similar result here, shown by the solid blue line in the lower panel of Fig. 5. This already indicates that normalization uncertainties are not playing a relevant role in the fit, which is expected since the effect from a sterile neutrino would be detected by looking at spectral distortions in the event rates. Indeed, once shape uncertainties are included at the 2% level (shown by the dashed blue lines in the lower panel) the corresponding result is worsened by two orders of magnitude with respect to the case where only normalization uncertainties are included. Nevertheless, we find that even in this case the DUNE ND data should be able to probe most of the regions favored by the LSND and MiniBooNE experiments at more than 90% CL, with the exception of the region at very low ∆m 2 41 (which is in any case already disfavored by the combination of MINOS/MINOS+ and Daya Bay/Bugey-3 data). Furthermore, notice that DUNE has the additional advantage of being able to consistently probe the sterile neutrino hypothesis with the same experimental setup using different and complementary neutrino oscillation channels: P ee , P µµ and P µe .

Non-Standard Interactions in production and detection
Finally, for completeness we also include an explicit recasting of our bounds from low scale NU to the NSI formalism in production and detection, following the prescription from Eq. (2.15). Our results are shown in Fig. 6, for NSI parameters affecting P µe (left panel) and P µτ (right panel).  101,103], while the shaded yellow/orange regions are favored at the 99% CL by the LSND [36] and MiniBooNE [37] anomalies. The colored lines indicate the sensitivity for the DUNE ND, for different choices of systematic uncertainties as indicated by the legend. In each case, the region to the right of the lines would be disfavored at 90% CL (2 d.o.f.). Finally, the shaded band to the left of the dashed lines indicate the increase in sensitivity due to the addition of 3.5 years of data taken in the HE mode.
As can be seen from Eq. (2.14), the oscillation probability in this case depends on the production and detection NSI including their possible CP-violating phases, which come in through an interference term that is proportional to cos(Φ s βγ − Φ d βγ ) ≡ cos ∆Φ βγ . For this reason we show in both panels in Fig. 6 the sensitivity to production and detection NSI in two limiting cases, depending on the value of ∆Φ. Figure 6. Sensitivity to NSI in production and detection. Results are shown for effects on the P µe (left panel) and P µτ (right panel). In both panels the sensitivity is shown for the two limiting cases ∆Φ = π, 0, which lead to a destructive and constructive interference between production and detection NSI respectively (see Eq. (2.14)).
Our results for NSI can be compared to those obtained in Ref. [105]. Regarding the sensitivity in the disappearance channels, we find that the results are essentially limited by the normalization uncertainties (as expected from naive theoretical arguments) and are therefore not shown here. On the other hand, for the appearance channels we have checked that, when no shape uncertainties are considered, our results for the P µe and P µτ channels agree relatively well with those obtained in Ref. [105]. For the P µτ channel we have checked that including in our analysis the additional events from the τ decay channel into electrons has a negligible impact on our results.

Summary and conclusions
Oscillation experiments using conventional neutrino beams rely on the use of a near detector (ND), located at typical distances of O(500) m from the target station, to measure the unoscillated neutrino event rates. In this work we have studied the potential of the ND at the future DUNE experiment to test new physics affecting oscillations in three scenarios: effects coming from non-unitarity (NU) of the leptonic mixing matrix due to new states at low scales, light sterile neutrino oscillations, and Non-Standard Interactions (NSI) affecting neutrino production and detection processes. While sterile neutrinos would induce a new oscillation at the ND, NU and NSI would lead to a so-called zero-distance effect on the oscillation probabilities.
Using the ND to search for new physics is, however, technically challenging due to the large systematic uncertainties associated to the neutrino flux and cross sections. While most works in previous literature include only normalization uncertainties, for the scenarios considered in this work shape uncertainties are most important. Our results show that shape uncertainties at the 2% level can lead to a reduction in the sensitivity to sterile neutrino oscillations in the P µe channel of two orders of magnitude, and around one order of magnitude in the case of oscillations in the P ee and P µµ channels. Their impact is also relevant for the NU scenario, leading to a decrease in sensitivity of about a factor of 2 (3) in the case of α µe for shape uncertainties at the 2% (5%) level. Our results qualitatively agree with those presented in Ref. [24], where the sensitivity to both scenarios (including shape uncertainties) via the P ee , P µµ and P µe channels was studied (albeit for different ND configurations).
We have also studied the sensitivity to new physics leading to ν µ → ν τ transitions. In this case we find that the sensitivity is limited by the low signal-to-background ratio. This stems from the fact that the charged-current ν τ cross section starts to grow for neutrino energies E ν 3.5 GeV, above the region where the nominal DUNE neutrino flux reaches its maximum. Our results show that this is specially relevant for the NU scenario, for which the DUNE ND will only be able to improve marginally over current constraints if shape uncertainties can be kept at the 2% level or below.
Finally, we have studied the impact on the results of additional data taken with the neutrino beam in the so-called High-Energy (HE) configuration. Being much more energetic, this leads to a slightly better sensitivity to new physics effects in the P µτ channel. However, for the rest of oscillation channels considered in this work, the impact of the HE data is always rather minor since the sensitivities are mostly limited by systematic errors instead of statistics.
In summary, our results show that the potential of the DUNE ND to constrain new physics scenarios affecting oscillations at short baselines is severely limited by shape uncertainties, an effect that is typically overlooked in the literature. This stresses the importance of a joint experimental and theoretical effort to improve our understanding of neutrino nucleus cross sections, as well as hadron production uncertainties and beam focusing effects. Nevertheless, even with our more conservative and realistic implementation of systematic uncertainties, our results indicate that an improvement over current bounds is generally expected.

A Total event rates
In this work the event rates are computed following Ref. [80]. While we refer the interested reader to the mentioned references for details, here we provide a summary of the total event rates considered in our analysis, which are useful to understand the final results from our simulations.
The total event rates are summarized in Tab. 4 for the nominal and HE beam configurations. They are given separately for the contributions from: intrinsic contamination (ν e andν e ) of the beam; ν µ andν µ CC events mis-identified as ν e orν e events; and NC events mis-identified as CC events (in all cases, adding up the neutrino and antineutrino contributions to the event rates). The total number of ν µ andν µ CC events are also provided, for P µµ = 1. The number of ν e andν e CC events can be estimated from this multiplying by the corresponding P µe probability (up to small differencies in detection efficiencies for the two flavors). The total number of ν τ andν τ events is also provided for P µτ = 1, and can be trivially rescaled down for a given value of the NU parameters, see Eqs. (2.11) and (2.12). All event rates are provided after efficiencies and smearing effects are taken into account. The numbers provided correspond to the total rates within a given energy window above 0.5 GeV and below the maximum observed energy given in the last column of the table. The binning is taken as in Ref. [80].

B Details on the χ 2 implementation used
In order to account for shape uncertainties in our fit, a set of bin-to-bin uncorrelated nuisance parameters has been included in the analysis. Given the large number of bins involved in the analysis and, in order to keep the total number of nuisance parameters as low as possible, these are applied to the final background event rates (adding up all the background contributions listed in Tab. 4) and to the overall signal event rates as well (adding up the neutrino and antineutrino contributions to the signal). The shape nuisance parameters in the i-th energy bin are denoted as ξ sig i or ξ bg i for the signal and background, respectively. An independent set of nuisance parameters (bin-to-bin correlated) is also included to allow for an overall change in normalization, in this case for each signal and background channel separately (ζ c , for a given contribution c ≡ s, b to the signal or background rates). Thus, for each sample a Poissonian χ 2 is defined as where O i stands for the simulated "observed" event rates in the i-th energy bin, which we always take to be those predicted in the standard three-family framework. On the other hand, N i ({Θ, ξ, ζ}) stands for the predicted event rates in the same bin for a new physics  Table 4. Expected total number of events for the different contributions to the event samples simulated in our analysis, for the nominal and high-energy (HE) beam configurations. Note that for the nominal beam we consider 3.5 years of data taking in neutrino and in antineutrino modes, while for the HE configuration we only consider 3.5 years in neutrino mode, see Tab. 2. The background contributions are separated into intrinsic beam contamination (intrinsic cont.), flavor mis-identification (flavor mis-ID) and neutral-current (NC) backgrounds, adding the contributions from neutrinos and antineutrinos. We also provide the expected number of charged-current (CC) ν µ andν µ events for P µµ = 1, as well as the ν τ andν τ CC events for P µτ = 1. The latter can be trivially rescaled to a given value of α, see Eqs. (2.11) and (2.12). In all cases, detection efficiencies and smearing effects have already been accounted for. Finally, the last column indicates the maximum energy considered for each of the samples (corresponding to observed energy). See text for additional details.
model which depends on a set of parameters generically denoted as {Θ}. It also depends on the nuisance parameters, as: where we have included a set of possible signal (s) and background (b) contributions to the event rate sample in the bin. Also, from Eq. (B.2) it is clear that new physics effects are included on the backgrounds event rates as well (and not just on the signal). A pull term is then added to Eq. (B.1) for each nuisance parameter, and the final χ 2 is obtained after minimization over all nuisance parameters included in the fit. The final χ 2 then reads: where σ norm,c indicates the prior uncertainties assumed for the normalization of each signal or background contribution, while σ shape stands for the shape uncertainties assumed for the overall background or signal event rates in the sample considered. A series of simplifications and approximations have been performed in order to keep the number of nuisance parameters as low as possible and to avoid numerical issues during minimization. First, we assume that the shape uncertainties are the same in each energy bin and therefore have removed the index i in the corresponding priors as can be seen from Eq. (B.2) and (B.3). Next, we have assumed that all normalization uncertainties are uncorrelated between different signal and background contributions as well as between different event samples. Although in practice some degree of correlation is expected, our approach avoids numerical difficulties during minimization, which is important given the large number of nuisance parameters involved in the fit. It is also a conservative approach, as including correlations will most likely lead to an effective reduction on the size of systematic errors and a corresponding increase in sensitivity. Finally, although in the most general case one should in principle assume that both the signal and background would be affected by shape uncertainties, this would lead to a large number of nuisance parameters in the fit. Therefore, shape uncertainties are only included for the background contributions to the ν e -and ν τ -like samples, while for the ν µ -like sample we do include it both for signal and background rates. With the aim of keeping the simulations feasible, we include shape uncertainties for the sum of all background contributions in all cases, while we do include a separate normalization uncertainty for each contribution separately. In this way, we allow the different background contributions to float independently in the fit, while the shape systematics would affect the overall background rates. We believe this is enough given that in any case the background is typically dominated by a single channel (intrinsic contamination in the case of ν e andν e samples, or NC events in the case of ν µ and ν τ -like samples, see Tab. 4). The values assumed for the prior uncertainties are summarized in Tab. 3.
As shown in Sec. 4, in the case of sterile neutrino oscillations in the ν µ → ν e channel, it is convenient to add information from the ν e → ν e and ν µ → ν µ event rates as well. While information from ν e → ν e oscillations is already included in the e-like event sample (since the intrinsic ν e background event rates would be affected by oscillations), it is necessary to combine this information with the one extracted from the µ-like event rates. In this case, a separate χ 2 contribution is computed for the two samples using Eq. (B.3), and the total χ 2 is obtained after adding up the two.

C Numerical implementation of oscillations in the average-out regime
GLoBES uses a sampling method to evaluate the number of events in each bin. To avoid fast oscillations leading to aliasing, a low-pass Gaussian filter is available [78,79]. However, by default this filter does not handle well sampling points 6 which are not equally separated. This is relevant for our DUNE implementation, due to the very different bin sizes used at low and high energies. In this appendix we present a formal description of the Gaussian filter, and we describe the way it has been implemented in our code to avoid this issue.
Our goal is to perform the average of the probability within each sampling bin. We will consider that the value of the E i in each sampling bin follows a normal distribution with mean given by the central value of the bin (E c i ) and standard deviation given by the size of the bin (∆E i ). The energy enters in the oscillatory part of the probability as sin 2 ( ∆m 2 41 L 4E ), therefore to perform the average we need to know the probability distribution of L E . This leads us to the next formal problem.
Let us consider two independent random variables, normally distributed, X ∼ N µ x , σ 2 x and Y ∼ N µ y , σ 2 y . Then the ratio Z = X Y follows the non-trivial distribution [106]: However, when δ y < 0.1 the density distribution in Eq. (C.1) also follows a normal distribution to a very good approximation, namely X Y ∼ N µ x /µ y , σ 2 X/Y , where: (C. 6) In our case, Y = L and X = E i , and therefore µ x = L, µ y = E c i , σ x = 0, σ y = ∆E i . The condition that allow us to approximate the distribution of L/E i by a normal distribution is now given by ∆E i < 0.1E c i , implying that the choice of the sampling bins cannot be arbitrary if we want to use a Gaussian filter. Under these conditions the density distribution for L/E i is given by: (C.7) Therefore in each sampling bin the average probability is given by: In fact, the result from this expression can be computed analytically [107], yielding the following result for the case of the sterile oscillation probability at the ND: which is the expression we use for our sterile neutrino analysis with GLoBES.