1 Introduction

The Deep Underground Neutrino Experiment (DUNE) is a next-generation, long-baseline (LBL) neutrino oscillation experiment, designed to be sensitive to \(\nu _\mu \) to \(\nu _e \) oscillation. The experiment consists of a high-power, broadband neutrino beam, a powerful precision near detector (ND) complex located at Fermi National Accelerator Laboratory, in Batavia, Illinois, USA, and a massive liquid argon time-projection chamber (LArTPC) far detector (FD) located at the 4850 ft level of Sanford Underground Research Facility (SURF), in Lead, South Dakota, USA. The baseline of 1285 km provides sensitivity, in a single experiment, to all parameters governing LBL neutrino oscillation. The deep underground location of the FD facilitates sensitivity to nucleon decay and other rare processes including low-energy neutrino detection enabling, for instance, observation of neutrinos from a core-collapse supernova.

Owing to the high-power proton beam facility, the ND consisting of precision detectors capable of off-axis data taking and the massive FD, DUNE provides enormous opportunities to probe phenomena beyond the SM traditionally difficult to reach in neutrino experiments. Of such vast, rich physics topics that profoundly expand those probed in the past neutrino experiments, this paper reports a selection of studies of DUNE’s sensitivity to a variety of BSM particles and effects, initially presented in the physics volume of the DUNE Technical Design Report (TDR) [1] recently made available. Some of these phenomena impact the LBL oscillation measurement, while others may be detected by DUNE using specific analyses.

Section 2 describes some of the common assumptions and tools used in these analyses. Section 3 discusses sensitivity to sterile neutrinos, Sect. 4 looks into the effect of non-unitary of the neutrino mixing matrix, Sect. 5 describes sensitivity to non-standard neutrino interactions, Sect. 6 discusses sensitivity to CPT and Lorentz violation, Sect. 7 describes the sensitivity to new physics by measuring neutrino trident production, Sect. 8 discusses various dark matter searches that could be performed by DUNE, Sect. 9 describes sensitivity to baryon number violation by one and two units, and Sect. 10 lists some other possible avenues for BSM physics searches.

These studies reveal that DUNE can probe a rich and diverse BSM phenomenology at the discovery level, as in the case of searches for dark matter created in the high-power proton beam interactions and from cosmogenic sources, or by significantly improving existing constraints, as in the cases of sterile neutrino mixing, non-standard neutrino interactions, CPT violation, new physics enhancing neutrino trident production, and nucleon decay.

2 Analysis details

The BSM searches presented in this paper span a wide variety of physics topics and techniques. The analyses rely on neutrino beam data taken at the ND and/or FD, atmospheric or other astrophysical sources of neutrinos, or signal from the detector material itself, as in nucleon decay searches. This section summarizes some of the common assumptions and tools used in the analyses, with more details provided in the following sections.

2.1 Detector assumptions

The DUNE FD will consist of four \(10 \, \text {kt}\) fiducial mass LArTPC modules with integrated photon detection systems (PD systems) [2,3,4]. In these analyses, we assume all four modules have identical responses. All of the analyses described will use data from the FD, except for the analyses presented in Sects. 78.1, and 10.3, which use data exclusively from the ND.

The ND will be located at a distance of \(574 \, \text {m}\) from the target. The ND concept consists of a modular LArTPC, a magnetized high-pressure gas argon TPC and a beam monitor. The combination of the first two detectors is planned to be movable to sample the off-axis neutrino spectrum to reduce flux uncertainties, a concept called DUNE-PRISM [1]. Since the ND configuration, however, was not yet finalized at the time these studies were performed, we adopted only the LArTPC component of the detector and its fiducial volume. In the analyses presented here, the LArTPC is assumed to be \({7} \, \text {m}\) wide, \({3} \, \text {m}\) high, and \({5} \, \text {m}\) long. The fiducial volume is assumed to include the detector volume up to 50 cm of each face of the detector. The ND properties are given in Table 1. The signal and background efficiencies vary with the physics model being studied. Detailed signal and background efficiencies for each physics topic are discussed along with each analysis.

Table 1 LArTPC ND properties used in some of the BSM physics analyses

2.2 Neutrino beam assumptions

The analyses described in Sects. 345, and 6 are based on analysis of neutrino beam data at both the ND and FD. The DUNE neutrino beam is produced using protons from Fermilab’s Main Injector and a traditional horn-focusing system [5]. The polarity of the focusing magnets may be reversed to produce a neutrino- or antineutrino-dominated beam. This optimized beam configuration includes a three-horn focusing system with a 1 m long target embedded within the first horn and a decay pipe with \(194 \, \text {m}\) length and \(4 \, \text {m}\) diameter. The neutrino flux produced by this beamline is simulated at a distance of \(574 \, \text {m}\) downstream of the neutrino target for the ND and \(1285 \, \text {km}\) for the FD. Fluxes have been generated for both neutrino mode and antineutrino mode using G4LBNF [1, 6], a Geant4-based simulation [7,8,9].

Results based on beam neutrino data are given for a \(300~\text {kt} \, \cdot \, \text {MW} \cdot \, \text {year} \) exposure. With the current deployment plan [1], this exposure will be achieved in approximately 7 years once the beam is operational. For results not based on beam data, the exposure is given in units of \(\text {kt} \, \cdot \, \text {year}\) in each relevant section.

2.3 Tools

In the analyses presented in Sects. 345, and 6, the simulation of the DUNE experimental setup was performed with the General Long-Baseline Experiment Simulator (GLoBES) software [10, 11]. Unless otherwise noted, the neutrino fluxes used in the BSM physics analysis are the same as those used in the DUNE LBL three-flavor analysis [1]. The configuration of the beam used in ND analyses is assumed to be a 120 GeV proton beam with 1.2 MW beam power at 56% uptime, providing \(1.1\times 10^{21}\) POT/year. Cross-section files describing neutral current (NC) and charged current (CC) interactions with argon are generated using Generates Events for Neutrino Interaction Experiments (GENIE) [12, 13] version 2.8.4. The true-to-reconstructed smearing matrices and the selection efficiency as a function of energy for various signal and background modes are generated using nominal DUNE MC simulation. A \(40 \, \text {kt}\) fiducial mass is assumed for the FD, exposed to a \(120 \, \text {GeV}\), \(1.2 \, \text {MW}\) beam. The \(\nu _{e}\) and \({{\bar{\nu }}}_{e}\) appearance signal modes have independent normalization uncertainties of \(2\%\) each, while \(\nu _{\mu }\) and \({\bar{\nu }}_{\mu }\) disappearance signal modes have independent normalization uncertainties of \(5\%\). The background normalization uncertainties range from 5 to \(20\%\) and include correlations among various sources of background. More details can be found in Ref. [1].

The neutrino trident search presented in Sect. 7 and the baryon number violation analyses presented in Sect. 9 use samples of simulated and reconstructed signal and background events, produced using standard DUNE detection simulation and reconstruction software. Further details are given in those sections.

For analyses that use neither GLoBES nor the standard DUNE simulation and reconstruction software, such as the dark matter analyses described in Sect. 8 and several of the analyses described in Sect. 10, details are given in the relevant sections.

3 Sterile Neutrino Mixing

Experimental results in tension with the three-neutrino-flavor paradigm, which may be interpreted as mixing between the known active neutrinos and one or more sterile states, have led to a rich and diverse program of searches for oscillations into sterile neutrinos [14, 15]. DUNE is sensitive over a broad range of potential sterile neutrino mass splittings by looking for disappearance of CC and NC interactions over the long distance separating the ND and FD, as well as over the short baseline of the ND. With a longer baseline, a more intense beam, and a high-resolution large-mass FD, compared to previous experiments, DUNE provides a unique opportunity to improve significantly on the sensitivities of the existing probes, and greatly enhance the ability to map the extended parameter space if a sterile neutrino is discovered. In the sterile neutrino mixing studies presented here, we assume a minimal 3+1 oscillation scenario with three active neutrinos and one sterile neutrino, which includes a new independent neutrino mass-squared difference, \({\varDelta }m^2_{41}\), and for which the mixing matrix is extended with three new mixing angles, \(\theta _{14}\), \(\theta _{24}\), \(\theta _{34}\), and two additional phases \(\delta _{14}\) and \(\delta _{24}\).

Disappearance of the beam neutrino flux between the ND and FD results from the quadratic suppression of the sterile mixing angle measured in appearance experiments, \(\theta _{\mu e}\), with respect to its disappearance counterparts, \(\theta _{\mu \mu }\approx \theta _{24}\) for LBL experiments, and \(\theta _{ee}\approx \theta _{14}\) for reactor experiments. These disappearance effects have not yet been observed and are in tension with appearance results [14, 15] when global fits of all available data are carried out. The exposure of DUNE’s high-resolution FD to the high-intensity LBNF beam will also allow direct probes of non-standard electron (anti)neutrino appearance.

DUNE will look for active-to-sterile neutrino mixing using the reconstructed energy spectra of both NC and CC neutrino interactions in the FD, and their comparison to the extrapolated predictions from the ND measurement. Since NC cross sections and interaction topologies are the same for all three active neutrino flavors, the NC spectrum is insensitive to standard neutrino mixing. However, should there be oscillations into a fourth light neutrino, an energy-dependent depletion of the neutrino flux would be observed at the FD, as the sterile neutrino would not interact in the detector volume. Furthermore, if sterile neutrino mixing is driven by a large mass-square difference \({\varDelta }m^2_{{\mathrm{41}}} \, \sim 1\, {\text {eV}}^{2}\), the CC spectrum will be distorted at energies higher than the energy corresponding to the standard oscillation maximum. Therefore, CC disappearance is also a powerful probe of sterile neutrino mixing at long baselines.

We assume the mixing matrix augmented with one sterile state is parameterized by \(U=R_{34}S_{24}S_{14}R_{23}S_{13}R_{12}\) [16], where \(R_{ij}\) is the rotational matrix for the mixing angle \(\theta _{ij}\), and \(S_{ij}\) represents a complex rotation by the mixing angle \(\theta _{ij}\) and the CP-violating phase \(\delta _{ij}\). At long baselines the NC disappearance probability to first order for small mixing angles is then approximated by:

$$\begin{aligned} 1 - P(\nu _{\mu } \rightarrow \nu _s)&\approx 1 - \cos ^4\theta _{14}\cos ^2\theta _{34}\sin ^{2}2\theta _{24}\sin ^2{\varDelta }_{41} \nonumber \\&\quad - \sin ^2\theta _{34}\sin ^22\theta _{23}\sin ^2{\varDelta }_{31} \nonumber \\&\quad + \frac{1}{2}\sin \delta _{24}\sin \theta _{24}\sin 2\theta _{23}\sin {\varDelta }_{31}, \end{aligned}$$
(1)

where \({\varDelta }_{ji} = \frac{{\varDelta }m^2_{ji}L}{4E}\). The relevant oscillation probability for \(\nu _\mu \)  CC disappearance is the \(\nu _\mu \)  survival probability, similarly approximated by:

$$\begin{aligned} P(\nu _{\mu } \rightarrow \nu _{\mu })&\approx 1 - \sin ^22\theta _{23}\sin ^2{\varDelta }_{31} \nonumber \\&\quad + 2\sin ^22\theta _{23}\sin ^2\theta _{24}\sin ^2{\varDelta }_{31} \nonumber \\&\quad - \sin ^22\theta _{24}\sin ^2{\varDelta }_{41}. \end{aligned}$$
(2)

Finally, the disappearance of \(\overset{(-)}{\nu }_e\) CC is described by:

$$\begin{aligned} P(\overset{(-)}{\nu }_e \rightarrow \overset{(-)}{\nu }_e)&\approx 1 - \sin ^22\theta _{13}\sin ^2{\varDelta }_{31} \nonumber \\&\quad - \sin ^22\theta _{14}\sin ^2{\varDelta }_{41}. \end{aligned}$$
(3)

Figure 1 shows how the standard three-flavor oscillation probability is distorted at neutrino energies above the standard oscillation peak when oscillations into sterile neutrinos are included.

Fig. 1
figure 1

Regions of L/E probed by the DUNE detector compared to 3-flavor and \(3+1\)-flavor neutrino disappearance and appearance probabilities. The gray-shaded areas show the range of true neutrino energies probed by the ND and FD. The top axis shows true neutrino energy, increasing from right to left. The top plot shows the probabilities assuming mixing with one sterile neutrino with \({\varDelta }m^2_{\mathrm{41}}=0.05 \, {\text {eV}}^2\), corresponding to the slow oscillations regime. The middle plot assumes mixing with one sterile neutrino with \({\varDelta }m^2_{\mathrm{41}}=0.5 \, {\text {eV}}^2\), corresponding to the intermediate oscillations regime. The bottom plot includes mixing with one sterile neutrino with \({\varDelta }m^2_{\mathrm{41}}=50 \, {\text {eV}}^2\), corresponding to the rapid oscillations regime. As an example, the slow sterile oscillations cause visible distortions in the three-flavor \(\nu _\mu \)  survival probability (blue curve) for neutrino energies \(\sim 10 \, {\text {GeV}}\), well above the three-flavor oscillation minimum

The sterile neutrino effects have been implemented in GLoBES via the existing plug-in for sterile neutrinos and non-standard interactions [17]. As described above, the ND will play a very important role in the sensitivity to sterile neutrinos both directly, for rapid oscillations with \({\varDelta }m_{41}^2 > 1 \, {\text {eV}}^2\) where the sterile oscillation matches the ND baseline, and indirectly, at smaller values of \({\varDelta }m_{41}^2\) where the ND is crucial to reduce the systematic uncertainties affecting the FD to increase its sensitivity. To include these ND effects in these studies, the most recent GLoBES DUNE configuration files describing the FD were modified by adding a ND with correlated systematic errors with the FD. As a first approximation, the ND is assumed to be an identical scaled-down version of the TDR FD, with identical efficiencies, backgrounds and energy reconstruction. The systematic uncertainties originally defined in the GLoBES DUNE conceptual design report (CDR) configuration already took into account the effect of the ND constraint. Thus, since we are now explicitly simulating the ND, larger uncertainties have been adopted but partially correlated between the different channels in the ND and FD, so that their impact is reduced by the combination of both data sets. The full set of systematic uncertainties employed in the sterile neutrino studies is listed in Table 2.

Table 2 List of systematic errors assumed in the sterile neutrino studies

Finally, for oscillations observed at the ND, the uncertainty on the production point of the neutrinos can play an important role. We have included an additional \(20\%\) energy smearing, which produces a similar effect given the L/E dependence of oscillations. We implemented this smearing in the ND through multiplication of the migration matrices provided with the GLoBES files by an additional matrix with the \(20\%\) energy smearing obtained by integrating the Gaussian

$$\begin{aligned} R^c(E,E')\equiv \frac{1}{\sigma (E)\sqrt{2\pi }}e^{-\frac{(E-E')^2}{2(\sigma (E))^2}}, \end{aligned}$$
(4)

with \(\sigma (E)=0.2 E\) in reconstructed energy \(E'\), where E is the true neutrino energy from simulation.

By default, GLoBES treats all systematic uncertainties included in the fit as normalization shifts. However, depending on the value of \({\varDelta }m^2_{41}\), sterile mixing will induce shape distortions in the measured energy spectrum beyond simple normalization shifts. As a consequence, shape uncertainties are very relevant for sterile neutrino searches, particularly in regions of parameter space where the ND, with virtually infinite statistics, has a dominant contribution. The correct inclusion of systematic uncertainties affecting the shape of the energy spectrum in the two-detector fit GLoBES framework used for this analysis posed technical and computational challenges beyond the scope of the study. Therefore, for each limit plot, we present two limits bracketing the expected DUNE sensitivity limit, namely: the black limit line, a best-case scenario, where only normalization shifts are considered in a ND + FD fit, where the ND statistics and shape have the strongest impact; and the grey limit line, corresponding to a worst-case scenario where only the FD is considered in the fit, together with a rate constraint from the ND.

Studying the sensitivity to \(\theta _{14}\), the dominant channels are those regarding \(\nu _e\) disappearance. Therefore, only the \(\nu _e\) CC sample is analyzed and the channels for NC and \(\nu _{\mu }\) CC disappearance are not taken into account, as they do not influence greatly the sensitivity and they slow down the simulations. The sensitivity at the 90% confidence level (CL), taking into account the systematic uncertainties mentioned above, is shown in Fig. 2, along with a comparison to current constraints.

For the \(\theta _{24}\) mixing angle, we analyze the \(\nu _{\mu }\) CC disappearance and the NC samples, which are the main contributors to the sensitivity. The results are shown in Fig. 2, along with comparisons with present constraints.

Fig. 2
figure 2

The top plot shows the DUNE sensitivities to \(\theta _{14}\) from the \(\nu _e\) CC samples at the ND and FD, along with a comparison with the combined reactor result from Daya Bay and Bugey-3. The bottom plot is adapted from Ref. [18] and displays sensitivities to \(\theta _{24}\) using the \(\nu _\mu \) CC and NC samples at both detectors, along with a comparison with previous and existing experiments. In both cases, regions to the right of the contours are excluded

In the case of the \(\theta _{34}\) mixing angle, we look for disappearance in the NC sample, the only contributor to this sensitivity. The results are shown in Fig. 3. Further, a comparison with previous experiments sensitive to \(\nu _\mu \), \(\nu _\tau \)  mixing with large mass-squared splitting is possible by considering an effective mixing angle \(\theta _{\mu \tau }\), such that \(\sin ^2{2\theta _{\mu \tau }}\equiv 4|U_{\tau 4}|^2|U_{\mu 4}|^2=\cos ^4\theta _{14}\sin ^22\theta _{24}\sin ^2\theta _{34}\), and assuming conservatively that \(\cos ^4\theta _{14}=1\), and \(\sin ^22\theta _{24}=1\). This comparison with previous experiments is also shown in Fig. 3. The sensitivity to \(\theta _{34}\) is largely independent of \({\varDelta }m^2_{41}\), since the term with \(\sin ^2\theta _{34}\) in Eq. (1), the expression describing \({P(\nu _{\mu }\rightarrow \nu _s)}\), depends solely on the \({\varDelta }m^2_{31}\) mass splitting.

Another quantitative comparison of our results for \(\theta _{24}\) and \(\theta _{34}\) with existing constraints can be made for projected upper limits on the sterile mixing angles assuming no evidence for sterile oscillations is found, and picking the value of \({\varDelta }m^2_{41} = 0.5 \, {\text {eV}}^2\) corresponding to the simpler counting experiment regime. For the \(3+1\) model, upper limits of \(\theta _{24}\, < \, 1.8^{\circ } \, (15.1^{\circ })\) and \(\theta _{34}\, < \, 15.0^{\circ } \, (25.5^{\circ })\) are obtained at the 90% CL from the presented best(worst)-case scenario DUNE sensitivities. If expressed in terms of the relevant matrix elements

$$\begin{aligned} \begin{aligned} |U_{\mu 4}|^2 =&\,\,\cos ^2\theta _{14}\sin ^2\theta _{24} \\ |U_{\tau 4}|^2=&\,\,\cos ^2\theta _{14}\cos ^2\theta _{24}\sin ^2\theta _{34}, \end{aligned} \end{aligned}$$
(5)

these limits become \(|U_{\mu 4}|^{2} \, < \,0.001\) (0.068) and \(|U_{\tau 4}|^{2} \, < \,0.067\) (0.186) at the 90% CL, where we conservatively assume \(\cos ^2\theta _{14} \,=\,1\) in both cases, and additionally \(\cos ^2\theta _{24} \,=\,1\) in the second case.

Finally, sensitivity to the \(\theta _{\mu e}\) effective mixing angle, defined as \(\sin ^2{2\theta _{\mu e}}\equiv 4|U_{e4}|^2|U_{\mu 4}|^2=\sin ^22\theta _{14}\sin ^2\theta _{24}\), is shown in Fig. 4, which also displays a comparison with the allowed regions from the Liquid Scintillator Neutrino Detector (LSND) and MiniBooNE, as well as with present constraints and projected constraints from the Fermilab Short-Baseline Neutrino (SBN) program.

Fig. 3
figure 3

Comparison of the DUNE sensitivity to \(\theta _{34}\) using the NC samples at the ND and FD with previous and existing experiments. Regions to the right of the contour are excluded

Fig. 4
figure 4

DUNE sensitivities to \(\theta _{\mu e}\) from the appearance and disappearance samples at the ND and FD are shown on the top plot, along with a comparison with previous existing experiments and the sensitivity from the future SBN program. Regions to the right of the DUNE contours are excluded. The plot is adapted from Ref. [18]. In the bottom plot, the ellipse displays the DUNE discovery potential assuming \(\theta _{\mu e}\) and \({\varDelta }m_{41}^2\) set at the best-fit point determined by LSND [19] (represented by the star) for the best-case scenario referenced in the text

As an illustration, Fig. 4 also shows DUNE’s discovery potential for a scenario with one sterile neutrino governed by the LSND best-fit parameters:

\(\left( {\varDelta }m_{41}^2= 1.2\;\text {eV}^2;\,\,\sin ^2{2\theta _{\mu e}}=0.003\right) \) [19]. A small 90% CL allowed region is obtained, which can be compared with the LSND allowed region in the same figure.

4 Non-unitarity of the neutrino mixing matrix

A generic characteristic of most models explaining the neutrino mass pattern is the presence of heavy neutrino states, additional to the three light states of the SM of particle physics [20,21,22]. These types of models imply that the \(3 \times 3\) Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix is not unitary due to mixing with additional states. Besides the type-I seesaw mechanism [23,24,25,26], different low-scale seesaw models include right-handed neutrinos that are relatively not-so-heavy, with mass of 1–10 TeV [27], and perhaps detectable at collider experiments.

These additional heavy leptons would mix with the light neutrino states and, as a result, the complete unitary mixing matrix would be a squared \(n \times n\) matrix, with n the total number of neutrino states. Therefore, the usual \(3 \times 3\) PMNS matrix, which we dub N to stress its non-standard nature, will be non-unitary. One possible general way to parameterize these unitarity deviations in N is through a triangular matrix [28]Footnote 1

$$\begin{aligned} N = \left\lgroup \begin{array}{ccc} 1-\alpha _{ee} &{} 0 &{} 0 \\ \alpha _{\mu e} &{} 1-\alpha _{\mu \mu } &{} 0 \\ \alpha _{\tau e} &{} \alpha _{\tau \mu } &{} 1-\alpha _{\tau \tau } \end{array} \right\rgroup U , \end{aligned}$$
(6)

with U representing the unitary PMNS matrix, and the \(\alpha _{ij}\) representing the non-unitary parameters.Footnote 2 In the limit where \(\alpha _{ij}=0\), N becomes the usual PMNS mixing matrix.

The triangular matrix in this equation accounts for the non-unitarity of the \(3 \times 3\) matrix for any number of extra neutrino species. This parameterization has been shown to be particularly well-suited for oscillation searches [28, 31] since, compared to other alternatives, it minimizes the departures of its unitary component U from the mixing angles that are directly measured in neutrino oscillation experiments when unitarity is assumed.

The phenomenological implications of a non-unitary leptonic mixing matrix have been extensively studied in flavor and electroweak precision observables as well as in the neutrino oscillation phenomenon [26, 28, 32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52]. For recent global fits to all flavor and electroweak precision data summarizing present bounds on non-unitarity see Refs. [46, 53].

Recent studies have shown that DUNE can constrain the non-unitarity parameters [31, 52]. The summary of the \(90 \%\) CL bounds on the different \(\alpha _{ij}\) elements profiled over all other parameters is given in Table 3.

Table 3 Expected \(90\%\) CL constraints on the non-unitarity parameters \(\alpha \) from DUNE

These bounds are comparable with other constraints from present oscillation experiments, although they are not competitive with those obtained from flavor and electroweak precision data. For this analysis, and those presented below, we have used the GLoBES software [10, 11] with the DUNE TDR configuration presented in Ref. [1] and assumed a data exposure of \(300~\text {kt} \, \cdot \, \text {MW} \cdot \, \text {year} \). The standard (unitary) oscillation parameters have also been treated as in [1]. The unitarity deviations have been included both by an independent code (used to obtain the results shown in Ref. [52]) and via the Monte Carlo Utility Based Experiment Simulator (MonteCUBES) [54] plug-in to cross validate our results.

Fig. 5
figure 5

The impact of non-unitarity on the DUNE CPV discovery potential. See the text for details

Fig. 6
figure 6

Expected frequentist allowed regions at the \(1 \sigma \), \(90\%\) and \(2\sigma \) CL for DUNE. All new physics parameters are assumed to be zero so as to obtain the expected non-unitarity sensitivities. A value \(\theta _{23} = 0.235 \pi \approx 0.738\) rad is assumed. The solid lines correspond to the analysis of DUNE data alone, while the dashed lines include the present constraints on non-unitarity. The values of \(\theta _{23}\) are shown in radians

Conversely, the presence of non-unitarity may affect the determination of the Dirac charge parity (CP)-violating phase \(\delta _{CP}\) in LBL experiments [50, 52, 53]. Indeed, when allowing for unitarity deviations, the expected CP discovery potential for DUNE could be significantly reduced. However, the situation is alleviated when a combined analysis with the constraints on non-unitarity from other experiments is considered. This is illustrated in Fig. 5. In the left panel, the discovery potential for charge-parity symmetry violation (CPV) is computed when the non-unitarity parameters introduced in Eq. (6) are allowed in the fit. While for the Asimov data all \(\alpha _{ij}=0\), the non-unitary parameters are allowed to vary in the fit with \(1 \sigma \) priors of \(10^{-1}\), \(10^{-2}\) and \(10^{-3}\) for the dotted green, dashed blue and solid black lines respectively. For the dot-dashed red line no prior information on the non-unitarity parameters has been assumed. As can be observed, without additional priors on the non-unitarity parameters, the capabilities of DUNE to discover CPV from \(\delta _{CP}\) would be seriously compromised [52]. However, with priors of order \(10^{-2}\) matching the present constraints from other neutrino oscillation experiments [31, 52], the sensitivity expected in the three-flavor model is almost recovered. If the more stringent priors of order \(10^{-3}\) stemming from flavor and electroweak precision observables are added [46, 53], the standard sensitivity is obtained.

The right panel of Fig. 5 concentrates on the impact of the phase of the element \(\alpha _{\mu e}\) in the discovery potential of CPV from \(\delta _{CP}\), since this element has a very important impact in the \(\nu _e\) appearance channel. In this plot the modulus of \(\alpha _{ee}\), \(\alpha _{\mu \mu }\) and \(\alpha _{\mu e}\) have been fixed to \(10^{-1}\), \(10^{-2}\), \(10^{-3}\) and 0 for the dot-dashed red, dotted green, dashed blue and solid black lines respectively. All other non-unitarity parameters have been set to zero and the phase of \(\alpha _{\mu e}\) has been allowed to vary both in the fit and in the Asimov data, showing the most conservative curve obtained. As for the right panel, it can be seen that a strong deterioration of the CP discovery potential could be induced by the phase of \(\alpha _{\mu e}\) (see Ref. [52]). However, for unitarity deviations of order \(10^{-2}\), as required by present neutrino oscillation data constraints, the effect is not too significant in the range of \(\delta _{CP}\) for which a \(3 \sigma \) exclusion of CP conservation would be possible and it becomes negligible if the stronger \(10^{-3}\) constraints from flavor and electroweak precision data are taken into account.

Similarly, the presence of non-unitarity worsens degeneracies involving \(\theta _{23}\), making the determination of the octant or even its maximality challenging. This situation is shown in Fig. 6 where an input value of \(\theta _{23} = 42.3^\circ \) was assumed. As can be seen, the fit in presence of non-unitarity (solid lines) introduces degeneracies for the wrong octant and even for maximal mixing [31]. However, these degeneracies are resolved upon the inclusion of present priors on the non-unitarity parameters from other oscillation data (dashed lines) and a clean determination of the standard oscillation parameters following DUNE expectations is again recovered.

Fig. 7
figure 7

Allowed regions of the non-standard oscillation parameters in which we see important degeneracies (top) and the complex non-diagonal ones (bottom). We conduct the analysis considering all the NSI parameters as non-negligible. The sensitivity regions are for 68% CL [red line (left)], 90% CL [green dashed line (middle)], and 95% CL [blue dotted line (right)]. Current bounds are taken from [78]

The sensitivity that DUNE would provide to the non-unitarity parameters is comparable to that from present oscillation experiments, while not competitive to that from flavor and electroweak precision observables, which are roughly an order of magnitude more stringent. On the other hand, the capability of DUNE to determine the standard oscillation parameters such as CPV from \(\delta _{CP}\) or the octant or maximality of \(\theta _{23}\) would be seriously compromised by unitarity deviations in the PMNS matrix. This negative impact is however significantly reduced when priors on the size of these deviations from other oscillation experiments are considered, and disappears altogether if the more stringent constraints from flavor and electroweak precision data are added instead.

5 Non-standard neutrino interactions

Non-standard neutrino interactions (NSI), affecting neutrino propagation through the Earth, can significantly modify the data to be collected by DUNE as long as the new physics parameters are large enough [55]. Leveraging its very long baseline and wide-band beam, DUNE is uniquely sensitive to these probes. NSI may impact the determination of current unknowns such as CPV [56, 57], mass hierarchy [58, 59] and octant of \(\theta _{23}\) [60]. If the DUNE data are consistent with the standard oscillation for three massive neutrinos, off-diagonal NC NSI effects of order 0.1 \(G_F\) can be ruled out at the 68 to 95% CL [61, 62]. We note that DUNE might improve current constraints on \(|\epsilon ^m_{e \tau }|\) and \(|\epsilon ^m_{e \mu }|\), the electron flavor-changing NSI intensity parameters (see Eq. 8), by a factor 2-5 [55, 63, 64]. New CC interactions can also lead to modifications in the production, at the beam source, and the detection of neutrinos. The findings on source and detector NSI studies at DUNE are presented in [65, 66], in which DUNE does not have sensitivity to discover or to improve bounds on source/detector NSI. In particular, the simultaneous impact on the measurement of \(\delta _{\mathrm{CP}}\) and \(\theta _{23}\) is investigated in detail. Depending on the assumptions, such as the use of the ND and whether NSI at production and detection are the same, the impact of source/detector NSI at DUNE may be relevant. We focus our attention on the propagation, based on the results from [65].

NC NSI can be understood as non-standard matter effects that are visible only in an FD at a sufficiently long baseline. They can be parameterized as new contributions to the matter potential in the Mikheyev–Smirnov–Wolfenstein effect (MSW) [67,68,69,70,71,72] matrix in the neutrino-propagation Hamiltonian:

$$\begin{aligned} H = U \left( \begin{array}{ccc} 0 &{} &{} \\ &{} {\varDelta }m_{21}^2/2E &{} \\ &{} &{} {\varDelta }m_{31}^2/2E \end{array} \right) U^\dag + {\tilde{V}}_{\mathrm{MSW}} , \end{aligned}$$
(7)

with

$$\begin{aligned} {\tilde{V}}_{\mathrm{MSW}} = \sqrt{2} G_F N_e \left( \begin{array}{ccc} 1 + \epsilon ^m_{ee} &{} \epsilon ^m_{e\mu } &{} \epsilon ^m_{e\tau } \\ \epsilon ^{m*}_{e\mu } &{} \epsilon ^m_{\mu \mu } &{} \epsilon ^m_{\mu \tau } \\ \epsilon ^{m*}_{e\tau } &{} \epsilon ^{m*}_{\mu \tau } &{} \epsilon ^m_{\tau \tau } \end{array} \right) \end{aligned}$$
(8)

Here, U is the standard PMNS leptonic mixing matrix, for which we use the standard parameterization found, e.g., in [73], and the \(\epsilon \)-parameters give the magnitude of the NSI relative to standard weak interactions. For new physics scales of a few hundred GeV, a value of \(|\epsilon |\) of the order 0.01 or less is expected [74,75,76]. The DUNE baseline provides an advantage in the detection of NSI relative to existing beam-based experiments with shorter baselines. Only atmospheric-neutrino experiments have longer baselines, but the sensitivity of these experiments to NSI is limited by systematic effects [77].

In this analysis, we use GLoBES with the MonteCUBES C library, a plugin that replaces the deterministic GLoBES minimizer by a Markov Chain Monte Carlo (MCMC) method that is able to handle higher dimensional parameter spaces. In the simulations we use the configuration for the DUNE TDR [1]. Each point scanned by the MCMC is stored and a frequentist \(\chi ^2\) analysis is performed with the results. The analysis assumes an exposure of \(300~\text {kt} \, \cdot \, \text {MW} \cdot \, \text {year} \).

In an analysis with all the NSI parameters free to vary, we obtain the sensitivity regions in Fig. 7. We omit the superscript m that appears in Eq. (8). The credible regions are shown for different confidence levels. We note, however, that constraints on \(\epsilon _{\tau \tau }-\epsilon _{\mu \mu }\) coming from global fit analysis [55, 64, 78, 79] can remove the left and right solutions of \(\epsilon _{\tau \tau }-\epsilon _{\mu \mu }\) in Fig. 7.

In order to constrain the standard oscillation parameters when NSI are present, we use the fit for three-neutrino mixing from [78] and implement prior constraints to restrict the region sampled by the MCMC. The sampling of the parameter space is explained in [62] and the priors that we use can be found in Table 4.

Table 4 Oscillation parameters and priors implemented in MCMC for calculation of Fig. 7

The effects of NSI on the measurements of the standard oscillation parameters at DUNE are explicit in Fig. 8, where we superpose the allowed regions with non-negligible NSI and the standard-only credible regions at 90% CL. In the blue filled areas we assume only standard oscillation. In the regions delimited by the red, black dashed, and green dotted lines we constrain standard oscillation parameters allowing NSI to vary freely.

An important degeneracy appears in the measurement of the mixing angle \(\theta _{23}\). Notice that this degeneracy appears because of the constraints obtained for \(\epsilon _{\tau \tau }-\epsilon _{\mu \mu }\) shown in Fig. 7. We also see that the sensitivity of the CP phase is strongly affected.

Fig. 8
figure 8

Projections of the standard oscillation parameters with nonzero NSI. The sensitivity regions are for 68, 90, and 95% CL. The allowed regions considering negligible NSI (standard oscillation (SO) at 90% CL) are superposed to the SO + NSI

The effects of matter density variation and its average along the beam path from Fermilab to SURF were studied considering the standard neutrino oscillation framework with three flavors [80, 81]. In order to obtain the results of Figs. 7 and 8, we use a high-precision calculation for the baseline of \(1285 \, \text {km}\) and the average density of \(2.848 \, {\text {g/cm}}^{3}\) [80].

The DUNE collaboration has been using the so-called PREM [82, 83] density profile to consider matter density variation. With this assumption, the neutrino beam crosses a few constant density layers. However, a more detailed density map is available for the USA with more than 50 layers and \(0.25 \times 0.25\) degree cells of latitude and longitude: The Shen–Ritzwoller or S.R. profile [80, 84]. Comparing the S.R. with the PREM profiles, Ref. [81] shows that in the standard oscillation paradigm, DUNE is not highly sensitive to the density profile and that the only oscillation parameter with its measurement slightly impacted by the average density true value is \(\delta _{\mathrm{CP}}\) . NSI, however, may be sensitive to the profile, particularly considering the phase \(\phi _{e\tau }\) [85], where \(\epsilon _{e\tau }=|\epsilon _{e\tau }|e^{i\phi _{e\tau }}\), to which DUNE will have a high sensitivity [55, 61,62,63,64], as we also see in Fig. 7.

In order to compare the results of our analysis predictions for DUNE with the constraints from other experiments, we use the results from [55]. There are differences in the nominal parameter values used for calculating the \(\chi ^2\) function and other assumptions. This is the reason why the regions in Fig. 9 do not have the same central values, but this comparison gives a good view of how DUNE can substantially improve the bounds on, for example, \(\varepsilon _{\tau \tau }-\varepsilon _{\mu \mu }\), \({\varDelta }m^2_{31}\), and the non-diagonal NSI parameters.

Fig. 9
figure 9

One-dimensional DUNE constraints compared with current constraints calculated in Ref. [55]. The left half of the figure shows constraints on the standard oscillation parameters, written in the bottom of each comparison. The five comparisons in the right half show constraints on non-standard interaction parameters

NSI can significantly impact the determination of current unknowns such as CPV and the octant of \(\theta _{23}\). Clean determination of the intrinsic CP phase at LBL experiments, such as DUNE, in the presence of NSI, is a formidable task [86]. A feasible strategy to disambiguate physics scenarios at DUNE using high-energy beams was suggested in [87]. The conclusion here is that, using a tunable beam, it is possible to disentangle scenarios with NSI. Constraints from other experiments can also solve the NSI induced degeneracy on \(\theta _{23}\).

6 CPT and Lorentz violation

Charge, parity, and time reversal symmetry (CPT) is a cornerstone of our model-building strategy. DUNE can improve the present limits on Lorentz and CPT violation by several orders of magnitude [88,89,90,91,92,93,94,95], contributing as a very important experiment to test these fundamental assumptions underlying quantum field theory.

CPT invariance is one of the predictions of major importance of local, relativistic quantum field theory. One of the predictions of CPT invariance is that particles and antiparticles have the same masses and, if unstable, the same lifetimes. To prove the CPT theorem one needs only three ingredients [88]: Lorentz invariance, hermiticity of the Hamiltonian, and locality.

Experimental bounds on CPT invariance can be derived using the neutral kaon system [96]:

$$\begin{aligned} \frac{|m(K^0) - m({\overline{K}}^0)|}{m_K} < 0.6 \times 10^{-18}. \end{aligned}$$
(9)

This result, however, should be interpreted very carefully for two reasons. First, we do not have a complete theory of CPT violation, and it is therefore arbitrary to take the kaon mass as a scale. Second, since kaons are bosons, the term entering the Lagrangian is the mass squared and not the mass itself. With this in mind, we can rewrite the previous bound as: \( |m^2(K^0) - m^2({\overline{K}}^0)| < 0.3~\text{ eV}^2 \, \). Modeling CPT violation as differences in the usual oscillation parameters between neutrinos and antineutrinos, we see here that neutrinos can test the predictions of the CPT theorem to an unprecedented extent and could, therefore, provide stronger limits than the ones regarded as the most stringent ones to date.Footnote 3

In the absence of a solid model of flavor, not to mention one of CPT violation, the spectrum of neutrinos and antineutrinos can differ both in the mass eigenstates themselves as well as in the flavor composition of each of these states. It is important to notice then that neutrino oscillation experiments can only test CPT in the mass differences and mixing angles. An overall shift between the neutrino and antineutrino spectra will be missed by oscillation experiments. Nevertheless, such a pattern can be bounded by cosmological data [97]. Unfortunately direct searches for neutrino mass (past, present, and future) involve only antineutrinos and hence cannot be used to draw any conclusion on CPT invariance on the absolute mass scale, either. Therefore, using neutrino oscillation data, we will compare the mass splittings and mixing angles of neutrinos with those of antineutrinos. Differences in the neutrino and antineutrino spectrum would imply the violation of the CPT theorem.

In Ref. [93] the authors derived the most up-to-date bounds on CPT invariance from the neutrino sector using the same data that was used in the global fit to neutrino oscillations in Ref. [98]. Of course, experiments that cannot distinguish between neutrinos and antineutrinos, such as atmospheric data from Super-Kamiokande [99], IceCube-DeepCore [100, 101] and ANTARES [102] were not included. The complete data set used, as well as the parameters to which they are sensitive, are (1) from solar neutrino data [103,104,105,106,107,108,109,110,111,112]: \(\theta _{12}\), \({\varDelta }m_{21}^2\), and \(\theta _{13}\); (2) from neutrino mode in LBL experiments K2K [113], MINOS [114, 115], T2K [116, 117], and NO\(\nu \)A [118, 119]: \(\theta _{23}\), \({\varDelta }m_{31}^2\), and \(\theta _{13}\); (3) from KamLAND reactor antineutrino data [120]: \({\overline{\theta }}_{12}\), \({\varDelta }{\overline{m}}_{21}^2\), and \({\overline{\theta }}_{13}\); (4) from short-baseline reactor antineutrino experiments Daya Bay [121], RENO [122], and Double Chooz [123]: \({\overline{\theta }}_{13}\) and \({\varDelta }{\overline{m}}_{31}^2\); and (5) from antineutrino mode in LBL experiments MINOS [114, 115] and T2K [116, 117]: \({\overline{\theta }}_{23}\), \({\varDelta }{\overline{m}}_{31}^2\), and \({\overline{\theta }}_{13}\).Footnote 4

From the analysis of all previous data samples, one can derive the most up-to-date (3\(\sigma \)) bounds on CPT violation:

$$\begin{aligned} |{\varDelta }m_{21}^2-{\varDelta }{\overline{m}}_{21}^2|&< 4.7\times 10^{-5} \, \text {eV}^2,\,\, \nonumber \\ |{\varDelta }m_{31}^2-{\varDelta }{\overline{m}}_{31}^2|&< 3.7\times 10^{-4} \, \text {eV}^2,\,\, \nonumber \\ |\sin ^2\theta _{12}-\sin ^2{\overline{\theta }}_{12}|&< 0.14,\,\, \nonumber \\ |\sin ^2\theta _{13}-\sin ^2{\overline{\theta }}_{13}|&< 0.03, \,\, \nonumber \\ |\sin ^2\theta _{23}-\sin ^2{\overline{\theta }}_{23}|&< 0.32. \end{aligned}$$
(10)

At the moment it is not possible to set any bound on \(|\delta -{\overline{\delta }}|\), since all possible values of \(\delta \) or \({\overline{\delta }}\) are allowed by data. The preferred intervals of \(\delta \) obtained in Ref. [98] can only be obtained after combining the neutrino and antineutrino data samples. The limits on \({\varDelta }({\varDelta }m_{31}^2)\) and \({\varDelta }({\varDelta }m_{21}^2)\) are already better than the one derived from the neutral kaon system and should be regarded as the best current bounds on CPT violation on the mass squared. Note that these results were derived assuming the same mass ordering for neutrinos and antineutrinos. If the ordering was different for neutrinos and antineutrinos, this would be an indication for CPT violation on its own. In the following we show how DUNE could improve this bound.

Sensitivity of the DUNE experiment to measure CPT violation in the neutrino sector is studied by analyzing neutrino and antineutrino oscillation parameters separately. We assume the neutrino oscillations being parameterized by the usual PMNS matrix \(U_{\text {PMNS}}\), with parameters \(\theta _{12},\theta _{13},\theta _{23},{\varDelta }m_{21}^2,{\varDelta }m_{31}^2, \mathrm{and}~\delta \), while the antineutrino oscillations are parameterized by a matrix \({\overline{U}}_{\text {PMNS}}\) with parameters \({\overline{\theta }}_{12},{\overline{\theta }}_{13},{\overline{\theta }}_{23},{\varDelta }{\overline{m}}_{21}^2,{\varDelta }{\overline{m}}_{31}^2, \mathrm{and}~{\overline{\delta }}\). Hence, antineutrino oscillation is described by the same probability functions as neutrinos with the neutrino parameters replaced by their antineutrino counterparts.Footnote 5 To simulate the expected neutrino data signal in DUNE, we assume the true values for neutrinos and antineutrinos to be as listed in Table 5. Then, in the statistical analysis, we vary freely all the oscillation parameters, except the solar ones, which are fixed to their best fit values throughout the simulations. Given the great precision in the determination of the reactor mixing angle by the short-baseline reactor experiments [121,122,123], in our analysis we use a prior on \({\overline{\theta }}_{13}\), but not on \(\theta _{13}\). We also consider three different values for the atmospheric angles, as indicated in Table 5. The exposure considered in the analysis corresponds to \(300~\text {kt} \, \cdot \, \text {MW} \cdot \, \text {year} \).

Table 5 Oscillation parameters used to simulate neutrino and antineutrino data for the DUNE CPT sensitivity analysis

Therefore, to test the sensitivity at DUNE we perform the simulations assuming \({\varDelta }x = |x-{\overline{x}}| = 0\), where x is any of the oscillation parameters. Then we estimate the sensitivity to \({\varDelta }x\ne 0\). To do so, we calculate two \(\chi ^2\)-grids, one for neutrinos and one for antineutrinos, varying the four parameters of interest, in this case the atmospheric oscillation parameters. After minimizing over all parameters except x and \({\overline{x}}\), we calculate

$$\begin{aligned} \chi ^2({\varDelta }x) = \chi ^2(|x-{\overline{x}}|) = \chi ^2(x)+\chi ^2({\overline{x}}), \end{aligned}$$
(11)

where we have considered all the possible combinations of \(|x-{\overline{x}}|\). The results are presented in Fig. 10, where we plot three different lines, labelled as “high”, “max” and “low.” These refer to the assumed value for the atmospheric angle: in the lower octant (low), maximal mixing (max) or in the upper octant (high). Here we can see that there is sensitivity neither to \({\varDelta }(\sin ^2\theta _{13})\), where the 3\(\sigma \) bound would be of the same order as the current measured value for \(\sin ^2{\overline{\theta }}_{13}\), nor to \({\varDelta }\delta \), where no single value of the parameter would be excluded at more than 2\(\sigma \).

On the contrary, interesting results for \({\varDelta }({\varDelta }m_{31}^2)\) and \({\varDelta }(\sin ^2\theta _{23})\) are obtained. First, we see that DUNE can put stronger bounds on the difference of the atmospheric mass splittings, namely \({\varDelta }({\varDelta }m_{31}^2) < 8.1\times 10^{-5}\), improving the current neutrino bound by one order of magnitude. For the atmospheric angle, we obtain different results depending on the true value assumed in the simulation of DUNE data. In the lower right panel of Fig. 10 we see the different behavior obtained for \(\theta _{23}\) with the values of \(\sin ^2\theta _{23}\) from Table 5, i.e., lying in the lower octant, being maximal, and lying in the upper octant. As one might expect, the sensitivity increases with \({\varDelta }\sin ^2\theta _{23}\) in the case of maximal mixing. However, if the true value lies in the lower or upper octant, a degenerate solution appears in the complementary octant.

Fig. 10
figure 10

The sensitivities of DUNE to the difference of neutrino and antineutrino parameters: \({\varDelta }\delta \), \({\varDelta }({\varDelta }m_{31}^2)\), \({\varDelta }(\sin ^2\theta _{13})\) and \({\varDelta }(\sin ^2\theta _{23})\) for the atmospheric angle in the lower octant (black line), in the upper octant (light gray line) and for maximal mixing (dark gray line)

In some types of neutrino oscillation experiments, e.g., accelerator experiments, neutrino and antineutrino data are obtained in separate experimental runs. The usual procedure followed by the experimental collaborations, as well as the global oscillation fits as for example Ref. [98], assumes CPT invariance and analyzes the full data sample in a joint way. However, if CPT is violated in nature, the outcome of the joint data analysis might give rise to what we call an “imposter” solution, i.e., one that does not correspond to the true solution of any channel.

Under the assumption of CPT conservation, the \(\chi ^2\) functions are computed according to

$$\begin{aligned} \chi ^2_{\text {total}}=\chi ^2(\nu )+\chi ^2({\overline{\nu }})\, , \end{aligned}$$
(12)

and assuming that the same parameters describe neutrino and antineutrino flavor oscillations. In contrast, in Eq. (11) we first profiled over the parameters in neutrino and antineutrino mode separately and then added the profiles. Here, we shall assume CPT to be violated in nature, but perform our analysis as if it were conserved. As an example, we assume that the true value for the atmospheric neutrino mixing is \(\sin ^2\theta _{23}=0.5\), while the antineutrino mixing angle is given by \(\sin ^2{\overline{\theta }}_{23}=0.43\). The rest of the oscillation parameters are set to the values in Table 5. Performing the statistical analysis in the CPT-conserving way, as indicated in Eq. (12), we obtain the profile of the atmospheric mixing angle presented in Fig. 11. The profiles for the individual reconstructed results (neutrino and antineutrino) are also shown in the figure for comparison. The result is a new best fit value at \(\sin ^2\theta ^\text {comb}_{23}=0.467\), disfavoring the true values for neutrino and antineutrino parameters at approximately 3\(\sigma \) and more than 5\(\sigma \), respectively.

Fig. 11
figure 11

DUNE sensitivity to the atmospheric angle for neutrinos (blue), antineutrinos (red), and to the combination of both under the assumption of CPT conservation (black)

Atmospheric neutrinos are a unique tool for studying neutrino oscillations: the oscillated flux contains all flavors of neutrinos and antineutrinos, is very sensitive to matter effects and to both \({\varDelta }m^2_{}\) parameters, and covers a wide range of L/E. In principle, all oscillation parameters could be measured, with high complementarity to measurements performed with a neutrino beam. Studying DUNE atmospheric neutrinos is also a promising approach to search for BSM effects such as Lorentz and CPT violation. The DUNE FD, with its large mass and the overburden to protect it from atmospheric muon background, is an ideal tool for these studies.

The effective field theory describing CPT violation is the Standard-Model Extension (SME) [124], where CPT violation is accompanied by Lorentz violation. This approach introduces a large set of neutrino coefficients governing corrections to standard neutrino-neutrino and antineutrino-antineutrino mixing probabilities, oscillations between neutrinos and antineutrinos, and modifications of oscillation-free propagation, all of which incorporate unconventional dependencies on the magnitudes and directions of momenta and spin. For DUNE atmospheric neutrinos, the long available baselines, the comparatively high energies accessible, and the broad range of momentum directions offer advantages that can make possible great improvements in sensitivities to certain types of Lorentz and CPT violation [90,91,92, 125,126,127,128]. To date, experimental searches for Lorentz and CPT violation with atmospheric neutrinos have been published by the IceCube and Super-Kamiokande collaborations [129,130,131]. Similar studies are possible with DUNE, and many SME coefficients can be measured that remain unconstrained to date.

Fig. 12
figure 12

Estimated sensitivity to Lorentz and CPT violation with atmospheric neutrinos in the non-minimal isotropic Standard Model Extension. The sensitivities are estimated by requiring that the Lorentz/CPT-violating effects are comparable in size to those from conventional neutrino oscillations

An example of the potential reach of studies with DUNE is shown in Fig. 12, which displays estimated sensitivities from atmospheric neutrinos in DUNE to a subset of SME coefficients controlling isotropic (rotation-invariant) violations in the Sun-centered frame [132]. The sensitivities are estimated by requiring that the Lorentz/CPT-violating effects are comparable in size to those from conventional neutrino oscillations. The eventual DUNE constraints will be determined by the ultimate precision of the experiment (which is set in part by the exposure). The gray bars in Fig. 12 show existing limits. These conservative sensitivity estimates show that DUNE can achieve first measurements (red) on some coefficients that have never previously been measured and improved measurements (green) on others, that have already been constrained in previous experiments but that can be measured with greater sensitivity with DUNE.

To illustrate an SME modification of oscillation probabilities, consider a measurement of the atmospheric neutrino and antineutrino flux as a function of energy. For definiteness, we adopt atmospheric neutrino fluxes [133], evaluated using the NRLMSISE-00 global atmospheric model [134], that result from a production event at an altitude of \(20 \, \text {km}\). Assuming conventional oscillations with standard three-flavor oscillation parameter values from the PDG [135], the fluxes at the FD are shown in Fig. 13. The sum of the \(\nu _e\) and \({{\bar{\nu }}}_e\) fluxes is shown as a function of energy as a red dashed line, while the sum of the \(\nu _\mu \) and \({{\bar{\nu }}}_\mu \) fluxes is shown as a blue dashed line. Adding an isotropic non-minimal coefficient for Lorentz violation of magnitude \({\mathring{c}}^{(6)}_{e \mu } = 1 \times 10^{-28} \, {\text {GeV}}^{-1}\) changes the fluxes from the dashed lines to the solid ones. This coefficient is many times smaller than the current experimental limit. Nonetheless, the flux spectrum is predicted to change significantly at energies over approximately \(100 \, \text {GeV}\), changing the expected number of events.

Fig. 13
figure 13

Atmospheric fluxes of neutrinos and antineutrinos as a function of energy for conventional oscillations (dashed line) and in the non-minimal isotropic Standard Model Extension (solid line)

7 Neutrino tridents at the near detector

Neutrino trident production is a weak process in which a neutrino, scattering off the Coulomb field of a heavy nucleus, generates a pair of charged leptons [136,137,138,139,140,141,142,143,144], as shown in Fig. 14.

Fig. 14
figure 14

Example diagrams for muon-neutrino-induced trident processes in the Standard Model. A second set of diagrams where the photon couples to the negatively charged leptons is not shown. Analogous diagrams exist for processes induced by different neutrino flavors and by antineutrinos. A diagram illustrating trident interactions mediated by a new \(Z'\) gauge boson, discussed in the text, is shown on the top right

Measurements of muonic neutrino tridents (\(\nu _\mu \rightarrow \nu _\mu ~\mu ^+~\mu ^-\)) were carried out at the CHARM-II [145], CCFR [146] and NuTeV [147] experiments:

$$\begin{aligned} \frac{\sigma (\nu _\mu \rightarrow \nu _\mu \mu ^+\mu ^-)_\text {exp}}{\sigma (\nu _\mu \rightarrow \nu _\mu \mu ^+\mu ^-)_\text {SM}} = {\left\{ \begin{array}{ll} 1.58 \pm 0.64 &{} \text {(CHARM-II)} \\ 0.82 \pm 0.28 &{} \text {(CCFR)} \\ 0.72 ^{+1.73}_{-0.72} &{} \text {(NuTeV)} \end{array}\right. } \end{aligned}$$

The high-intensity muon-neutrino flux at the DUNE ND will lead to a sizable production rate of trident events (see Table 6), offering excellent prospects to improve the above measurements [148,149,150]. A deviation from the event rate predicted by the SM could be an indication of new interactions mediated by the corresponding new gauge bosons [151].

Table 6 Expected number of SM \(\nu _\mu \) and \({{\bar{\nu }}}_\mu \)-induced trident events at the LArTPC of the DUNE ND per metric ton of argon and year of operation
Fig. 15
figure 15

Event kinematic distributions of signal and background considered for the selection of muonic trident interactions in the ND LArTPC: number of tracks (top left), angle between the two main tracks (top right), length of the shortest track (bottom left), and the difference in length between the two main tracks (bottom right). The dashed, black vertical lines indicate the optimal cut values used in the analysis

The main challenge in obtaining a precise measurement of the muonic trident cross section will be the copious backgrounds, mainly consisting of CC single-pion production events, \(\nu _\mu N \rightarrow \mu \pi N^\prime \), as muon and pion tracks can be easily confused in LArTPC detectors. The discrimination power of the DUNE ND LArTPC was evaluated using large simulated data sets of signal and background. Each simulated event represents a different neutrino-argon interaction in the active volume of the detector. Signal events were generated using a standalone code [148] that simulates trident production of muons and electrons through the scattering of \(\nu _{\mu }\) and \(\nu _e\) on argon nuclei. The generator considers both the coherent scattering on the full nucleus (the dominant contribution) and the incoherent scattering on individual nucleons. Background events, consisting of several SM neutrino interactions, were generated using GENIE. Roughly \(38\%\) of the generated events have a charged pion in the final state, leading to two charged tracks with muon-like energy deposition pattern (\(\mathrm {d}E/\mathrm {d}x\)), as in the trident signal. All final-state particles produced in the interactions were propagated through the detector geometry using the Geant4-based simulation of the DUNE ND. Charge collection and readout were not simulated, and possible inefficiencies due to mis-reconstruction effects or event pile-up were disregarded for simplicity.

Figure 15 shows the distribution (area normalized) for signal and background of the different kinematic variables used in our analysis for the discrimination between signal and background. As expected, background events tend to contain a higher number of tracks than the signal. The other distributions also show a clear discriminating power: the angle between the two tracks is typically much smaller in the signal than in the background. Moreover, the signal tracks (two muons) tend to be longer than tracks in the background (mainly one muon plus one pion).

The sensitivity of neutrino tridents to heavy new physics (i.e., heavy compared to the momentum transfer in the process) can be parameterized in a model-independent way using a modification of the effective four-fermion interaction Hamiltonian. Focusing on the case of muon neutrinos interacting with muons, the vector and axial-vector couplings can be written as

$$\begin{aligned} g_{\mu \mu \mu \mu }^V= & {} 1 + 4 \sin ^2\theta _W + {\varDelta }g_{\mu \mu \mu \mu }^V \quad \mathrm {and} \nonumber \\ \quad g_{\mu \mu \mu \mu }^A= & {} -1 + {\varDelta }g_{\mu \mu \mu \mu }^A ~, \end{aligned}$$
(13)

where \({\varDelta }g_{\mu \mu \mu \mu }^V\) and \({\varDelta }g_{\mu \mu \mu \mu }^A\) represent possible new physics contributions. Couplings involving other combinations of lepton flavors can be modified analogously. Note, however, that for interactions that involve electrons, very strong constraints can be derived from LEP bounds on electron contact interactions [152]. The modified interactions of the muon-neutrinos with muons alter the cross section of the \(\nu _\mu N \rightarrow \nu _\mu \mu ^+\mu ^- N\) trident process. In Fig. 16 we show the regions in the \({\varDelta }g^V_{\mu \mu \mu \mu }\) vs. \({\varDelta }g^A_{\mu \mu \mu \mu }\) plane that are excluded by the existing CCFR measurement \(\sigma _\text {CCFR} / \sigma _\text {CCFR}^\text {SM} = 0.82 \pm 0.28\) [146] at the 95% CL in gray. A measurement of the \(\nu _\mu N \rightarrow \nu _\mu \mu ^+\mu ^- N\) cross section with \(40\%\) uncertainty (obtained after running for \(\sim 6\)  years in neutrino mode or, equivalently, 3 years in neutrino mode and 3 years in antineutrino mode) at the DUNE ND could cover the blue hashed regions (95% CL). These numbers show that a measurement of the SM di-muon trident production at the 40% level could be possible. Our baseline analysis does not extend the sensitivity into parameter space that is unconstrained by the CCFR measurement. However, it is likely that the use of a magnetized spectrometer, as it is being considered for the DUNE ND, able to identify the charge signal of the trident final state, along with a more sophisticated event selection (e.g., deep-learning-based), will significantly improve separation between neutrino trident interactions and backgrounds. Therefore, we also present the region (blue dashed line) that could be probed by a 25% measurement of the neutrino trident cross section at DUNE, which would extend the coverage of new physics parameter space substantially.

Fig. 16
figure 16

95% CL sensitivity of a 40% (blue hashed regions) and a 25% (dashed contours) uncertainty measurement of the \(\nu _\mu N \rightarrow \nu _\mu \mu ^+\mu ^- N\) cross section at the DUNE near detector to modifications of the vector and axial-vector couplings of muon-neutrinos to muons. The gray regions are excluded at 95% CL by existing measurements of the cross section by the CCFR Collaboration. The intersection of the thin black lines indicates the SM point. A 40% precision measurement could be possible with 6 years of data taking in neutrino mode

We consider a class of models that modify the trident cross section through the presence of an additional neutral gauge boson, \(Z'\), that couples to neutrinos and charged leptons. A consistent way of introducing such a \(Z'\) is to gauge an anomaly-free global symmetry of the SM. Of particular interest is the \(Z'\) that is based on gauging the difference of muon-number and tau-number, \(L_\mu - L_\tau \) [153, 154]. Such a \(Z'\) is relatively weakly constrained and can for example address the longstanding discrepancy between SM prediction and measurement of the anomalous magnetic moment of the muon, \((g-2)_\mu \) [155, 156]. The \(L_\mu - L_\tau \, Z'\) has also been used in models to explain B physics anomalies [157] and as a portal to dark matter (DM) [158, 159]. The \(\nu _\mu N \rightarrow \nu _\mu \mu ^+\mu ^- N\) trident process has been identified as an important probe of gauged \(L_\mu - L_\tau \) models over a broad range of \(Z^\prime \) masses [151, 157].

Fig. 17
figure 17

Existing constraints and projected DUNE sensitivity in the \(L_\mu - L_\tau \) parameter space. Shown in green is the region where the \((g-2)_\mu \) anomaly can be explained at the \(2\sigma \) level. The parameter regions already excluded by existing constraints are shaded in gray and correspond to a CMS search for \(pp \rightarrow \mu ^+\mu ^- Z' \rightarrow \mu ^+\mu ^-\mu ^+\mu ^-\) [160] (“LHC”), a BaBar search for \(e^+e^- \rightarrow \mu ^+\mu ^- Z' \rightarrow \mu ^+\mu ^-\mu ^+\mu ^-\) [161] (“BaBar”), a previous measurement of the trident cross section [146, 151] (“CCFR”), a measurement of the scattering rate of solar neutrinos on electrons [162,163,164] (“Borexino”), and bounds from Big Bang Nucleosynthesis [165, 166] (“BBN”). The DUNE sensitivity shown by the solid blue line assumes 6 years of data running in neutrino mode, leading to a measurement of the trident cross section with 40% precision

In Fig. 17 we show the existing CCFR constraint on the model parameter space in the \(m_{Z'}\) vs. \(g'\) plane, where \(g'\) is the \(L_\mu - L_\tau \) gauge coupling, and compare it to the region of parameter space where the anomaly in \((g-2)_\mu = 2 a_\mu \) can be explained. The green region shows the \(1\sigma \) and \(2\sigma \) preferred parameter space corresponding to a shift \({\varDelta }a_\mu = a_\mu ^\text {exp}-a_\mu ^\text {SM} = (2.71 \pm 0.73) \times 10^{-9}\) [167]. In addition, constraints from LHC searches for the \(Z'\) in the \(pp \rightarrow \mu ^+\mu ^- Z' \rightarrow \mu ^+\mu ^-\mu ^+\mu ^-\) process [160] (see also [151]) and direct searches for the \(Z'\) at BaBar using the \(e^+e^- \rightarrow \mu ^+\mu ^- Z' \rightarrow \mu ^+\mu ^-\mu ^+\mu ^-\) process [161] are shown. A Borexino bound on non-standard contributions to neutrino-electron scattering [162,163,164] has also been used to constrain the \(L_\mu - L_\tau \) gauge boson [166, 168, 169]. Our reproduction of the Borexino constraint is shown in Fig. 17. For very light \(Z'\) masses of O(few MeV) and below, strong constraints from measurements of the effective number of relativistic degrees of freedom during Big Bang Nucleosynthesis (BBN) apply [165, 166]. Taking into account all relevant constraints, parameter space to explain \((g-2)_\mu \) is left below the di-muon threshold \(m_{Z'} \lesssim 210~\text {MeV}\). The DUNE sensitivity shown by the solid blue line assumes a measurement of the trident cross section with \(40\%\) precision.

8 Dark matter probes

Dark matter is a crucial ingredient to understand the cosmological history of the universe, and the most up-to-date measurements suggests the existence of DM with a density parameter (\({\varOmega }_{c}\)) of 0.264 [170]. In light of this situation, a tremendous amount of experimental effort has gone into the search for DM-induced signatures, for example, DM direct and indirect detections and collider searches. However, no “smoking-gun” signals have been discovered thus far while more parameter space in relevant DM models is simply ruled out. It is noteworthy that most conventional DM search strategies are designed to be sensitive to signals from the weakly-interacting massive particle (WIMP), one of the well-motivated DM candidates, whose mass range is from a few GeV to tens of TeV. The non-observation of DM via non-gravitational interactions actually motivates unconventional or alternative DM search schemes. One such possibility is a search for experimental signatures induced by boosted, hence relativistic, DM for which a mass range smaller than that of the weak scale is often motivated.

One of the possible ways to produce and then detect relativistic DM particles can be through accelerator experiments, for example, neutrino beam experiments [3, 171,172,173,174]. Due to highly intensified beam sources, large signal statistics is usually expected so that this sort of search strategy can allow for significant sensitivity to DM-induced signals despite the feeble interaction of DM with SM particles. DUNE will perform a search for the relativistic scattering of light-mass dark matter (LDM), whose lowest mass particle is denoted as \(\chi \) throughout this section, at the ND, as it is close enough to the beam source to sample a substantial level of DM flux, assuming that DM is produced.

Alternatively, it is possible that boosted dark matter (BDM) particles are created in the universe under non-minimal dark-sector scenarios [175, 176], and can reach terrestrial detectors. For example, one can imagine a two-component DM scenario in which a lighter component (\(\chi \)) is usually a subdominant relic with direct coupling to SM particles, while the heavier (denoted as \(\psi \) throughout this section) is the cosmological DM that pair-annihilates directly to a lighter DM pair, not to SM particles. Other mechanisms such as semi-annihilation in which a DM particle pair-annihilates to a (lighter) DM particle and a dark sector particle that may decay away are also possible [177,178,179]. In typical cases, the BDM flux is not large and thus large-volume neutrino detectors are desirable to overcome the challenge in statistics (for an exception, see [180,181,182,183]).

Indeed, a (full-fledged) DUNE FD with a fiducial mass of \(40 \, \text {kt}\) and quality detector performance is expected to possess competitive sensitivity to BDM signals from various sources in the current universe such as the galactic halo [175, 181, 184,185,186,187,188], the sun [178, 179, 184, 187, 189], and dwarf spheroidal galaxies [188]. Furthermore, the ProtoDUNE detectors have taken data, and we anticipate preliminary studies with their cosmic data. Interactions of BDM with electrons [175] and with hadrons (protons) [179], were investigated for Cherenkov detectors, such as Super-Kamiokande, which recently published a dedicated search for BDM in the electron channel [190]. However, in such detectors the BDM signal rate is shown to often be significantly attenuated due to Cherenkov threshold, in particular for hadronic channels. LAr detectors, such as DUNE’s, have the potential to greatly improve the sensitivity for BDM compared to Cherenkov detectors. This is due to improved particle identification techniques, as well as a significantly lower energy threshold for proton detection. Earlier studies have shown an improvement with DUNE for BDM-electron interaction [188].

We consider several benchmark DM models. These describe only couplings of dark-sector states including LDM particles. We consider two example models: (i) a vector portal-type scenario where a (massive) dark-sector photon V mixes with the SM photon and (ii) a leptophobic \(Z'\) scenario. DM and other dark-sector particles are assumed to be fermionic for convenience.

Benchmark Model (i) The relevant interaction Lagrangian is given by [185]

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\mathrm{int}}&\supset -\frac{\epsilon }{2}V_{\mu \nu }F^{\mu \nu } +g_D {\bar{\chi }}\gamma ^\mu \chi V_\mu \\&\quad +g'_D {\bar{\chi }}'\gamma ^\mu \chi V_\mu +h.c. , \end{aligned} \end{aligned}$$
(14)

where \(V^{\mu \nu }\) and \(F^{\mu \nu }\) are the field strength tensors for the dark-sector photon and the SM photon, respectively. Here we have introduced the kinetic mixing parameter \(\epsilon \), while \(g_D\) and \(g'_D\) parameterize the interaction strengths for flavor-conserving (second operator) and flavor-changing (third operator) couplings, respectively. Here \(\chi \) and \(\chi '\) denote a dark matter particle and a heavier, unstable dark-sector state, respectively (i.e., \(M_{\chi '}>M_{\chi }\)), and the third term allows (boosted) \(\chi \) transition to \(\chi '\) after a scattering (i.e., an “inelastic” scattering process).

This model introduces six new free parameters that may be varied for our sensitivity analysis: dark photon mass \(M_V\), DM mass \(M_{\chi }\), heavier dark-sector state mass \(M_{\chi '}\), kinetic mixing parameter \(\epsilon \), dark-sector couplings \(g_D\) and \(g'_D\). We shall perform our analyses with some of the parameters fixed to certain values for illustration.

Benchmark Model (ii) This model employs a leptophobic \(Z^\prime \) mediator for interactions with the nucleons. The interaction Lagrangian for this model is [179]

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_{\mathrm{int}}&\supset - g_{\mathrm{Z^\prime }} \sum _f Z^\prime _\mu {\bar{q}}_f \gamma ^\mu \gamma ^5 q_f - g_{\mathrm{Z^\prime }} Z^\prime _\mu {\bar{\chi }} \gamma ^\mu \gamma ^5 \chi \\&\quad - Q_\psi g_{\mathrm{Z^\prime }} Z^\prime _\mu {\bar{\psi }} \gamma ^\mu \gamma ^5 \psi . \end{aligned} \end{aligned}$$
(15)

Here, all couplings are taken to be axial. f denotes the quark flavors in the SM sector. The dark matter states are denoted by \(\chi \) and \(\psi \) with \(M_\chi < M_\psi \). The coupling \(g_{\mathrm{Z^\prime }}\) and the masses of the dark matter states are free parameters. The DM flux abundance parameter, \(Q_\psi \) is taken to be less than 1 and determines the abundance of dark matter in the universe. The hadronic interaction model study presented here is complementary to and has different phenomenology compared to others such as Benchmark Model (i).

Table 7 A summary of the three different studies in this section

We summarize key information for the three different studies in this section in Table 7. The \(e^-\) (p) outside (inside) the parentheses in the third column imply the electron (proton) scattering channel. N in the last column denotes a nucleon, while X stands for particle(s) created via the \(\chi -N\) scattering process.

8.1 Search for low-mass dark matter at the near detector

Here, we focus on Benchmark Model (i) from Eq. (14), specifically where only one DM particle \(\chi \) is relevant. We also define the dark fine structure constant \(\alpha _D \equiv g_D^2/(4\pi )\). We assume that \(\chi \) is a fermionic thermal relic – in this case, the DM/dark photon masses and couplings will provide a target for which the relic abundance matches the observed abundance in the universe. Here, the largest flux of dark photons V and DM to reach the DUNE ND will come from the decays of light pseudoscalar mesons (specifically \(\pi ^0\) and \(\eta \) mesons) that are produced in the DUNE target, as well as proton bremsstrahlung processes \(p + p \rightarrow p + p + V\). For the entirety of this analysis, we will fix \(\alpha _D = 0.5\) and assume that the DM mass \(M_{\chi }\) is lighter than half the mass of a pseudoscalar meson \({\mathfrak {m}}\) that is produced in the DUNE target. In this scenario, \(\chi \) is produced via two decays, those of on-shell V and those of off-shell V. This production is depicted in Fig. 18.

The flux of DM produced via meson decays – via on-shell V – may be estimated byFootnote 6

$$\begin{aligned} N_\chi= & {} 2 N_{\mathrm {POT}} c_{\mathfrak {m}} \{ {\mathrm {Br}}({\mathfrak {m}}\rightarrow \gamma \gamma ) \times 2 \varepsilon ^2 \left( 1 - \frac{M_{V}^2}{m_{\mathrm {m}}^2}\right) ^3 \nonumber \\&\times {\mathrm {Br}}(V \rightarrow \chi {\bar{\chi }}) \} g(M_\chi , M_{V}), \end{aligned}$$
(16)

where \(N_{\mathrm {POT}}\) is the number of protons on target delivered by the beam, \(c_{\mathfrak {m}}\) is the average number of meson \({\mathfrak {m}}\) produced per POT, the term in braces is the relative branching fraction of \({\mathfrak {m}} \rightarrow \gamma V\) relative to \(\gamma \gamma \), and g(xy) characterizes the geometrical acceptance fraction of DM reaching the DUNE ND. g(xy) is determined given model parameters using Monte Carlo techniques. For the range of dark photon and DM masses in which DUNE will set a competitive limit, the DM flux due to meson decays will dominate over the flux due to proton bremsstrahlung. Considering DM masses in the \(\sim \)1–300 MeV range, this will require production via the \(\pi ^0\) and \(\eta \) mesons. Our simulations using Pythia determine that \(c_{\pi ^0} \approx 4.5\) and \(c_\eta \approx 0.5\).

Fig. 18
figure 18

Production of fermionic DM via two-body pseudoscalar meson decay \({\mathfrak {m}} \rightarrow \gamma V\), when \(M_{V} < m_{\mathfrak {m}}\) (top) or via three-body decay \({\mathfrak {m}} \rightarrow \gamma \chi {\overline{\chi }}\) (center) and DM-electron elastic scattering (bottom)

If the DM reaches the near detector, it may scatter elastically off nucleons or electrons in the detector, via a t-channel dark photon. Due to its smaller backgrounds, we focus on scattering off electrons, depicted in the bottom panel of Fig. 18. The differential cross section of this scattering, as a function of the recoil energy of the electron \(E_e\), is

$$\begin{aligned} \frac{d\sigma _{{\chi }e}}{dE_{e}}= & {} 4\pi \epsilon ^{2}\alpha _D\alpha _{EM} \nonumber \\&\times \frac{2m_{e}E_{\chi }^{2} - (2m_{e}E_{\chi } + M_{\chi }^{2})(E_e-m_{e})}{(E_e^{2}-M_{\chi }^{2})(M_{V}^{2} +2m_{e}E_{e}-2m_{e}^{2})^{2}}, \end{aligned}$$
(17)

where \(E_{\chi }\) is the incoming DM \(\chi \) energy. The signal is an event with only one recoil electron in the final state. We can exploit the difference between the scattering angle and the energy of the electron to distinguish between signal and the background from neutrino-electron scattering (discussed in the following) events.

The background to the process shown in the bottom panel of Fig. 18 consists of any processes involving an electron recoil. As the ND is located near the surface, background events, in general, can be induced by cosmic rays as well as by neutrinos generated from the beam. Since the majority of cosmic-induced events, however, will be vetoed by triggers and timing information, the dominant background will be from neutrinos coming in the DUNE beam.

The two neutrino-related backgrounds are \(\nu _\mu -e^-\) scattering, which looks nearly identical to the signal, and \(\nu _e\) CCQE scattering, which does not. The latter has a much larger rate (\(\sim \) 10 times higher) than the former, however, we expect that using the kinematical variable \(E_e \theta _e^2\) of the final state, where \(\theta _e\) is the direction of the outgoing electron relative to the beam direction, will enable us to exploit the differences in the scattering angle of the electron from the DM interactions to reduce a substantial fraction of the \(\nu _e\) CCQE background [192].

While spectral information regarding \(E_e\) could allow a search to distinguish between \(\chi e\) and \(\nu _\mu e\) scattering, we expect that uncertainties in the \(\nu _\mu \) flux (both in terms of overall normalization and shape as a function of neutrino energy) will make such an analysis very complicated. For this reason, we include a normalization uncertainty of \(10\%\) on the expected background rate and perform a counting analysis. Studies are ongoing to determine how such an analysis may be improved.

For this analysis we have assumed 3.5 years of data collection each in neutrino and antineutrino modes, analyzing events that occur within the fiducial volume of the DUNE near detector. We compare results assuming either all data is collected with the ND on-axis, or data collection is divided equally among all off-axis positions, 0.7 year at each position i, between 0 and 24 m transverse to the beam direction (in steps of 6 meters). We assume three sources of uncertainty: statistical, correlated systematic, and an uncorrelated systematic in each bin. For a correlated systematic uncertainty, we include a nuisance parameter A that modifies the number of neutrino-related background events in all bins – an overall normalization uncertainty across all off-axis locations.

We further include an additional term in our test statistic for A, a Gaussian probability with width \(\sigma _A = 10\%\). We also include an uncorrelated uncertainty in each bin, which we assume to be much narrower than \(\sigma _A\). We assume this uncertainty to be parameterized by a Gaussian with width \(\sigma _{f_i} = 1\%\). After marginalizing over the corresponding uncorrelated nuisance parameters, the test statistic reads

$$\begin{aligned} -2{\varDelta }{\mathcal {L}} =&\sum _i \frac{r_i^m\left( \left( \frac{\varepsilon }{\varepsilon _0}\right) ^4 N_i^\chi + (A-1)N_i^\nu \right) ^2}{A\left( N_i^\nu + (\sigma _{f_i} N_i^\nu )^2 \right) } \nonumber \\&+ \frac{\left( A-1\right) ^2}{\sigma _A^2}. \end{aligned}$$
(18)

In Eq. (18), \(N_i^\chi \) is the number of DM scattering events, calculated assuming \(\varepsilon \) is equal to some reference value \(\varepsilon _0 \ll 1\). \(N_i^\nu \) is the number of \(\nu _\mu e^-\) scattering events expected in detector position i, and \(r_i^m\) is the number of years of data collection in detector position i during beam mode m (neutrino or antineutrino mode). If data are only collected on-axis, then this test statistic will be dominated by the systematic uncertainty associated with \(\sigma _A\). If on- and off-axis measurements are combined, then the resulting sensitivity will improve significantly.

We present results in terms of the DM or dark photon mass and the parameter Y, where

$$\begin{aligned} Y \equiv \varepsilon ^2 \alpha _D \left( \frac{M_\chi }{M_V}\right) ^4. \end{aligned}$$
(19)

Assuming \(M_V \gg M_\chi \), this parameter determines the relic abundance of DM in the universe today, and sets a theoretical goal in terms of sensitivity reach. We present the 90% CL sensitivity reach of the DUNE ND in Fig. 19. We assume \(\alpha _D = 0.5\) in our simulations and we display the results fixing \(M_V = 3M_\chi \) (left panel) and \(M_\chi = 20 \, \text {MeV}\) (right panel). We also compare the sensitivity reach of this analysis with other existing experiments, shown as grey shaded regions. We further show for comparison the sensitivity curve expected for a proposed dedicated experiment to search for LDM, LDMX-Phase I [193] (solid blue).

Fig. 19
figure 19

Expected DUNE On-axis (solid red) and PRISM (dashed red) sensitivity using \(\chi e^- \rightarrow \chi e^-\) scattering. We assume \(\alpha _D = 0.5\) in both panels, and \(M_V = 3M_\chi \) (\(M_\chi = 20 \, \text {MeV}\)) in the left (right) panel, respectively. Existing constraints are shown in grey, and the relic density target is shown as a black line. We also show for comparison the sensitivity curve expected for LDMX-Phase I (solid blue) [193]

From our estimates, we see that DUNE can significantly improve the constraints from LSND [194] and the MiniBooNE-DM search [195], as well as BaBar [196] if \(M_V \lesssim 200 \, \text {MeV}\). We also show limits in the right panel from beam-dump experiments (where the dark photon is assumed to decay visibly if \(M_V < 2 M_\chi \)) [197,198,199,200,201,202], as well as the lower limits obtained from matching the thermal relic abundance of \(\chi \) with the observed one (black).

The features in the sensitivity curve in the right panel can be understood by looking at the DM production mechanism. For a fixed \(\chi \) mass, as \(M_V\) grows, the DM production goes from off-shell to on-shell and back to off-shell. The first transition explains the strong feature near \(M_V=2M_\chi = 40~\text {MeV}\), while the second is the source for the slight kink around \(M_V=m_{\pi ^0}\) (which appears also in the left panel).

8.2 Inelastic boosted dark matter search at the DUNE FD

We consider an annihilating two-component DM scenario [176] in this study. The heavier DM (denoted \({\varPsi }\)) plays a role of cosmological DM and pair-annihilates to a pair of lighter DM particles (denoted \(\chi \)) in the universe today. The expected flux near the earth is given by [175, 181, 187]

$$\begin{aligned} {\mathcal {F}}_1&= 1.6 \times 10^{-6} \mathrm{cm}^{-2}\mathrm{s}^{-1}\times \left( \frac{\langle \sigma v\rangle _{{\varPsi }\rightarrow \chi }}{5\times 10^{-26}\mathrm{cm}^3\mathrm{s}^{-1}}\right) \nonumber \\&\quad \times \left( \frac{10\, \mathrm{GeV}}{M_{{\varPsi }}}\right) ^2, \end{aligned}$$
(20)

where \(m_{{\varPsi }}\) is the mass of \({\varPsi }\) and \(\langle \sigma v\rangle _{{\varPsi }\rightarrow \chi }\) stands for the velocity-averaged annihilation cross section of \({\varPsi }{\bar{{\varPsi }}} \rightarrow \chi {\bar{\chi }}\) in the current universe. To evaluate the reference value shown as the first prefactor, we take \(M_{{\varPsi }} = 10 \, \text {GeV}\) and \(\langle \sigma v\rangle _{{\varPsi }\rightarrow \chi }=5\times 10^{-26}~\mathrm{cm}^3\mathrm{s}^{-1}\), the latter of which is consistent with the current observation of DM relic density assuming \({\varPsi }\) and its anti-particle \({\bar{{\varPsi }}}\) are distinguishable. To integrate all relevant contributions over the entire galaxy, we assume the Navarro–Frenk–White (NFW) DM halo profile [203, 204]. In this section we assume the BDM flux with a \(M_{{\varPsi }}\) dependence given by Eq. (20) for the phenomenological analysis.

The BDM that is created, e.g., at the galactic center, reaches the DUNE FD detectors and scatters off either electrons or protons energetically. In this study, we focus on electron scattering signatures for illustration, under Benchmark Model (i) defined in Eq. (14). The overall process is summarized as follows:

$$\begin{aligned}&\chi + e^-~(\mathrm{or}~p) \nonumber \\&\quad \rightarrow e^-~(\mathrm{or}~p) + \chi ' (\rightarrow \chi + V^{(*)} \rightarrow \chi + e^+ +e^-), \end{aligned}$$
(21)

where \(\chi '\) is a dark-sector unstable particle that is heavier than \(\chi \) as described earlier. A diagrammatic description is shown in Fig. 20 where particles visible by the detector are circled in blue. In the final state of the e-scattering case, there exist three visible particles that usually leave sizable (e-like) tracks in the detectors. On the other hand, for the p-scattering case we can replace \(e^-\) in the left-hand side and the first \(e^-\) in the right-hand side of the above process by p. In the basic model, Eq. (14), and given the source of BDM at the galactic center, the resulting signature accompanies a quasi-elastic proton recoil [205] together with a pair of \(e^+e^-\) tracks.

Fig. 20
figure 20

The inelastic BDM signal under consideration

As we have identified a possible inelastic BDM (iBDM) signature, we are now in a position to discuss potential SM background events. For the DUNE detector modules located \(\sim 1480\) m deep underground, the cosmic-induced backgrounds are not an issue except the background induced by atmospheric neutrinos. The most plausible scenario for background production is that an atmospheric neutrino event involves the creation of multiple pions that subsequently decay to electrons, positrons, photons, and neutrinos. Relevant channels are the resonance production and/or deep inelastic scattering (DIS) by the CC \(\nu _e\) or \({{\bar{\nu }}}_e\) scattering with a nucleon in the LAr target. Summing up all the resonance production and DIS events that are not only induced by \(\nu _e\) or \({\bar{\nu }}_e\) but relevant to production of a few pions, we find that the total number of multi-pion production events is at most \(\sim 20 \, (\hbox {kt} \, \cdot \, \hbox {year})^{-1}\) [206], based on the neutrino flux calculated in Ref. [133] and the cross section in Ref. [207]. In addition, the charged pions often leave long enough tracks inside the detector so that the probability of misidentifying the \(e^\pm \) from the decays of \(\pi ^\pm \) with the iBDM signal events would be very small. Some quasi-elastic scattering events by atmospheric neutrinos may involve a detectable proton recoil together with a single e-like track, which might behave like backgrounds in the proton scattering channel. However, this class of events can be rejected by requiring two separated e-like tracks. Hence, we conclude that it is fairly reasonable to assume that almost no background events exist. See also Ref. [206] for a more systematic background consideration for the iBDM signals.

We finally present the expected experimental sensitivities of DUNE, in the searches for iBDM. We closely follow the strategies illustrated in Refs. [181, 205] to represent phenomenological interpretations. In displaying the results, we separate the signal categories into

  • Scenario 1: \(M_V > 2 M_{\chi }\), experimental limits for \(V \rightarrow \) invisible applied.

  • Scenario 2: \(M_V \le 2 M_{\chi }\), experimental limits for \(V \rightarrow e^+ e^-\) applied.

We develop an event simulation code using the ROOT package with the matrix elements for the \(\chi \) scattering and the \(\chi '\) decays implemented. Once an event is generated, we require that all the final state particles should pass the (kinetic) energy threshold (30 MeV for electrons and protons) and their angular separation from the other particles should be greater than the angular resolution (\(1^\circ \) for electrons and \(5^\circ \) for protons) [206].

Fig. 21
figure 21

The experimental sensitivities in terms of reference model parameters \(M_V - \epsilon \) for \(M_{{\varPsi }} = 0.4 \, \text {GeV}\), \(M_{\chi } = 5 \, \text {MeV}\), and \(\delta M = M_{\chi '} - M_{\chi } = 10 \, \text {MeV}\) (top-left panel) and \(M_{{\varPsi }} = 2 \, \text {GeV}\), \(M_{\chi '} = 50 \, \text {MeV}\), and \(\delta M = 10 \, \text {MeV}\) (top-right panel). The left panels are for Scenario 1 and the right ones are for Scenario 2. The bottom panels compare different reference points in the p-scattering channel. See the text for the details

We first show the results for Scenario 1 in the left panels of Fig. 21, taking a parameter set, \(M_{{\varPsi }} = 0.4~\text {GeV}\), \(M_{\chi } = 5~\text {MeV}\), \(\delta M \equiv M_{\chi '}-M_\chi = 10~\text {MeV}\) with \(g'_D=1\). The brown-shaded region shows the latest limits set by various experiments such as the fixed-target experiment NA64 [208] at the CERN SPS and the B-factory experiment BaBar [209]. Note that some of the limits are from ongoing experiments such as NA64 which will collect more data in the next years and improve their sensitivity reaches. The blue solid and the green solid lines describe the experimental sensitivityFootnote 7 of DUNE FD to the e-scattering and p-scattering signals, respectively, under a zero background assumption. The associated exposure is \(40~\text {kt} \, \cdot \, \text {year} \), i.e., a total fiducial volume of \(40 \, \text {kt}\) times one year of running time.

For Scenario 2 (the right panels of Fig. 21), we choose a different reference parameter set: \(M_{{\varPsi }} = 2~\text {GeV}\), \(M_{\chi } = 50~\text {MeV}\), \(\delta M = 10~\text {MeV}\) with \(g'_D=1\). The current limits (brown shaded regions), from various fixed target experiments, B-factory experiments, and astrophysical observations, are taken from Refs. [210, 211].

In both scenarios, the proton scattering channel enables us to explore different regions of parameter space as it allows heavier \(\chi '\) to be accessible which would be kinematically forbidden to access in the electron scattering channel. Inspired by this potential of the proton scattering channel, we study other reference parameters and compare them with the original ones in the top panels of Fig. 21, and show the results in the bottom panels. We see that different parameter choices in the proton scattering channel allow us to cover a wider or different range of parameter space.

We next discuss model-independent experimental sensitivities. The experimental sensitivities are determined by the number of signal events excluded at 90% CL in the absence of an observed signal. The expected number of signal events, \(N_{\mathrm{sig}}\), is given by

$$\begin{aligned} N_{\mathrm{sig}} = \sigma _\epsilon {\mathcal {F}} A(\ell _{\mathrm{lab}}) t_{\mathrm{exp}} N_T, \end{aligned}$$
(22)

where \(N_T\) is the number of target particles T, \(\sigma _\epsilon \) is the cross section of the primary scattering \(\chi T \rightarrow \chi ' T\), \({\mathcal {F}}\) is the flux of \(\chi \), \(t_{\mathrm{exp}}\) is the exposure time, and \(A(\ell _{\mathrm{lab}})\) is the acceptance that is defined as 1 if the event occurs within the fiducial volume and 0 otherwise. Here we determine the acceptance for an iBDM signal by the distance between the primary and secondary vertices in the laboratory frame, \(\ell _{\mathrm{lab}}\), so \(A(\ell _{\mathrm{lab}}) = 1\) when both the primary and secondary events occur inside the fiducial volume. (Given this definition, obviously, \(A(\ell _{\mathrm{lab}}) = 1\) for elastic BDM.) Our notation \(\sigma _\epsilon \) includes additional realistic effects from cuts, threshold energy, and the detector response, hence it can be understood as the fiducial cross section.

The 90% CL exclusion limit, \(N_s^{90}\), can be obtained with a modified frequentist construction [212, 213]. We follow the methods in Refs. [214,215,216] in which the Poisson likelihood is assumed. An experiment becomes sensitive to the signal model independently if \(N_{\mathrm{sig}} \ge N_s^{90}\). Plugging Eq. (22) here, we find the experimental sensitivity expressed by

$$\begin{aligned} \sigma _\epsilon {\mathcal {F}} \ge \frac{N_s^{90}}{A(\ell _{\mathrm{lab}}) t_{\mathrm{exp}} N_T}. \end{aligned}$$
(23)

Since \(\ell _{\mathrm{lab}}\) differs event-by-event, we take the maximally possible value of laboratory-frame mean decay length, i.e., \({\bar{\ell }}_{\mathrm{lab}}^{\mathrm{max}} \equiv \gamma _{\chi '}^{\max } {\bar{\ell }}_{\mathrm{rest}}\) where \(\gamma _{\chi '}^{\max }\) is the maximum boost factor of \(\chi '\) and \({\bar{\ell }}_{\mathrm{rest}}\) is the rest-frame mean decay length. We emphasize that this is a rather conservative approach, because the acceptance A is inversely proportional to \(\ell _{\mathrm{lab}}\). We then show the experimental sensitivity of any kind of experiment for a given background expectation, exposure time, and number of targets, in the plane of \({\bar{\ell }}_{\mathrm{lab}}^{\mathrm{max}} - \sigma _\epsilon \cdot {\mathcal {F}}\). The top panel of Fig. 22 demonstrates the expected model-independent sensitivities at the DUNE experiment. The green (blue) line is for the DUNE FD with a background-free assumption and 20 (40) \(\text {kt} \, \cdot \, \text {year}\) exposure.

Fig. 22
figure 22

Top: model-independent experimental sensitivities of iBDM search in \({\bar{\ell }}_{\mathrm{lab}}^{\mathrm{max}} - \sigma _\epsilon \cdot {\mathcal {F}}\) plane. The reference experiments are DUNE \({20} \, \text {kt}\) (green), and DUNE \({40} \, \text {kt}\) (blue) with zero-background assumption for 1-year time exposure. Bottom: Experimental sensitivities of iBDM search in \(M_{{\varPsi }} - \sigma _\epsilon \) plane. The sensitivities for \({\bar{\ell }}_{\mathrm{lab}}^{\mathrm{max}} = 0\) and 100 m are shown as solid and dashed lines for each reference experiment in the top panel

The bottom panel of Fig. 22 reports model-dependent sensitivities for \({\bar{\ell }}_{\mathrm{lab}}^{\mathrm{max}} = 0\) m and 100 m corresponding to the experiments in the top panel. Note that this method of presentation is reminiscent of the widely known scheme for showing the experimental reaches in various DM direct detection experiments, i.e., \(M_{\mathrm{DM}} - \sigma _{\mathrm{DM - target}}\) where \(M_{\mathrm{DM}}\) is the mass of DM and \(\sigma _{\mathrm{DM - target}}\) is the cross section between the DM and target. For the case of non-relativistic DM scattering in the direct-detection experiments, \(M_{\mathrm{DM}}\) determines the kinetic energy scale of the incoming DM, just like \(M_{{\varPsi }}\) sets out the incoming energy of boosted \(\chi \) in the iBDM search.

8.3 Elastic boosted dark matter from the sun

In this section, we focus on Benchmark Model (ii) described by Eq. (15). This study uses DUNE’s full FD event generation and detector simulation. We focus on BDM flux sourced by DM annihilation in the core of the sun. DM particles can be captured through their scattering with the nuclei within the sun, mostly hydrogen and helium. This makes the core of the sun a region with concentrated DM distribution. The BDM flux is

$$\begin{aligned} {\varPhi }= f \frac{A}{4\pi D^2}, \end{aligned}$$
(24)

where A is the annihilation rate, and \(D = 1\,\mathrm{AU}\) is the distance from the sun. f is a model-dependent parameter, where \(f = 2\) for two-component DM as considered here.

For the parameter space of interest, assuming that the DM annihilation cross section is not too small, the DM distribution in the sun has reached an equilibrium between capture and annihilation. This helps to eliminate the annihilation cross section dependence in our study. The chain of processes involved in giving rise to the boosted DM signal from the sun is illustrated in Fig. 23.

Fig. 23
figure 23

The chain of processes leading to boosted DM signal from the sun. The semi-annihilation and two-component DM models refer to the two examples of the non-minimal dark-sector scenarios introduced in the beginning of Sect. 8. DM\('\) denotes the lighter DM in the two-component DM model. X is a lighter dark sector particle that may decay away

Two additional comments are in order. First, the DM particles cannot be too light, i.e., lighter than 4 GeV [217, 218], otherwise we will lose most of the captured DM through evaporation rather than annihilation; this would dramatically reduce the BDM flux. Additionally, one needs to check that BDM particles cannot lose energy and potentially be recaptured by scattering with the solar material when they escape from the core region after production. Rescattering is found to be rare for the benchmark models considered in this study and we consider the BDM flux to be monochromatic at its production energy.

The event rate to be observed at DUNE is

$$\begin{aligned} R = {\varPhi }\times \sigma _{{\mathrm{SM}} - \chi } \times \varepsilon \times N_T, \end{aligned}$$
(25)

where \({\varPhi }\) is the flux given by Eq. (24), \(\sigma _{\mathrm{SM} - \chi }\) is the scattering cross section of the BDM off of SM particles, \(\varepsilon \) is the efficiency of the detection of such a process, and \(N_T\) is the number of target particles in DUNE. The computation of the flux of BDM from the sun can be found in [179].

The processes of typical BDM scattering in argon are illustrated in Fig. 24. We generate the signal events and calculate interaction cross sections in the detector using a newly developed BDM module [12, 13, 219] that includes elastic and deep inelastic scattering, as well as a range of nuclear effects. This conservative event generation neglects the dominant contributions from baryon resonances in the final state hadronic invariant mass range of 1.2–1.8 GeV, which should not have a major effect on our main results. The interactions are taken to be mediated by an axial, flavor-universal \(Z^\prime \) coupling to both the BDM and with the quarks. The axial charge is taken to be 1. The events are generated for the \(10 \, \text {kt}\) DUNE detector module [220], though we only study the dominant scattering off of the \({}^{40}\)Ar atoms therein. The method for determining the efficiency \(\varepsilon \) is described below. The number of target argon atoms is \(N = 1.5 \times 10^{32}\) assuming a target mass of \(10 \, \text {kt}\) .

Fig. 24
figure 24

Diagram illustrating each of the three processes contributing to dark matter scattering in argon: elastic (left), baryon resonance (middle), and deep inelastic (right)

The main background in this process comes from the NC interactions of atmospheric neutrinos and argon, as they share the features that the timing of events is unknown in advance (unlike events of neutrinos produced by the accelerator), and that the interactions with argon produce hadronic activity in the detector. We use GENIE to generate the NC atmospheric neutrino events. This simulation predicts 845 events in a \(10 \, \text {kt}\) module for one year of exposure.

The finite detector resolution is taken into account by smearing the direction of the stable final state particles, including protons, neutrons, charged pions, muons, electrons, and photons, with the expected angular resolution, and by ignoring the ones with kinetic energy below detector threshold, using the parameters reported in the DUNE CDR [3]. We form as the observable the total momentum from all the stable final state particles, and obtain its angle with respect to the direction of the sun. The sun position is simulated with the SolTrack package [221] including the geographical coordinates of the DUNE FD. We consider both the scenarios in which we can reconstruct neutrons, according to the parameters described in the DUNE CDR, and in which neutrons will not be reconstructed at all. Figure 25 shows the angular distributions of the BDM signals with mass of 10 GeV and different boost factors, and of the background events.

Fig. 25
figure 25

Angular distribution of the BDM signal events for a BDM mass of 10 GeV and different boosted factors, \(\gamma \), and of the atmospheric neutrino NC background events. \(\theta \) represents the angle of the sum over all the stable final state particles as detailed in the text. The amount of background represents 1-year data collection, magnified by a factor 100, while the amount of signal reflects the detection efficiency of 10,000 MC events. The top plot shows the scenario where neutrons can be reconstructed, while the bottom plot represents the scenario without neutrons

To increase the signal fraction in our samples, we select events with \(\cos \theta > 0.6\), and obtain the selection efficiency \(\varepsilon \) for different BDM models. We predict that \(104.0 \pm 0.7\) and \(79.4 \pm 0.6\) background events per year, in the scenarios with and without neutrons respectively, survive the selection in a DUNE \(10 \, \text {kt}\) module.

Fig. 26
figure 26

Expected \(5\sigma \) discovery reach with one year of DUNE livetime for one \(10 \, \text {kt}\) module including neutrons in reconstruction (top) and excluding neutrons (bottom)

The resulting expected sensitivity is presented in Fig. 26 in terms of the DM mass and the \(Z^\prime \) gauge coupling for potential DM boosts of \(\gamma = 1.25,2,10\) and for a fixed mediator mass of \(M_{Z^\prime } = 1~\mathrm{GeV}\). We assume a DUNE livetime of one year for one \(10 \, \text {kt}\) module. The models presented here are currently unconstrained by direct detection searches if the thermal relic abundance of the DM is chosen to fit current observations. Figure 27 compares the sensitivity of 10 years of data collected in DUNE (\({40} \, \text {kt}\)) to re-analyses of the results from other experiments, including Super Kamiokande [222] and DM direct detection, PICO-60 [223] and PandaX [224]. An extension to this study can be found in Ref. [225].

Fig. 27
figure 27

Comparison of sensitivity of DUNE for 10 years of data collection and \({40} \, \text {kt}\) of detector mass with Super Kamiokande, assuming 10 and 100% of the selection efficiency on the atmospheric neutrino analysis in Ref. [222], and with the reinterpretations of the current results from PICO-60 [223] and PandaX [224]. The samples with two boosted factors, \(\gamma = 1.25\) (top) and \(\gamma = 10\) (bottom), are also presented

8.4 Summary of dark matter detection prospects

We have conducted simulation studies of the dark matter models described in Eqs. (14) and (15) in terms of their detection prospects at the DUNE ND and FD. Thanks to its relatively low threshold and strong particle identification capabilities, DUNE presents an opportunity to significantly advance the search for LDM and BDM beyond what has been possible with water Cherenkov detectors.

In the case of the ND, we assumed that the relativistic DM is being produced directly at the target and leaves an experimental signature through an elastic electron scattering. Using two constrained parameters of the light DM model and a range of two free parameters, a sensitivity map was produced. Within the context of the vector portal DM model and the chosen parameter constraints along with the electron scattering as the signal event, this result sets stringent limits on DM parameters that are comparable or even better than recent experimental bounds in the sub-GeV mass range.

By contrast, in the case of the FD modules, we assumed that the signal events are due to DM coming from the galactic halo and the sun with a significant boost factor. For the inelastic scattering case, the DM scatters off either an electron or proton in the detector material into a heavier unstable dark-sector state. The heavier state, by construction, decays back to DM and an electron-positron pair via a dark-photon exchange. Therefore, in the final state, a signal event comes with an electron or proton recoil plus an electron-positron pair. This distinctive signal feature enabled us to perform (almost) background-free analyses.

As ProtoDUNE detectors are prototypes of DUNE FD modules, the same study was conducted [186] and corresponding results were compared with the ones of the DUNE FD modules. We first investigated the experimental sensitivity in a dark-photon parameter space, dark-photon mass \(M_V\) versus kinetic mixing parameter \(\epsilon \). The results are shown separately for Scenarios 1 and 2 in Fig. 21. They suggest that DUNE FD modules would probe a broad range of unexplored regions; they would allow for reaching \(\sim 1-2\) orders of magnitude smaller \(\epsilon \) values than the current limits along MeV to sub-GeV-range dark photons. We also examined model-independent reaches at DUNE FD modules, providing limits for models that assume the existence of iBDM (or iBDM-like) signals (i.e., a target recoil and a fermion pair).

For the elastic scattering case, we considered the case in which BDM comes from the sun. With one year of data, the \(5\sigma \) sensitivity is expected to reach a coupling of \(g_{Z^\prime }^4 = 9.57 \times 10^{-10}\) for a boost of 1.25 and \(g_{Z^\prime }^4 = 1.49 \times 10^{-10}\) for a boost of 10 at a DM mass of \(10 \, \text {GeV}\) without including neutrons in the reconstruction.

9 Baryon number violating processes

Unifying three of the fundamental forces in the universe, the strong, electromagnetic, and weak interactions, is a shared goal for the current world-wide program in particle physics. Grand unified theories (GUTs), extending the SM to include a unified gauge symmetry at very high energies (more than \(10^{15}~\text {GeV}\)), predict a number of observable effects at low energies, such as nucleon decay [226,227,228,229,230]. Since the early 1980s, supersymmetric GUT models were preferred for a number of reasons, including gauge-coupling unification, natural embedding in superstring theories, and their ability to solve the fine-tuning problem of the SM. Supersymmetric GUT models [231,232,233,234,235,236,237,238,239] generically predict that the dominant proton decay mode is \(p\rightarrow K^+ {\overline{\nu }}\), in contrast to non-supersymmetric GUT models, which typically predict the dominant decay mode to be \(p \rightarrow e^+ \pi ^0\). Although the LHC did not find any evidence for supersymmetry (SUSY) at the electroweak scale, as was expected if SUSY were to solve the gauge hierarchy problem in the SM, the appeal of a GUT still remains. In particular, gauge-coupling unification can still be achieved in non-supersymmetric GUT models by the introduction of one or more intermediate scales (see, for example, [240]). Several experiments have sought signatures of nucleon decay, with the best limits for most decay modes set by the Super-Kamiokande experiment [241,242,243], which features the largest sensitive mass and exposure to date.

The excellent imaging, as well as calorimetric and particle identification capabilities, of the LArTPC technology implemented for the DUNE FD will exploit a number of complementary signatures for a broad range of baryon-number violating processes. Should nucleon decay rates lie just beyond current limits, observation of even one or two candidate events with negligible background could constitute compelling evidence. In the DUNE era, two other large detectors, Hyper-Kamiokande [244] and JUNO [245] will be conducting nucleon decay searches. Should a signal be observed in any single experiment, confirmation from experiments using different detector technologies and nuclear targets, and therefore subject to different backgrounds, would be very powerful.

Neutron-antineutron (\(n-{\bar{n}}\)) oscillation is a baryon number violating process that has never been observed but is predicted by a number of BSM theories [246]. In this context, baryon number conservation is an accidental symmetry rather than a fundamental one, which means baryon number violation does not stand against the fundamental gauge symmetries. Discovering baryon number violation would have implications on the source of matter-antimatter symmetry in our universe given Sakharov’s conditions for such asymmetry to arise [247]. In particular, the neutron-antineutron oscillation (\(n-{\bar{n}}\)) process violates baryon number by two units and, therefore, could also have further implications for the smallness of neutrino masses [246]. Since the \(n-{\bar{n}}\) transition operator is a six-quark operator, of dimension 9, with a coefficient function of dimension (mass)\(^{-5}\), while the proton decay operator is a four-fermion operator, of dimension 6, with a coefficient function of dimension (mass)\(^{-2}\), one might naively assume that \(n-{\bar{n}}\) oscillations would always be suppressed relative to proton decay as a manifestation of baryon number violation. However, this is not necessarily the case; indeed, there are models [248,249,250,251] in which proton decay is very strongly suppressed down to an unobservably small level, while \(n-{\bar{n}}\) oscillations occur at a level comparable to present limits. This shows the value of a search for \(n-{\bar{n}}\) transitions at DUNE. Searches for this process using both free neutrons and nucleus-bound neutron states have been carried out since the 1980s. The current best 90% CL limits on the (free) neutron oscillation lifetime are \(8.6 \times 10^{7} \, \text {s}\) from free \(n-{\bar{n}}\) searches and \(2.7 \times 10^{8} \, \text {s}\) from nucleus-bound \(n-{\bar{n}}\) searches [252, 253]. As with nucleon decay, searches for \(n-{\bar{n}}\) oscillations performed by DUNE and those performed by Super-Kamiokande, Hyper-Kamiokande, and the European Spallation Source [246] are highly complementary. Should a signal be observed in any one experiment, confirmation from another experiment with a different detector technology and backgrounds would be very powerful.

9.1 Event simulation and reconstruction

To estimate the sensitivity to baryon number violation in DUNE, simulation of both signal and background events is performed using GENIE version 2.12.10. For nucleon decay, a total of 68 single-nucleon exclusive decay channels listed in the 2016 update of the PDG [135] are available in GENIE. The list includes two-, three-, and five-body decays. If a bound nucleon decays, the remaining nucleus can be in an excited state and will typically de-excite by emitting nuclear fission fragments, nucleons, and photons. At present, de-excitation photon emission is simulated only for oxygen. The simulation of neutron-antineutron oscillation was developed [254] and implemented in GENIE. Implementing this process in GENIE used GENIE’s existing modeling of Fermi momentum and binding energy for both the oscillating neutron and the nucleon with which the resulting antineutron annihilates. Once a neutron has oscillated to an antineutron in a nucleus, the antineutron has a 18/39 chance of annihilating with a proton in argon, and a 21/39 chance of annihilating with a neutron. The energies and momenta of the annihilation products are assigned randomly but consistently with four-momentum conservation. The products of the annihilation process follow the branching fractions (shown in Table 9) measured in low-energy antiproton annihilation on hydrogen [254].

The default model in GENIE for the propagation of particles inside the nucleus is hA2015, an empirical, data-driven model that does not simulate the cascade of hadronic interactions step by step, but instead uses one effective interaction to represent the effect of final-state interactions (FSI). Hadron-nucleus scattering data is used to tune the predictions.

The dominant background for these searches is from atmospheric neutrino interactions. Backgrounds from neutrino interactions are simulated with GENIE, using the Bartol model of atmospheric neutrino flux [255]. To estimate the event rate, we integrate the product of the neutrino flux and interaction cross section. Table 8 shows the event rate for different neutrino species for an exposure of \(10~\text {kt} \, \cdot \, \text {year} \), where oscillation effects are not included. To suppress atmospheric neutrino background to the level of one event per \(\text {Mt} \, \cdot \, \text {year}\), which would yield 0.4 events after ten years of operation with a \({40} \, \text {kt}\) fiducial volume, the necessary background rejection is \(1 - (1/288600) = 1 - 3\times 10^{-6} = 0.999997\), where background rejection is defined as the fraction of background that is not selected.

Table 8 Expected rate of atmospheric neutrino interactions in \({}^{40}\)Ar for a \({10} \, \text {kt} \, \cdot \, \text {year}\) exposure (not including oscillations)

These analyses assume that the detector is successfully triggered on all signal events, and that the PD system correctly determines the event start time (\(t_0\)). Two distinct methods of reconstruction and event selection have been applied in these analyses. One employs 3D track and vertex reconstruction provided by Projection Matching Algorithm (PMA) [1], a standard DUNE reconstruction algorithm. PMA was designed to address transformation from a set of independently reconstructed 2D projections of objects into a 3D representation. This algorithm uses clusters of hits from 2D pattern recognition as its input. The other reconstruction method involves image classification of 2D images of reconstructed hits using a Convolutional Neural Network (CNN). The two methods are combined in the form of a multivariate analysis, which uses the image classification score with other physical observables extracted from traditional reconstruction.

9.2 Nucleon decay

Because of the already stringent limits set by Super-Kamiokande on \(p \rightarrow e^+ \pi ^0\) and the unique ability to track and identify kaons in a LArTPC, the initial nucleon decay studies in DUNE have focused on nucleon decay modes featuring kaons, in particular \(p\rightarrow K^+ {\overline{\nu }}\). The experimental signature of this channel is a single \(K^{+}\) originating inside the fiducial volume of the detector. The kaon typically stops and decays at rest with a lifetime of 12 ns. The most common decay mode, \(K^{+} \rightarrow \mu ^{+}\nu _{\mu }\), results in a monoenergetic muon with momentum of 236 MeV/c. In the next most probable decay, \(K^{+} \rightarrow \pi ^{+}\pi ^{0}\), the two pions are produced back to back. In a water Cherenkov detector, the kaon is typically below Cherenkov threshold, and only the kaon decay products are observed. In DUNE’s LArTPC, the kaon can be detected and identified by its distinctive dE/dx signature, as well as by its decay [256].

For a proton decay at rest, the outgoing kaon is monoenergetic with kinetic energy of \({105} \, \text {MeV}\) and momentum of \({339} \, \text {MeV}/\hbox {c}\). In bound proton decay, the momentum of the kaon is smeared by the Fermi motion of the protons inside the nucleus. FSI between the outgoing kaon and the residual nucleus may reduce the kaon momentum, and may also modify the final state, by ejecting nucleons for example. Protons ejected from the nucleus can obscure the dE/dx measurement of the kaon if the tracks overlap. The \(K^{+}\) may also charge exchange, resulting in a \(K^{0}\) in the final state. The \(K^{+}\) cannot be absorbed due to strangeness conservation and the lack of \(S=1\) baryons. The residual nucleus may also be in an excited state, producing de-excitation photons.

The main backgrounds in nucleon decay searches are interactions of atmospheric neutrinos. For \(p\rightarrow K^+ {\overline{\nu }}\), the background is neutrino interactions that mimic a single \(K^{+}\) and its decay products. Because the kaon is not detected in a water Cherenkov detector, neutrino interactions that produce a single \(K^{+}\) and no other particles above Cherenkov threshold are an irreducible background. This includes charged-current reactions like the Cabibbo-suppressed \(\nu _{\mu } n \rightarrow \mu ^{-} K^{+} n\), where the final-state muon and kaon are below threshold, as well as neutral-current processes like \(\nu p \rightarrow \nu K^{+} {\varLambda }\) followed by \({\varLambda }\rightarrow p \pi ^{-}\) where the \({\varLambda }\) decay products are below threshold. Strangeness is always conserved in neutral-currents, so kaons produced in NC interactions are always accompanied by a hyperon or another kaon. Water Cherenkov detectors and liquid scintillator detectors like JUNO can also detect neutron captures, which provide an additional handle on backgrounds, many of which have final-state neutrons. However, neutrons can also be present in \(p\rightarrow K^+ {\overline{\nu }}\) signal due to FSI, and the rate of nucleon ejection in kaon-nucleus interactions is not well understood. Nuclear de-excitation photons are also typically produced, but these are similar in both proton decay and atmospheric neutrino events. In the Super-Kamiokande analysis of \(p\rightarrow K^+ {\overline{\nu }}\) the time difference between the de-excitation photons from the oxygen nucleus and the muon from kaon decay was found to be an effective way to reduce backgrounds [241]. In JUNO, the three-fold time coincidence between the kaon, the muon from the kaon decay, and the electron from the muon decay is expected to be an important discriminant between signal and background [245].

The possibility of using the time difference between the kaon scintillation signal and the scintillation signal from the muon from the kaon decay has been investigated in DUNE. Studies indicate that measuring time differences on the scale of the kaon lifetime (12 ns) is difficult in DUNE, independent of photon detector acceptance and timing resolution, due to both the scintillation process in argon - consisting of fast (ns-scale) and slow (\(\mu \)s-scale) components - and Rayleigh scattering over long distances.

In a LArTPC, a charged particle traveling just a few cm can be detected, and the other particles produced in association with a kaon by atmospheric neutrinos are generally observed. However, with FSI the signal process can also include final-state protons, so requiring no other final-state particles will reject some signal events. Furthermore, \(\nu _{\mu }\) charged-current quasi-elastic scattering (CCQE), \(\nu _{\mu } n \rightarrow \mu ^{-} p\), can mimic the \(K^{+} \rightarrow \mu ^{+} \nu _{\mu }\) decay when the proton is mis-reconstructed as a kaon.

The kaon reconstruction is especially challenging for very short tracks, which may traverse only a few wires. The dE/dx signature in signal events can be obscured by additional final-state protons that overlap with the start of the kaon track. Without timing resolution sufficient to resolve the 12 ns kaon lifetime, the dE/dx profile is the only distinguishing feature. The background from atmospheric neutrino events without true final-state kaons, which is important given the presence of FSI, was neglected in previous estimates of \(p\rightarrow K^+ {\overline{\nu }}\) sensitivity in LArTPC [257].

Other backgrounds, such as those initiated by cosmic-ray muons, can be controlled by requiring no activity close to the edges of the time projection chambers (TPCs) and by stringent single kaon identification within the energy range of interest [77, 258].

FSI significantly modify the observable distributions in the detector. For charged kaons, the hA2015 model includes only elastic scattering and nucleon knock-out, tuned to \(K^{+}-\)C data [259, 260]. Charge exchange is not included, nor are strong processes that produce kaons inside the nucleus, such as \(\pi ^{+}n \rightarrow K^{+} {\varLambda }\). Figure 28 shows the kinetic energy of a kaon from \(p\rightarrow K^+ {\overline{\nu }}\) before and after FSI as simulated with hA2015. Kaon interactions always reduce the kaon energy, and the kaon spectrum becomes softer on average with FSI. Of the kaons, 31.5% undergo elastic scattering resulting in events with very low kinetic energy; 25% of kaons have a kinetic energy of \(\le {50} \, \text {MeV}\). When the kaon undergoes elastic scattering, a nucleon can be knocked out of the nucleus. Of decays via this channel, 26.7% have one neutron coming from FSI, 15.3% have at least one proton, and 10.3% have two protons coming from FSI. These secondary nucleons are detrimental to reconstructing and selecting \(K^{+}\).

Fig. 28
figure 28

Kinetic energy of kaons in simulated proton decay events, \(p\rightarrow K^+ {\overline{\nu }}\), in DUNE. The kinetic energy distribution is shown before and after final state interactions in the argon nucleus

Other FSI models include the full cascade, and predict slightly different final states, but existing data lack power to favor one model over another. MINERvA has measured the differential cross section for charged-current \(K^{+}\) production by neutrinos on plastic scintillator (CH) as a function of kaon energy, which is sensitive to FSI, and shows a weak preference for the GENIE hA2015 FSI model over a prediction with no FSI [261]. Compared to the kaon energy spectrum measured by MINERVA, FSI have a much larger impact on \(p\rightarrow K^+ {\overline{\nu }}\) in argon, and the differences between models are less significant than the overall effect.

The kaon FSI in Super-Kamiokande ’s simulation of \(p\rightarrow K^+ {\overline{\nu }}\) in oxygen seem to have a smaller effect on the outgoing kaon momentum distribution [241] than is seen here with the GENIE simulation on argon. Some differences are expected due to the different nuclei, but differences in the FSI models are under investigation.

Kaon FSI have implications on the ability to identify \(p\rightarrow K^+ {\overline{\nu }}\) events in DUNE. Track reconstruction efficiency for a charged particle \(x^{\pm }\) is defined as

$$\begin{aligned} \epsilon _{x^{\pm }} = \frac{x^{\pm } \text{ particles } \text{ with } \text{ a } \text{ reconstructed } \text{ track }}{\text{ events } \text{ with } x^{\pm } \text{ particle } }. \end{aligned}$$
(26)

The denominator includes events in which an \(x^{\pm }\) particle was created and has deposited energy within any of the TPCs. The numerator includes events in which an \(x^{\pm }\) particle was created and has deposited energy within any of the TPCs, and a reconstructed track can be associated to the \(x^{\pm }\) particle based on the number of hits generated by that particle along the track. This efficiency can be calculated as a function of true kinetic energy and true track length.

Fig. 29
figure 29

Tracking efficiency for kaons in simulated proton decay events, \(p\rightarrow K^+ {\overline{\nu }}\), as a function of kaon kinetic energy (top) and true path length (bottom)

Figure 29 shows the tracking efficiency for \(K^{+}\) from proton decay via \(p\rightarrow K^+ {\overline{\nu }}\) as a function of true kinetic energy and true path length. The overall tracking efficiency for kaons from proton decay is 58.0%, meaning that 58.0% of all the simulated kaons are associated with a reconstructed track in the detector. From Fig. 29, the tracking threshold is approximately \(\sim {40} \, \text {MeV}\) of kinetic energy, which translates to \(\sim {4.0} \, \text {cm}\) in true path length. The biggest loss in tracking efficiency is due to kaons with \(<{40} \, \text {MeV}\) of kinetic energy due to scattering inside the nucleus. The efficiency levels off to approximately 80% above \({80} \, \text {MeV}\) of kinetic energy; this inefficiency even at high kinetic energy is due mostly to kaons that decay in flight. Both kaon scattering in the liquid argon (LAr) and charge exchange are included in the detector simulation but are relatively small effects (4.6% of kaons scatter in the LAr and 1.2% of kaons experience charge exchange). The tracking efficiency for muons from the decay of the \(K^{+}\) in \(p\rightarrow K^+ {\overline{\nu }}\) is 90%.

Hits associated with a reconstructed track are used to calculate the energy loss of charged particles, which provides valuable information on particle energy and species. If the charged particle stops in the LArTPC active volume, a combination of dE/dx and the reconstructed residual range (R, the path length to the end point of the track) is used to define a parameter for particle ID (PID). The parameter, PIDA, is defined as [262]

$$\begin{aligned} PIDA = \left\langle \left( \frac{dE}{dx}\right) _{i}R^{0.42}_{i}\right\rangle , \end{aligned}$$
(27)

where the median is taken over all track points i for which the residual range \(R_i\) is less than \({30} \, \text {cm}\).

Fig. 30
figure 30

Particle identification using PIDA for muons and kaons in simulated proton decay events, \(p\rightarrow K^+ {\overline{\nu }}\), and protons in simulated atmospheric neutrino background events. The curves are normalized by area

Figure 30 shows the PIDA performance for kaons (from proton decay), muons (from kaon decay), and protons produced by atmospheric neutrino interactions. The tail with lower values in each distribution is due to cases where the decay/stopping point was missed by the track reconstruction. The tail with higher values is caused when a second particle overlaps at the decay/stopping point causing higher values of dE/dx and resulting in higher values of PIDA. In addition, ionization fluctuations smear out these distributions.

PID via dE/dx becomes complicated when the reconstructed track direction is ambiguous, in particular if additional energy is deposited at the vertex in events where FSI is significant. The dominant background to \(p\rightarrow K^+ {\overline{\nu }}\) in DUNE is atmospheric neutrino CC quasielastic (QE) scattering, \(\nu _{\mu } n \rightarrow \mu ^{-} p\). When the muon happens to have very close to the \(236 \text {MeV}/c\) momentum expected from a \(K^{+}\) decay at rest and is not captured, it is indistinguishable from the muon resulting from \(p\rightarrow K^+ {\overline{\nu }}\) followed by \(K^{+} \rightarrow \mu ^{+}\nu _{\mu }\). When the proton is also mis-reconstructed as a kaon, this background mimics the signal process.

The most important difference between signal and this background source is the direction of the hadron track. For an atmospheric neutrino, the proton and muon originate from the same neutrino interaction point, and the characteristic Bragg rise occurs at the end of the proton track farthest from the muon-proton vertex. In signal, the kaon-muon vertex is where the \(K^{+}\) stops and decays at rest, so its ionization energy deposit is highest near the kaon-muon vertex. To take advantage of this difference, a log-likelihood ratio discriminator is used to distinguish signal from background. Templates are formed by taking the reconstructed and calibrated energy deposit as a function of the number of wires from both the start and end of the \(K^{+}\) candidate hadron track. Two log-likelihood ratios are computed separately for each track. The first begins at the hadron-muon shared vertex and moves along the hadron track (the “backward” direction). The second begins at the other end of the track, farthest from the hadron-muon shared vertex, moves along the hadron track the other way (the “forward” direction). For signal events, this effectively looks for the absence of a Bragg rise at the \(K^{+}\) start, and the presence of one at the end, and vice versa for background. At each point, the probability density for signal and background, \(P^{sig}\) and \(P^{bkg}\), are determined from the templates. Forward and backward log-likelihood ratios are computed as

$$\begin{aligned} {\mathcal {L}}_{fwd(bkwd)} = \sum _{i} \log \frac{P^{sig}_i}{P^{bkg}_i}, \end{aligned}$$
(28)

where the summation is over the wires of the track, in either the forward or backward direction. Using either the forward or backward log-likelihood ratio alone gives some discrimination between signal and background, but using the sum gives better discrimination. While the probability densities are computed based on the same samples, defining one end of the track instead of the other as the vertex provides more information. The discriminator is the sum of the forward and backward log-likelihood ratios:

$$\begin{aligned} {\mathcal {L}} = {\mathcal {L}}_{fwd} + {\mathcal {L}}_{bkwd}. \end{aligned}$$
(29)

Applying this discriminator to tracks with at least ten wires gives a signal efficiency of roughly 0.4 with a background rejection of 0.99.

A Boosted Decision Tree (BDT) classifier is used for event selection in the analysis presented here. The software package Toolkit for Multivariate Data Analysis with ROOT (TMVA4) [263] is used with AdaBoost as the boosted algorithm. The BDT is trained on a sample of MC events (50,000 events for signal and background) that is statistically independent from the sample of MC events used in the analysis (approximately 100,000 events for signal and 600,000 events for background). Image classification using a CNN is performed using 2D images of DUNE MC events. The image classification provides a single score value as a metric of whether any given event is consistent with a proton decay, and this score can be used as a powerful discriminant for event identification. In the analysis presented here, the CNN technique alone does not discriminate between signal and background as well as a BDT, so the CNN score is used as one of the input variables to the BDT in this analysis. The other variables in the BDT include numbers of reconstructed objects (tracks, showers, vertices), variables related to visible energy deposition, PID variables [PIDA, Eq. (27), and \({\mathcal {L}}\), Eq. (29)], reconstructed track length, and reconstructed momentum. Figure 31 shows the distribution of the BDT output for signal and background. Backgrounds from atmospheric neutrinos are weighted by the oscillation probability in the BDT input distributions.

Fig. 31
figure 31

Boosted Decision Tree response for \(p\rightarrow K^+ {\overline{\nu }}\) for signal (blue) and background (red)

Figure 32 shows a \(p\rightarrow K^+ {\overline{\nu }}\) signal event. The event display shows the reconstructed kaon track in green and the reconstructed muon track from the kaon decay in red; hits from the Michel electron coming from the muon decay can be seen at the end of the muon track. Figure 33 shows an event with a similar topology produced by an atmospheric neutrino interaction, \(\nu _{\mu } n \rightarrow \mu ^{-} p\). This type of event can be selected in the \(p\rightarrow K^+ {\overline{\nu }}\) sample if the proton is misidentified as a kaon. Hits associated with the reconstructed muon track are shown in red, and hits associated with the reconstructed proton track are shown in green. Hits from the decay electron can be seen at the end of the muon track.

Fig. 32
figure 32

Event display for an easily recognizable \(p\rightarrow K^+ {\overline{\nu }}\) signal event. The vertical axis is TDC value, and the horizontal axis is wire number. The bottom view is induction plane one, the middle is induction plane two, and the top is the collection plane. Hits associated with the reconstructed muon track are shown in red, and hits associated with the reconstructed kaon track are shown in green. Hits from the decay electron can be seen at the end of the muon track

Fig. 33
figure 33

Event display for an atmospheric neutrino interaction, \(\nu _{\mu } n \rightarrow \mu ^{-} p\), which might be selected in the \(p\rightarrow K^+ {\overline{\nu }}\) sample if the proton is misidentified as a kaon. The vertical axis is TDC value, and the horizontal axis is wire number. The bottom view is induction plane one, the middle is induction plane two, and the top is the collection plane. Hits associated with the reconstructed muon track are shown in red, and hits associated with the reconstructed proton track are shown in green. Hits from the decay electron can be seen at the end of the muon track

The proton decay signal and atmospheric neutrino background events are processed using the same reconstruction chain and subject to the same selection criteria. There are two preselection cuts to remove obvious background. One cut requires at least two tracks, which aims to select events with a kaon plus a kaon decay product (usually a muon). The other cut requires that the longest track be less than \({100} \, \text {cm}\); this removes backgrounds from high energy neutrino interactions. After these cuts, 50% of the signal and 17.5% of the background remain in the sample. The signal inefficiency at this stage of selection is due mainly to the kaon tracking efficiency. Optimal lifetime sensitivity is achieved by combining the preselection cuts with a BDT cut that gives a signal efficiency of 0.15 and a background rejection of 0.999997, which corresponds to approximately one background event per \(\text {Mt} \, \cdot \, \text {year}\).

The limiting factor in the sensitivity is the kaon tracking efficiency. The reconstruction is not yet optimized, and the kaon tracking efficiency should increase with improvements in the reconstruction algorithms. To understand the potential improvement, a visual scan of simulated decays of kaons into muons was performed. For this sample of events, with kaon momentum in the \(150 \, {\text {MeV}}/c\) to \(450 \, {\text {MeV}}/c\) range, scanners achieved greater than 90% efficiency at recognizing the \(K^{+} \rightarrow \mu ^{+} \rightarrow e^{+}\) decay chain. The inefficiency came mostly from short kaon tracks (momentum below \(180 \, \text {MeV}/c\)) and kaons that decay in flight. Note that the lowest momentum kaons (\(< 150 \, \text {MeV}/c\)) were not included in the study; the path length for kaons in this range would also be too short to track. Based on this study, the kaon tracking efficiency could be improved to a maximum value of approximately 80% with optimized reconstruction algorithms, where the remaining inefficiency comes from low-energy kaons and kaons that charge exchange, scatter, or decay in flight. Combining this tracking performance improvement with some improvement in the K/p separation performance for short tracks, the overall signal selection efficiency improves from 15% to approximately 30%.

The analysis presented above is inclusive of all possible modes of kaon decay; however, the current version of the BDT preferentially selects kaon decay to muons, which has a branching fraction of roughly 64%. The second most prominent kaon decay is \(K^{+} \rightarrow \pi ^{+}\pi ^0\), which has a branching fraction of 21%. Preliminary studies that focus on reconstructing a \(\pi ^{+}\pi ^0\) pair with the appropriate kinematics indicate that the signal efficiency for kaons that decay via the \(K^{+} \rightarrow \pi ^{+}\pi ^0\) mode is approximately the same as the signal efficiency for kaons that decay via the \(K^{+} \rightarrow \mu ^{+}\nu _{\mu }\) mode. This assumption is included in our sensitivity estimates below.

Because the DUNE efficiency to reconstruct a kaon track is strongly dependent on the kaon kinetic energy as seen in Fig. 29, the FSI model is an important source of systematic uncertainty. To account for this uncertainty, kaon-nucleon elastic scattering (\(K^{+}p(n)\rightarrow K^{+}p(n)\)) is re-weighted by \(\pm {50}\%\) in the simulation. The absolute uncertainty on the efficiency with this re-weighting is 2%, which is taken as the systematic uncertainty on the signal efficiency. The dominant uncertainty in the background is due to the absolute normalization of the atmospheric neutrino rate. The Bartol group has carried out a detailed study of the systematic uncertainties, where the absolute neutrino fluxes have uncertainties of approximately 15% [264]. The remaining uncertainties are due to the cross section models for neutrino interactions. The uncertainty on the CC0\(\pi \) cross section in the energy range relevant for these backgrounds is roughly 10% [265]. Based on these two effects, a conservative 20% systematic uncertainty in the background is estimated.

With a 30% signal efficiency and an expected background of one event per \(\hbox {Mt} \, \cdot \, \hbox {year}\), a 90% CL lower limit on the proton lifetime in the \(p\rightarrow K^+ {\overline{\nu }}\) channel of \(1.3 \times 10^{34} \, \text {years}\) can be set, assuming no signal is observed over ten years of running with a total of \({40} \, \text {kt}\) of fiducial mass. This calculation assumes constant signal efficiency and background rejection over time and for each of the FD modules. Additional running improves the sensitivity proportionately if the experiment remains background-free.

Another potential mode for a baryon number violation search is the decay of the neutron into a charged lepton plus meson, i.e., \(n\rightarrow e^{-}K^{+}\). In this mode, \({\varDelta }B = -{\varDelta }L\), where B is baryon number and L is lepton number. The current best limit on this mode is \(3.2 \times 10^{31} \, \text {years}\) from the FREJUS collaboration [266]. The reconstruction software for this analysis is the same as for the \(p\rightarrow K^+ {\overline{\nu }}\) analysis; the analysis again uses a BDT that includes an image classification score as an input. To calculate the lifetime sensitivity for this decay mode the same systematic uncertainties and procedures are used. The selection efficiency for this channel including the expected tracking improvements is 0.47 with a background rejection of 0.99995, which corresponds to 15 background events per \(\hbox {Mt} \, \cdot \, \hbox {year}\). The lifetime sensitivity for a \({400} \, \text {kt} \, \cdot \, \text {year}\) exposure is \(1.1 \times 10^{34} \, \text {years}\).

9.3 Neutron–antineutron oscillation

Neutron-antineutron oscillations can be detected via the subsequent antineutron annihilation with a neutron or a proton. Table 9 shows the effective branching ratios for the antineutron annihilation modes applicable to intranuclear searches, modified from [253]. It is known that other, more fundamentally consistent branching fractions exist [267, 268], but the effects of these on final states is believed to be minimal. The annihilation event will have a distinct, roughly spherical signature of a vertex with several emitted light hadrons (a so-called “pion star”), with total energy of twice the nucleon mass and roughly zero net momentum. Reconstructing these hadrons correctly and measuring their energies is key to identifying the signal event. The main background for these \(n-{\bar{n}}\) annihilation events is caused by atmospheric neutrinos. As with nucleon decay, nuclear effects and FSI make the picture more complicated. As shown in Table  9, every decay mode contains at least one charged pion and one neutral pion. The pion FSI in the hA2015 model in GENIE include pion elastic and inelastic scattering, charge exchange and absorption.

Table 9 Effective branching ratios for antineutron annihilation in \({}^{40}\)Ar, as implemented in GENIE

Figure 34 shows the momentum distributions for charged and neutral pions before FSI and after FSI. These distributions show the FSI makes both charged and neutral pions less energetic. The effect of FSI on pion multiplicity is also rather significant; 0.9% of the events have no charged pions before FSI, whereas after FSI 11.1% of the events have no charged pions. In the case of the neutral pion, 11.0% of the events have no neutral pions before FSI, whereas after FSI, 23.4% of the events have no neutral pions. The decrease in pion multiplicity is primarily due to pion absorption in the nucleus. Another effect of FSI is nucleon knockout from pion elastic scattering. Of the events, 94% have at least one proton from FSI and 95% of the events have at least one neutron from FSI. Although the kinetic energy for these nucleons peak at a few tens of \(\text {MeV}\), the kinetic energy can be as large as hundreds of \(\text {MeV}\). In summary, the effects of FSI in \(n-{\bar{n}}\) become relevant because they modify the kinematics and topology of the event. For instance, even though the decay modes of Table 9 do not include nucleons in their decay products, nucleons appear with high probability after FSI.

Fig. 34
figure 34

Top: momentum of an individual charged pion before and after final state interactions. Bottom: momentum of an individual neutral pion before and after final state interactions

A BDT classifier is used. Ten variables are used in the BDT as input for event selection, including number of reconstructed tracks and showers, variables related to visible energy deposition, PIDA and dE/dx, reconstructed momentum, and CNN score. Figure 35 shows the distribution of the BDT output for signal and background.

Fig. 35
figure 35

Boosted Decision Tree response for \(n-{\bar{n}}\) oscillation for signal (blue) and background (red)

Figure 36 shows an \(n-{\bar{n}}\) signal event, \(n {\bar{n}} \rightarrow n \pi ^0 \pi ^0 \pi ^{+} \pi ^{-}\). Hits associated with the back-to-back tracks of the charged pions are shown in red. The remaining hits are from the showers from the neutral pions, neutron scatters, and low-energy de-excitation gammas. The topology of this event is consistent with charged pion and neutral pion production. Figure 37 shows an event with a similar topology produced by a NC DIS atmospheric neutrino interaction. This background event mimics the signal topology by having multi-particle production and an electromagnetic shower.

Fig. 36
figure 36

Event display for an \(n-{\bar{n}}\) signal event, \(n {\bar{n}} \rightarrow n \pi ^0 \pi ^0 \pi ^{+} \pi ^{-}\). The vertical axis is TDC value, and the horizontal axis is wire number. The bottom view is induction plane one, the middle is induction plane two, and the top is the collection plane. Hits associated with the back-to-back tracks of the charged pions are shown in red. The remaining hits are from the showers from the neutral pions, neutron scatters, and low-energy de-excitation gammas

Fig. 37
figure 37

Event display for a NC DIS interaction initiated by an atmospheric neutrino. The vertical axis is TDC value, and the horizontal axis is wire number. The bottom view is induction plane one, the middle is induction plane two, and the top is the collection plane. This event mimics the \(n-{\bar{n}}\) signal topology by having multi-particle production and electromagnetic showers

The sensitivity to the \(n-{\bar{n}}\) oscillation lifetime can be calculated for a given exposure, the efficiency of selecting signal events, and the background rate along with their uncertainties. The lifetime sensitivity is obtained at 90% CL for the bound neutron. Then, the lifetime sensitivity for a free neutron is acquired using the conversion from nucleus bounded neutron to free neutron \(n-{\bar{n}}\) oscillation [269]. The uncertainties on the signal efficiency and background rejection are conservatively estimated to be 25%. A detailed evaluation of the uncertainties is in progress.

The free \(n-{\bar{n}}\) oscillation lifetime, \(\tau _{n-{\bar{n}}}\), and bounded \(n-{\bar{n}}\) oscillation lifetime, \(T_{n-{\bar{n}}}\), are related to each other through the intranuclear suppression factor R as

$$\begin{aligned} \tau ^{2}_{n-{\bar{n}}} = \frac{T_{n-{\bar{n}}}}{R}. \end{aligned}$$
(30)

The suppression factor R varies for different nuclei. This suppression factor was calculated for \(^{16}\)O and \(^{56}\)Fe [269]. The R for \(^{56}\)Fe, \(0.666 \times 10^{23} \, {\text {s}}^{-1}\), is used in this analysis for \({}^{40}\)Ar nuclei. More recent work [268] gives a value of R for \({}^{40}\)Ar of \(0.56 \times 10^{23} \, {\text {s}}^{-1}\), which will be applied in future analyses.

The best bound neutron lifetime limit is achieved using a signal efficiency of 8.0% at the background rejection probability of 99.98%, which corresponds to approximately 23 atmospheric neutrino background events for a \({400} \, \text {kt} \, \cdot \, \text {year}\) exposure. The 90% CL limit of a bound neutron lifetime is \(6.45 \times 10^{32} \, \text {years}\) for a \({400} \, \text {kt} \, \cdot \, \text {year}\) exposure. The corresponding limit for the oscillation time of free neutrons is calculated to be \(5.53 \times 10^{8} \, \text {s}\). This is approximately an improvement by a factor of two from the current best limit, which comes from Super-Kamiokande [253].

10 Other BSM physics opportunities

10.1 BSM constraints with tau neutrino appearance

With only 19 \(\nu _{\tau }\)-CC and \({\bar{\nu }}_{\tau }\)-CC candidates detected with high purity, we have less direct experimental knowledge of tau neutrinos than of any other SM particle. Of these, nine \(\nu _{\tau }\)-CC and \({\bar{\nu }}_{\tau }\)-CC candidate events with a background of 1.5 events, observed by the DONuT experiment [270, 271], were directly produced though \(D_S\) meson decays. The remaining 10 \(\nu _{\tau }\)-CC candidate events with an estimated background of two events, observed by the OPERA experiment [272, 273], were produced through the oscillation of a muon neutrino beam. From this sample, a 20% measurement of \({\varDelta }m^{2}_{32}\) was performed under the assumption that \(\sin ^22\theta _{23} = 1\). The Super-Kamiokande and IceCube experiments developed methods to statistically separate samples of \(\nu _{\tau }\)-CC and \({\bar{\nu }}_{\tau }\)-CC events in atmospheric neutrinos to exclude the no-tau-neutrino appearance hypothesis at the 4.6\(\sigma \) level and 3.2\(\sigma \) level respectively [274,275,276], but limitations of Cherenkov detectors constrain the ability to select a high-purity sample and perform precision measurements.

The DUNE experiment has the possibility of significantly improving the experimental situation [277]. Tau-neutrino appearance can potentially improve the discovery potential for sterile neutrinos, NC NSI, and non-unitarity. This channel could also be used as a probe of secret couplings of neutrinos to new light bosons [278]. For model independence, the first goal should be measuring the atmospheric oscillation parameters in the \(\nu _{\tau }\) appearance channel and checking the consistency of this measurement with those performed using the \(\nu _{\mu }\) disappearance channel. A truth-level study of \(\nu _{\tau }\) selection in atmospheric neutrinos in a large, underground LArTPC detector suggested that \(\nu _{\tau }\)-CC interactions with hadronically decaying \(\tau \)-leptons, which make up 65% of total \(\tau \)-lepton decays [135], can be selected with high purity [279]. This analysis suggests that it may be possible to select up to 30% of \(\nu _{\tau }\)-CC events with hadronically decaying \(\tau \)-leptons with minimal neutral-current background. Under these assumptions, we expect to select \({\sim } 25 \,\nu _{\tau }\)-CC candidates per year using the CPV optimized beam. The physics reach of this sample has been studied in Refs. [280, 281]. As shown in Fig. 38 (top), this sample is sufficient to simultaneously constrain \({\varDelta }m^2_{31}\) and \(\sin ^22\theta _{23}\). Independent measurements of \({\varDelta }m^2_{31}\) and \(\sin ^22\theta _{23}\) in the \(\nu _{e}\) appearance, \(\nu _{\mu }\) disappearance, and \(\nu _{\tau }\) appearance channels should allow DUNE to constrain \(|U_{e3}|^2+|U_{\mu 3}|^2+|U_{\tau 3}|^2\) to 6% [280], a significant improvement over current constraints [49].

Fig. 38
figure 38

The 1\(\sigma \) (dashed) and 3\(\sigma \) (solid) expected sensitivity for measuring \({\varDelta }m^2_{31}\) and \(\sin ^2\theta _{23}\) using a variety of samples. Top: The expected sensitivity for 7 years of beam data collection, assuming 3.5 years each in neutrino and antineutrino modes, measured independently using \(\nu _{e}\) appearance (blue), \(\nu _{\mu }\) disappearance (red), and \(\nu _{\tau }\) appearance (green). Adapted from Ref. [280]. Bottom: The expected sensitivity for the \(\nu _{\tau }\) appearance channel using 350 \(\text {kt} \, \cdot \, \text {year}\) of atmospheric exposure

However, all of the events in the beam sample occur at energies higher than the first oscillation maximum due to kinematic constraints. Only seeing the tail of the oscillation maximum creates a partial degeneracy between the measurement of \({\varDelta }m^2_{31}\) and \(\sin ^22\theta _{23}\). Atmospheric neutrinos, due to sampling a much larger L/E range, allow for measuring both above and below the first oscillation maximum with \(\nu _{\tau }\) appearance. Although we only expect to select \({\sim } 70 \, \nu _{\tau }\)-CC and \({\bar{\nu }}_{\tau }\)-CC candidates in \(350~\text {kt} \, \cdot \, \text {year} \) in the atmospheric sample, as shown in Fig. 38 (bottom), a direct measurement of the oscillation maximum breaks the degeneracy seen in the beam sample. The complementary shapes of the beam and atmospheric constraints combine to reduce the uncertainty on \(\sin ^2\theta _{23}\), directly leading to improved unitarity constraints. Finally, a high-energy beam option optimized for \(\nu _{\tau }\) appearance should produce \(\sim \)150 selected \(\nu _{\tau }\)-CC candidates in one year [3]. These higher energy events are further in the tail of the first oscillation maximum, but they will permit a simultaneous measurement of the \(\nu _{\tau }\) cross section. When analyzed within the non-unitarity framework described in Sect. 4, the high-energy beam significantly improves constraints on the parameter \(\alpha _{\tau \tau }\) due to increased matter effects [280].

10.2 Large extra-dimensions

DUNE can search for or constrain the size of large extra-dimensions (LED) by looking for distortions of the oscillation pattern predicted by the three-flavor paradigm. These distortions arise through mixing between the right-handed neutrino Kaluza–Klein modes, which propagate in the compactified extra dimensions, and the active neutrinos, which exist only in the four-dimensional brane [282,283,284]. Such distortions are determined by two parameters in the model, specifically R, the radius of the circle where the extra-dimension is compactified, and \(m_0\), defined as the lightest active neutrino mass (\(m_1\) for normal mass ordering, and \(m_3\) for inverted mass ordering). Searching for these distortions in, for instance, the \(\nu _\mu \) CC disappearance spectrum, should provide significantly enhanced sensitivity over existing results from the MINOS/MINOS+ experiment [285].

Fig. 39
figure 39

Sensitivity to the LED model in Refs. [282,283,284] through its impact on the neutrino oscillations expected at DUNE. For comparison, the MINOS sensitivity [285] is also shown

Figure 39 shows a comparison between the DUNE and MINOS [285] sensitivities to LED at \(90\%\) CL for 2 d.o.f represented by the solid and dashed lines, respectively. In the case of DUNE, an exposure of \(300~\text {kt} \, \cdot \, \text {MW} \cdot \, \text {year} \) was assumed and spectral information from the four oscillation channels, (anti)neutrino appearance and disappearance, were included in the analysis. The muon (anti)neutrino fluxes, cross sections for the neutrino interactions in argon, detector energy resolutions, efficiencies and systematical errors were taken into account by the use of GLoBES files prepared for the DUNE LBL studies. In the analysis, we assumed DUNE simulated data as compatible with the standard three neutrino hypothesis (which corresponds to the limit \(R\rightarrow 0\)) and we have tested the LED model. The solar parameters were kept fixed, and also the reactor mixing angle, while the atmospheric parameters were allowed to float free. In general, DUNE improves over the MINOS sensitivity for all values of \(m_0\) and this is more noticeable for \(m_0\sim 10^{-3}~\text {eV}\), where the most conservative sensitivity limit to R is obtained.

Fig. 40
figure 40

The 90% CL sensitivity regions for dominant mixings \(|U_{e N}|^2\) (top left), \(|U_{\mu N}|^2\) (top right), and \(|U_{\tau N}|^2\) (bottom) are presented for DUNE ND (black) [287]. The regions are a combination of the sensitivity to HNL decay channels with good detection prospects.These are \(N\rightarrow \nu e e\), \(\nu e \mu \), \(\nu \mu \mu \), \(\nu \pi ^0\), \(e \pi \), and \(\mu \pi \). The study is performed for Majorana neutrinos (solid) and Dirac neutrinos (dashed), assuming no background. The region excluded by experimental constraints (grey/brown) is obtained by combining the results from PS191 [288, 289], peak searches [290,291,292,293,294], CHARM [295], NuTeV [296], DELPHI [297], and T2K [298]. The sensitivity for DUNE ND is compared to the predictions of future experiments, SBN [299] (blue), SHiP [300] (red), NA62 [301] (green), MATHUSLA [302] (purple), and the Phase II of FASER [303]. For reference, a band corresponding to the contribution light neutrino masses between 20 and 200 meV in a single generation see-saw type I model is shown (yellow). Larger values of the mixing angles are allowed if an extension to see-saw models is invoked, for instance, in an inverse or extended see-saw scheme

10.3 Heavy neutral leptons

The high intensity of the LBNF neutrino beam and the production of charm mesons in the beam enables DUNE to search for a wide variety of lightweight long-lived, exotic particles, by looking for topologies of rare event interactions and decays in the fiducial volume of the DUNE ND. These particles include weakly interacting heavy neutral leptons (HNLs) as right-handed partners of the active neutrinos, light super-symmetric particles, or vector, scalar, and/or axion portals to a Hidden Sector containing new interactions and new particles. Assuming the heavy neutral leptons are the lighter particles of their hidden sector, they will only decay into SM particles. The parameter space that can be explored by the DUNE ND extends into the cosmologically relevant region, and will be complementary to the LHC heavier mass searches.

Thanks to small mixing angles, the particles can be stable enough to travel from the production in the proton target to the detector and decay inside the active region. It is worth noting that, differently from a light neutrino beam, an HNL beam is not polarized, due to the large mass of the HNLs. The correct description of the helicity components in the beam is important for predicting the angular distributions of HNL decays, as they might depend on the initial helicity state. More specifically, there is a different phenomenology if the decaying HNL is a Majorana or a Dirac fermion [286, 287]. Typical decay channels are two-body decays into a charged lepton and a pseudo-scalar meson, or a vector meson if the mass allows it, and three-body leptonic decays.

A recent study illustrates the potential sensitivity for HNL searches with the DUNE ND [287]. The sensitivity for HNL particles with masses in the range of 10 MeV to 2 GeV, from decays of mesons produced in the proton beam dump that produces the pions for the neutrino beam production, was studied. The production of \(D_s\) mesons gives access to high mass part of the HNL production. The dominant HNL decay modes to SM particles have been included, as well as the basic detector constraints, and dominant background processes have been considered.

The experimental signature for these decays is a decay-in-flight event with no interaction vertex, typical of neutrino–nucleon scattering, and a rather forward direction with respect to the beam. The main background to this search comes from SM neutrino–nucleon scattering events in which the hadronic activity at the vertex is below threshold. Charged-current quasi-elastic events with pion emission from resonances are background to the semi-leptonic decay channels, whereas misidentification of long pion tracks as muons can constitute a background to three-body leptonic decays. Neutral pions are often emitted in neutrino scattering events and can be a challenge for decays into a neutral meson or channels with electrons in the final state.

We report in Fig. 40 the physics reach of the DUNE ND in its current configuration without backgrounds for a Majorana and a Dirac HNL. The sensitivity was estimated assuming a total of \(1.32 \,\times \, 10^{22}\) POT, i.e., for a running scenario with 6 years with a 80 GeV proton beam of 1.2 MW, followed by six years of a beam with 2.4 MW, but using only the neutrino mode configuration, which corresponds to half of the total runtime. As a result, a search can be conducted for HNLs with masses up to 2 GeV in all flavor-mixing channels.

The results show that DUNE will have an improved sensitivity to small values of the mixing parameters \(|U_{\alpha N}|^2\), where \(\alpha =e,\,\mu ,\,\tau \), compared to the presently available experimental limits on mixing of HNLs with the three lepton flavors. At 90% CL sensitivity, DUNE can probe mixing parameters as low as \(10^{-9}{-}10^{-10}\) in the mass range of 300–500 MeV for mixing with the electron or muon neutrino flavors. In the region above 500 MeV the sensitivity is reduced to \(10^{-8}\) for eN mixing and \(10^{-7}\) for \(\mu N\) mixing. The \(\tau N\) mixing sensitivity is weaker but still covering a new unexplored regime. A large fraction of the covered parameter space for all neutrino flavors falls in the region that is relevant for explaining the baryon asymmetry in the universe.

Studies are ongoing with full detector simulations to validate these encouraging results.

10.4 Dark matter annihilation in the sun

DUNE’s large FD LArTPC modules provide an excellent setting to conduct searches for neutrinos arising from DM annihilation in the core of the sun. These would typically result in a high-energy neutrino signal almost always accompanied by a low-energy neutrino component, which has its origin in a hadronic cascade that develops in the dense solar medium and produces large numbers of light long-lived mesons, such as \(\pi ^+\) and \(K^+\) that then stop and decay at rest. The decay of each \(\pi ^+\) and \(K^+\) will produce monoenergetic neutrinos with an energy 30 or \({236} \, \text {MeV}\), respectively. The \(236 \, \hbox {MeV}\) flux can be measured with the DUNE FD, thanks to its excellent energy resolution, and importantly, will benefit from directional information. By selecting neutrinos arriving from the direction of the sun, large reduction in backgrounds can be achieved. This directional resolution for sub-GeV neutrinos will enable DUNE to be competitive with experiments with even larger fiducial masses, but less precise angular information, such as Hyper-K [304].

11 Conclusions and outlook

DUNE will be a powerful discovery tool for a variety of physics topics under very active exploration today, from the potential discovery of new particles beyond those predicted in the SM, to precision neutrino measurements that may uncover deviations from the present three-flavor mixing paradigm and unveil new interactions and symmetries. The ND alone will offer excellent opportunities to search for light DM and to measure rare processes such as neutrino trident interactions. Besides enabling its potential to place leading constraints on deviations from the three-flavor oscillation paradigm, such as light sterile neutrinos and non-standard interactions, DUNE’s massive high-resolution FD will probe the possible existence of baryon number violating processes and BDM. The flexibility of the LBNF beamline opens prospects for high-energy beam running, providing access to probing and measuring tau neutrino physics with unprecedented precision. Through the ample potential for BSM physics, DUNE offers an opportunity for strong collaboration between theorists and experimentalists and will provide significant opportunities for breakthrough discoveries in the coming decades.