Baseline and other effects for a sterile neutrino at DUNE

We analyze the sensitivity of the Deep Underground Neutrino Experiment (DUNE) to a sterile neutrino, combining information from both near and far detectors. We quantify often-neglected effects which may impact the event rate estimation in a 3+1 oscillation scenario. In particular, we find that taking into account the information on the neutrino production point, in contrast to assuming a pointlike neutrino source, affects DUNE's sterile exclusion reach. Visible differences remain after the inclusion of energy bin-to-bin uncorrelated systematics. Instead, implementing exact oscillation formulae for near detector events, including a two slab density profile, does not result in any visible change in the sensitivity.


Introduction
Sterile neutrinos are hypothetical neutral leptons that interact with other particles via gravity, missing all the gauge interactions of the Standard Model. They can only be indirectly detected through their mixing with standard weakly-interacting neutrinos.
The possible existence of light sterile neutrinos with a mass of order 1 eV (for recent reviews see [1][2][3]) has been strongly hinted at in the mid-1990s with the results of the LSND accelerator experiment [4,5], which found evidence in favor of short-baseline (SBL) ν µ →ν e oscillations (ν e appearance). Theν e excess in the data indicated an oscillation to a sterile neutrino with a mass-squared difference to the interacting ones of ∆m 2 41 ≳ 0.1 eV 2 . The MiniBooNE experiment [6], developed later with the main purpose of testing the LSND signal and initially operating in neutrino mode, reached results which could not confirm the effect but, on the other hand, appeared to hint at additional effects in the low-energy part of the spectrum. More recently, new data from MiniBooNE operating in neutrino and antineutrino modes are consistent in energy and magnitude with the excess of events reported by LSND and lead to a combined LSND/MiniBooNE best fit at ∆m 2 41 = 0.041 eV 2 and sin 2 2θ eµ = sin 2 2θ 14 sin 2 θ 24 = 0.96 [7]. This result is, however, not confirmed by the OPERA [8] accelerator experiment. Moreover, large values of the mixing angles between active and sterile are in strong tension with solar neutrino data [9][10][11], as well as with the results from disappearance experiments, as we shall soon discuss. 1 Owing to the unclear situation, a new set of experiments -the Short Baseline Neutrino program at Fermilab -was proposed some time ago and is being developed to check the LSND and MiniBooNEν e excesses [15,16].
Since the release of the initial MiniBooNE results, the light sterile neutrino issue remained dormant, until it was revived in 2011 with the discovery of the reactor antineutrino anomaly [17][18][19]. In particular, reassessedν e fluxes were found to imply larger detection rates with respect to what was measured in several SBL reactor neutrino experiments (ν e disappearance), indicating a deficit in the data. This deficit could be made consistent with oscillations generated by ∆m 2 41 ≳ 0.5 eV 2 [1]. Further attention was also given to the previously-known Gallium neutrino anomaly [20][21][22][23]. In this context, a deficit of neutrinos originating from the intense 51 Cr and 37 Ar neutrino sources was found in the GALLEX [24] and SAGE [25] detectors (ν e disappearance) during calibration. The Gallium anomaly can be explained through oscillations generated by a sterile with ∆m 2 41 ≳ 1 eV 2 [1]. While recent updates to the flux calculation may have possibly resolved the reactor anomaly [26][27][28], the case for a Gallium anomaly has been strengthened by the latest results from BEST [29]. In a 3+1 interpretation, the combination of all Gallium data yields the best-fit point ∆m 2 41 = 1.25 eV 2 and sin 2 θ 14 = 9.4 × 10 −2 (see also [30]).
In order to investigate whether neutrino oscillations are the source of the aforementioned anomalies, one may measure the antineutrino flux from the reactor core at different distances, thus obtaining ratios of event rates. In this way, flux normalization and energy-shape uncertainties could be reduced. Along these lines, interesting results have been obtained by the NEOS experiment [31]. Their detector is, however, located at a fixed distance from the reactor and their data were compared with those from a similar reactor at Daya Bay at a much larger distance. This was followed by the DANSS experiment [32], measuring the flux from one single reactor core at different locations. Initially, there appeared to be a remarkable overlap of allowed regions of the two sets of data, indicating [33] ∆m 2 41 ≃ 1.3 eV 2 , sin 2 θ 14 ≃ 0.01 . (1.1) More recent results from DANSS have, however, exhibited a smaller significance [34,35] (see also [36]). The point of eq. (1.1) is not only in disagreement with the LSND and MiniBooNE preferred regions but also with the BEST joint Gallium fit. Other experiments [37][38][39][40] are also dedicated to measuring theν e reactor flux as a function of a varying (short) baseline. These, however, expect much smaller inverse-beta-decay event rates and signal-to-background ratios with respect to NEOS and DANSS (see, e.g., [41,42]). Of these experiments, Neutrino-4 has found an unexpected indication of oscillations corresponding to ∆m 2 41 ≃ 7 eV 2 and sin 2 θ 14 ≃ 0.1 [39], which does not overlap with the regions preferred by reactor data. The Neutrino-4 result has been reexamined [43], and its statistical significance remains a source of controversy.
If light sterile neutrinos exist and mix with the active ones, then they are expected to provide a signature in ν µ andν µ disappearance. So far, no such effect has been observed. A stringent upper bound on the corresponding 3+1 mixing angle θ 24 was obtained by combining MINOS and MINOS+ data [44], for ∆m 2 41 ≳ 10 −2 eV 2 . The preferred region for the sterile neutrino parameters is still unclear. There is, in fact, a strong appearance-disappearance tension that can be quantified through the violation of the approximate relation where sin 2 2θ eµ governs (-) ν e appearance, while sin 2 2θ ee = sin 2 2θ 14 and sin 2 2θ µµ ≃ sin 2 2θ 24 (for small angles) control ν e experiments which include LSND [4], KARMEN [45], OPERA [8], NOMAD [46], BNL-E776 [47] and ICARUS [48], namely [1]  Although these are approximate values, they are sufficient to illustrate the existing tension between appearance and disappearance data. Altogether, the question of the possible existence of sterile neutrinos and their intrinsic parameters is therefore far from settled. To this end, the upcoming Deep Underground Neutrino Experiment (DUNE) [49][50][51][52] may play an important clarifying role. In the DUNE setup, whose nominal mission is to perform precise measurements of neutrino properties and oscillations, a neutrino beam will be obtained from meson decays within a 194 m long decay pipe, which is followed by a hadron absorber. Originating from a proton beam impinging on a graphite target, positive and negative mesons are to be focused by a three-horn system [53], which can be operated either in forward (FHC) or reverse (RHC) horn current mode, leading to the production of a beam of mostly neutrinos or antineutrinos, respectively. Downstream from the graphite target, at a distance of approximately 574 m, a liquid-argon detector is to be installed as part of the near detector (ND) complex [54]. Further, a second liquid-argon far detector (FD) with a fiducial mass of around 40 kt is to be located at the Sanford Underground Research Facility in South Dakota, at a distance of 1297 km from Fermilab.
The purpose of this paper is to analyze the prospects of DUNE in constraining the sterile-neutrino parameter space, quantifying effects which may impact event rate estimations. We focus on the range ∆m 2 41 ∈ [0.1, 100] eV 2 for the sterile mass-squared difference. For this range of the new mass-squared difference and depending on the (anti)neutrino energy, the associated oscillation length may be of the order of the decay pipe length and of the distance between the graphite target and the ND complex. Thus, the location of the neutrino production point within the decay pipe can play an important role in determining event rates at the ND, as it influences the oscillation probabilities by affecting the oscillation baseline. This smeared-source effect goes beyond taking into account geometrical effects in the computation of unoscillated fluxes at the ND (present also in the 3+0 case) and has not been quantitatively assessed in previous works. 2 Note that the prospects for detecting light sterile neutrinos in the DUNE ND have been examined in refs. [55][56][57][58][59][60][61]. 3 Our analysis improves on the existing literature as it takes into account ND and FD event rates, shape uncertainties, exact 3+1 oscillation formulae in matter and, crucially, the aforementioned baseline dependence, while making use of the latest DUNE configurations [53].
The paper is organized as follows. In section 2 we present the theoretical framework to be used. In section 3 we expound the details of the DUNE simulation and present our main results. These are represented by discrepancies between exclusion curves in the sterile parameter space, which primarily arise as one takes into account the effect of a spatially-smeared source. It is found that such discrepancies persist after taking into account energy bin-to-bin uncorrelated systematic uncertainties. In section 4 we summarize our conclusions.

Neutrino oscillations in vacuum and matter
We work within the 3+1 framework, in which the 3ν-flavor paradigm is extended by a single sterile neutrino species ν s . In this framework, the propagation of a neutrino in the Earth's matter is described by the Hamiltonian (α, β = e, µ, τ, s) where U is the vacuum Pontecorvo-Maki-Nakagawa-Sakata (PMNS) lepton mixing matrix, ∆m 2 ij ≡ m 2 i − m 2 j are vacuum mass-squared differences, and E denotes the energy of the neutrino. The coherent forward elastic charged current (CC) scattering contribution A CC = 2 √ 2 G F N e E depends on the density N e of electrons in the medium. The net effect of neutral current (NC) scatterings is encoded in In the Earth, only the neutron density N n ≃ N e is relevant as the NC potentials of protons and electrons cancel. One has To study the propagation of antineutrinos in matter, it is sufficient to replace U → U * , A CC → −A CC and A NC → −A NC in eq. (2.1), see, e.g., [67]. The vacuum Hamiltonian is recovered when ρ = 0, 3) The eigenvalues of H vac are proportional to the vacuum mass-squared differences and read ∆m 2 i1 /2E. The matter Hamiltonian of eq. (2.1) can be brought to the form of eq. (2.3), up to the constant term∆m 2 11 /2E, which, while non-zero in general, does not impact oscillation probabilities. We denote the eigenvalues of the matter Hamiltonian by∆m 2 i1 /2E (i = 1, ..., 4), with∆m 2 i1 = ∆m 2 i1 +∆m 2 11 (see also appendix A). The oscillation parameters in matter, denoted with tildes, are defined by eq. (2.4). Neutrino oscillations in matter are sensitive to the matter mass-squared differences ∆m 2 ij ≡ ∆m 2 i1 − ∆m 2 j1 . Both vacuum and matter PMNS mixing matrices U andŨ can be parametrized in terms of six mixing angles and three CP violation (CPV) phases, 4 and eqs. (2.5) and (2.6) reduce to the standard parametrization [71] in the limits of vanishing sterile mixing angles, θ i4 = 0 andθ i4 = 0 (i = 1, 2, 3). The probability P αβ of transition between neutrino flavors α and β (or of survival for a given flavor, α = β) after traversing a length L under the influence of a Hamiltonian H is given by For propagation in vacuum, the eigenstates |ν k ⟩ of H = H vac are related to the flavor states via |ν α ⟩ = U * αk |ν k ⟩. By inserting the completeness relation k |ν k ⟩⟨ν k | = 1 in eq. (2.8), one obtains the known formula for the vacuum oscillation probabilities, where i, j = 1, 2, 3, 4, and the upper (lower) sign choice refers to (anti)neutrinos. This follows since, for antineutrinos, |ν α ⟩ = U αk |ν k ⟩. One has here defined (2.10) Instead, for propagation in a medium of constant density (ρ ̸ = 0), the eigenstates |ν k ⟩ of H = H mat obey |ν α ⟩ =Ũ * αk |ν k ⟩. The probabilities P mat αβ are then given by eq. (2.9) with the straightforward replacements U →Ũ and ∆ ij →∆ ij ≡ ∆m 2 ij L/4E. We are interested in the DUNE scenario, with an emphasis on ND events. The relevant setup is shown in Figure 1. One sees that the distance traveled by the neutrinos from the production point to the ND consists of a first segment of (event-dependent) length L 1 < 230 m in a low-density medium (ρ ≃ 0), which includes the 194 m decay pipe, and a second segment of fixed length L 2 ≃ 344 m in the Earth's crust with an approximately constant density, ρ ≃ 2.6 g cm −3 . In a two slab description of the matter It is clear that oscillation probabilities depend on the point of production of the (active) α-flavor neutrino, located at a distance L = L 1 + L 2 from the ND. As discussed in the introduction, our main goal is to quantitatively assess the impact of this dependence in the DUNE setup.
Event rates at the FD, on the other hand, are expected to be sensitive only to the average matter density [72], taken here to be ρ ≃ 2.6 g cm −3 . Furthermore, given its large distance to the FD, the source can be approximated as pointlike for the computation of FD events.
For both ND and FD event rates, a low-pass filter is applied at the probability level to appropriately average out unresolvable fast oscillations produced by large sterile mass-squared differences. Averaging the oscillation probabilities over a Gaussian energy distribution with standard deviation σ E (see also section 7.6 of [73]), we find 5 where the rightmost exponential factor is responsible for phasing out the high frequencies.
Eq. (2.12) is appropriate for the simulation of ND events, while for the FD, one has in the limit of a single baseline L in matter. Note that this filter is only considered due to the General Long Baseline Experiment Simulator (GLoBES) setup and sampling constraints. It should not overshadow the detector energy resolution (see section 3.2 for details on the chosen value). In a different setup, it could be exchanged by the increase in the precision of energy integrals -although with no expected improvement in the accuracy of the results, given the limitation from the detector resolution.
We further employ analytical expressions for the matter-dependent quantities ∆m 2 ij andŨ αiŨ * βi , making use of the exact diagonalization of the matter Hamiltonian performed in [75,76]. 6 The relevant results are given in appendix A, in eqs. (A.4), (A.5), (A.9) and (A.10). We have implemented the corresponding analytic probability engine within the GLoBES simulation software [78,79]. The relevant code is provided in the ancillary files exact 4nu.c and exact 4nu.h [80].
Accounting for matter effects is crucial for an accurate estimation of FD event rates. For ND event rates, however, the level of detail contained in eq. (2.12) turns out to be excessive in practice. As we numerically verify in what follows, it will be enough to consider the approximate result, (2.14) obtained in the limit of vanishing matter density and vanishing standard-neutrino masssquared differences (see appendix B for a derivation), where L = L 1 + L 2 is the eventdependent baseline for sterile oscillations.
In the next section, we make use of eqs. (2.13) and (2.14) to assess the sterile exclusion reach of DUNE. We contrast the usual scenario where the event dependence of the baseline is neglected, i.e., where P SBL αβ depends on a fixed baseline L = 574 m, with the more accurate scenario where information on the neutrino production point is incorporated in the event rate estimation.

Simulation details
The collisions of the primary proton beam on the graphite target and the resulting neutrino beam downstream from the hadron absorber have been simulated using version 6 Approximate results based on these works have been derived in [77]. v3r5p7 of the neutrino beamline simulator g4lbne [81], built against version 4.10.3.p03 of Geant4 [82,83]. We have assumed a 1.2 MW, 120 GeV proton beam (1.1×10 21 POT/y) and a cylindrical target, 1.5 m long and 1.6 cm in diameter, followed by a system of three magnetic focusing horns, operated at a current of ±300 kA for the FHC and RHC modes, respectively.
The g4lbne code allows one to compute unoscillated fluxes at the ND and the FD, given the positions of the detectors. The flux files are passed to the GLoBES software, which expects as input a single baseline L for each detector (or "experiment," in GLoBES terminology), as it assumes the source to be pointlike. Since the distance of the target to the ND is of the order of the decay volume length, the point-source approximation is not valid and source-volume geometry effects need to be taken into account -as they indeed are, even in 3+0 analyses. A routinely-used weighting procedure is enough to circumvent the GLoBES point-source restriction and encode these effects into a single ND flux file. This file, however, contains only the total unoscillated flux for each energy bin. Without further action, the information on the production point of each neutrino is lost to the subsequent stages of the simulation and an accurate simulation of the 3+1 scenario is impossible.
In Figure 2 we show the distribution of meson decays along the decay volume in FHC mode (left panel), as well as the neutrino flux effectively reaching the ND as a function of the location of the parent meson decay, obtained following the weighting procedure (right panel). This information on the distance traveled by the neutrinos may only enter GLoBES via i) the flux files produced by g4lbne or ii) the baseline L, which is used by the probability engine. To incorporate into GLoBES the desired baseline dependence of oscillation probabilities for ND event rates in the 3+1 scenario, we conceptually divide the decay volume into sections and associate each to a different GLoBES point-source "experiment", with its own fixed (average) baseline L. Thus, information on the distance traveled by the neutrinos enters the simulation at the level of the probability engine and via a flux normalization factor. 7 The flux arriving at the ND from each section is then passed to GLoBES as the flux of an independent "experiment". Of course, these do not represent actual independent experiments or detectors, as there is no way of discriminating the neutrino section of origin at detection. In this setup, one must therefore be careful to appropriately modify internal GLoBES functions, especially χ 2 functions (see section 3.2), in order to sum over all the virtual sections when computing ND event rates.
In the GLoBES simulation, we have taken the ND and FD to correspond to fiducial liquid argon masses of 50 ton and 40 kt respectively, following [54]. We have assumed 3.5 years of operation in FHC mode, as well as 3.5 years for the RHC mode. The ND is taken to be a scaled-down version of the FD and GLoBES AEDL files have been adapted from those in ref. [53]. We make use of the efficiencies, the smearing matrices, and the data on CC and NC neutrino interactions in argon from [53]. The latter have been obtained using version 2.8.4 of the GENIE event generator [84]. Downstream from the target, 20 sections are considered, each 2.5 m deep, followed by 10 additional ones, each 18 m deep, for a total of 30 sections covering a distance of 230 m. In order to collect enough statistics for each section, we have simulated 10 7 proton events within g4lbne for each mode of operation, FHC and RHC.

Chi-squared analysis
We are interested in producing exclusion plots in the sin 2 θ 14 , ∆m 2 41 and sin 2 θ 24 , ∆m 2 41 planes, comparing the single L = 574 m baseline computation with the more realistic one where the decay volume is sectioned, as described above. Given DUNE's capabilities and following ref. [53], we consider four "rules" altogether, corresponding to CC events associated with 1. ν e +ν e appearance in FHC mode, 2. ν e +ν e appearance in RHC mode, 3. ν µ +ν µ disappearance in FHC mode, and 4. ν µ +ν µ disappearance in RHC mode.
A "rule," as defined in the GLoBES language (see [74]), encompasses a number of signal and background channels and their associated systematic uncertainties. Each rule provides a separate contribution to the χ 2 . Taken together, the rules constitute the final link between the event computation and the statistical analysis.
Recall that in FHC and RHC modes, positive and negative meson decays along the decay pipe are responsible for the main component of the neutrino or antineutrino flux  Contaminations in FHC mode are due to theν µ , ν e , andν e beam components, whereas in RHC they arise from the ν µ , ν e , andν e components. Note that the (-) ν µ contamination channels are included as signal in the disappearance analyses, which group neutrino and antineutrino events for either mode (in contrast with, e.g., [85]). Instead, the (-) ν e contamination channels are one of the major backgrounds to ν e +ν e appearance. The remaining background channels correspond to misidentifications. These arise from (-) ν µ mistakenly accepted as (-) ν e and from NC events wrongly classified as CC events. The latter may arise from any component of the flux, and, in particular, from the main ν µ andν µ components. 8 If sterile neutrinos exist, then ν e andν e can appear at the ND as signal due to sterile-induced oscillations. In the 3+0 case, ν e andν e events at the ND can only originate from contaminations or misidentifications.
For illustrative purposes, we show in Figure 3 the distribution of neutrino and antineutrino events at the ND, after the assumed seven years of operation. We compare the 3+0 case and the 3+1 case for a set of sterile parameters to which DUNE will be sensitive (cf. results in the following section). Differences between numbers of events are also displayed. The left-and right-hand panels correspond to ν e +ν e appearance and ν µ +ν µ disappearance, respectively. Whereas  Figure 3. The lower limit on E guarantees the validity of the Gaussian averaging leading to eq. (2.12) (see also appendix C of [60]). It can be checked that this value of σ E is below the energy resolution of the detectors, see e.g. Figure 1 of ref. [86].
Once the binned event rates are obtained for each rule at both ND and FD, we proceed with the statistical analysis based on the minimization of the total χ 2 function, defined as The minimization is carried out over the parameter test values ω (see below) and over the normalization and shape nuisance parameters, ζ and ζ ′ . The index k goes over the set of 26 independent normalization systematic errors. These are obtained from ref. [55], together with the associated prior uncertainties σ k and are collected in Table 1. Energy bin-to-bin uncorrelated systematic errors (see also Table 1), which we consider at a later stage, are labeled by the indices r and i. These indices refer to the four rules mentioned earlier and to the 60 energy bins, respectively. A common 5% uncertainty (σ ′ = 0.05) is used for this set of 240 shape systematics. The first term in eq. (3.1) compares the event rates T d r,i (ω, ζ, ζ ′ ) obtained in the "test" 3+1 scenario with the event rates O d r,i (ω 0 ) obtained in the "true" or "observed" 3+0 scenario. It reads Here, the index d refers to the detector (ND or FD). It is essential to sum the contributions of all the different decay volume sections (interpreted as GLoBES "experiments") in computing T ND r,i and O ND r,i , before these quantities are introduced into χ 2 stat . This is accomplished by a nontrivial modification of GLoBES internal functions, mainly at the level of the source files glb minimize.c and glb sys.c, so that one effectively works with two experiments/detectors and the two corresponding nonzero matter densities.
The 3+1 event rates T d r,i depend on the test values ω and additionally on the relevant nuisance parameters ζ and ζ ′ . The vector ω contains the three mass-squared differences,   the six mixing angles, the three CPV phases, and the nonzero matter densities, ρ ND and ρ FD . The value of the nonstandard mass-squared difference ∆m 2 41 is fixed during the χ 2 minimization and, depending on the analysis, we also fix either θ 14 or θ 24 in order to obtain exclusion plots in the corresponding planes. For simplicity, we focus on a neutrino spectrum with normal ordering and also fix δ 13 to a best-fit value δ b.f. 13 = 1.28 π [87]. We marginalize over the remaining 11 parameters contained in ω, i.e. they are free to vary in the minimization. The 3+0 event rates O d r,i are evaluated using the parameter vector ω 0 , containing the central values of the standard oscillation parameters and densities, with sterile mixing angles being set to zero. We have Here, the N d r,c,i are the event rates for a given detector d, rule r, channel c, and energy bin i. The sum in c sequentially goes over the relevant signal and background channels for each rule, while the index l goes over all normalization systematics ζ contributing to a given channel. Nuisance parameters are correlated across rules, channels, and detectors, as indicated in Table 1. 10 To avoid a proliferation of shape systematics, these are introduced only for the expected dominant component of each rule, namely for the background events in (-) ν e appearance rules and for signal events in the (-) ν µ disappearance rules. Thus, nonzero shape systematics ζ ′ are effectively labeled only by the indices r and i, as previously indicated.
The second term in eq. (3.1) is the so-called prior term and includes information on measured neutrino parameters and their uncertainties, obtained from ref. [87], as well as a possible 5% uncertainty on the Earth density parameters. It reads where the sum extends only over the seven measured parameters ω j with central values (ω 0 ) j and standard deviations σ(ω j ), i.e., over two mass-squared differences, three mixing angles and two average densities. Finally, the last terms in eq. (3.1) correspond to penalty terms, appropriately controlling the magnitude of each nuisance parameter.

Results
In this section we present the results of minimizing the χ 2 function of eq. (3.1) independently in the sin 2 θ 14 , ∆m 2 41 and sin 2 θ 24 , ∆m 2 41 planes. At a first stage, we switch off shape systematics (ζ ′ = 0) and vary the 11 + 26 = 37 free parameters contained in ω and ζ. Our results are shown in Figure 4, where the usual computation of test events using a single baseline L = 574 m at the oscillation probability level (dashed red lines) is compared to the more realistic one, with a different L for each section of the decay volume (solid green lines). For each plane, the exclusion curves correspond to the locus of all points for which χ 2 min = 4.61 (90% C.L. for 2 d.o.f.). It is seen that for most of the values of the sterile mass-squared difference ∆m 2 41 , the inclusion of the information on the neutrino production zone decreases the sensitivity of DUNE to the sterile mixing angles θ 14 and θ 24 . This is manifested in the fact that the green curve in both graphs is distorted and displaced to the right relative to the red one. This effect is present in the case of θ 14 for ∆m 2 41 ≳ 0.3 eV 2 , while it is more pronounced for 0.1 ≲ ∆m 2 41 ≲ 4 eV 2 in the case of θ 24 .
We have repeated the minimization procedure including the bin-to-bin uncorrelated systematic errors ζ ′ , which implies marginalizing over 240 additional parameters (for a total of 277). The result is shown in Figure 5. There is a dramatic loss of sensitivity to sterile neutrino mixing parameters for both curves when the ζ ′ are switched on, as one may expect (cf. ND-only analyses [54,60] including shape systematics). However, visible differences persist between the point-source (red) and smeared-source (green)   exclusion curves even in this case. They indicate a decreased sensitivity that arises for ∆m 2 41 ≳ 1 eV 2 in the θ 14 plane and mainly for ∆m 2 41 ∼ few eV 2 in the θ 24 plane. The possible relevance of taking into account the neutrino production point in the estimation of 3+1 ND events has previously been recognized in ref. [55]. Given that oscillation probabilities depend on the baseline via the ratio L/E, a heuristic additional energy smearing has been used in [55] to mimic the smeared-source effect (while keeping a single baseline), in place of a more precise computation like the one considered here. In particular, the GLoBES migration matrices for the DUNE ND have been multiplied by an extra smearing matrix, obtained from the integration of a Gaussian with a 20% standard deviation in energy. Reproducing this procedure, we compute and present the corresponding 90% C.L. exclusion curves for the case ζ ′ = 0 in Figure 6 (blue lines), comparing them with the ones obtained previously. From this figure, one can see that the heuristic procedure does not, in general, provide a valid approximation to the more realistic case. Overall, for both sterile mixing angles, the correction seems to be underestimated for ∆m 2 41 ≲ 3 eV 2 and overestimated for ∆m 2 41 ≳ 10 eV 2 . As mentioned in section 2, one may use the approximate result of eq. (2.14) to assess the sterile exclusion potential of DUNE. We have numerically checked that this formula is sufficient for our purposes, i.e. that taking into account matter effects at the ND using exact expressions [76] and a two slab matter profile (with a nonzero density ρ ND ) via eq. (2.12) does not produce noticeable differences with respect to the considered limit of vanishing matter density and vanishing standard-neutrino mass-squared differences. This is evidenced by the overlap of the solid and dotted green lines in Figure 6.   To conclude this section, we comment on the behavior of the exclusion curves in the limits of small and large ∆m 2 41 . In both these limits, the probabilities ⟨P SBL αβ ⟩ become independent of L, the effects of a smeared source disappear at the probability level, and the curves in figs. 4 and 5 are expected to approach each other. In fact, for sufficiently large ∆m 2 41 , oscillations should rapidly average out due to the presence of the low-pass filter exponential factor (see eq. (2.14)). On the other hand, in the limit of small ∆m 2 41 , the sterile oscillation length becomes much larger than the distance traveled by the neutrinos from their production point to the ND. Therefore, while sterile effects may be detected at the FD, no effect will be present at the ND, and event rates become independent of the point of neutrino production within the decay volume. Whereas the curves do approach each other for small sterile mass-squared differences, it is not so in the large ∆m 2 41 limit since a large energy E in the low-pass filter exponential can counteract the largeness of ∆m 2 41 for the considered energy window and range of sterile masses. Note finally that the analyses presented here involve both detectors (ND and FD) and are thus able to provide bounds on sterile mixing for small values of ∆m 2 41 (see also [55]), unlike ND-only analyses [54,60]. Indeed, we find that the obtained exclusion curves for light sterile masses (∆m 2 41 ≲ 0.8 eV 2 ) are driven by FD event rates, which dominate the χ 2 in this regime. As the sterile oscillation lengths approach from above the scale of ND baselines, the sensitivity reach becomes fully attributable to ND event rates (observed for ∆m 2 41 ≳ 0.8 eV 2 ). Importantly, the differences observed between exclusion curves with the pointlike source approximation and without it (smeared source) remain even in the limit of no FD contribution to the χ 2 .    fig. 4 (no shape systematics considered), with additional exclusion curves superimposed. The solid blue curves are obtained following the heuristic procedure described in ref. [55]. Dotted dark green curves are obtained using eq. (2.12), which takes into account matter effects at the ND, and overlap with those obtained using the approximate SBL formula of eq. (2.14) (solid green).

Summary and Conclusions
The possible existence of light sterile neutrinos has been hinted at by several experiments. The upcoming Deep Underground Neutrino Experiment is expected to play a clarifying role. The DUNE near detector, whose main purpose is to measure the unoscillated neutrino flux in a 3+0 scenario, can play a crucial role in constraining the 3+1 parameter space since the presence of a sterile neutrino may induce short-baseline oscillations of the active neutrinos. These oscillations depend on the baseline effectively traveled by the neutrino, which can vary considerably between events, given the length of the decay volume when compared to the distance from the target to the near detector.
In this work, we have analyzed the sterile exclusion potential of DUNE, taking into account both near and far detector event rates and shape uncertainties. Unlike previous studies, which consider a single sterile oscillation baseline when computing oscillation probabilities, we have sought to quantify the effect of the event-dependent baseline. We find that taking into account such information on the neutrino production point, in contrast to assuming a pointlike source, affects DUNE's sterile exclusion reach. In most of the parameter space, the inclusion of this smeared-source effect leads to a decrease in sensitivity (see fig. 4). Furthermore, this effect is seen not to be appropriately modeled by an additional 20% Gaussian energy smearing (see fig. 6). Finally, we have verified that differences between the pointlike source and smeared-source computations persist even if one takes into account energy bin-to-bin uncorrelated systematic errors (see fig. 5), for which we have assumed a common 5% standard deviation.
Throughout our study, we have used the exact formulae of ref. [76] (see appendix A) for the evaluation of the matter oscillation probabilities affecting far detector event rates. As for the near detector, we have verified that the approximate short-baseline formula of eq. (2.14) and the exact one of eq. (2.12), which includes a two slab density profile, in practice, lead to the same results (see fig. 6).

B Approximate short-baseline probabilities
In this appendix, we explicitly derive the approximate oscillation probability formula of eq. (2.14), used for the estimation of ND event rates. We start by neglecting matter effects in short-baseline oscillations. Then, for a single baseline L in vacuum, which matches eq. (2.13) after dropping tildes, one has with the replacement U → U * for antineutrinos. Here, l.p.f. denotes the low-pass filter factor, Neglecting the solar and atmospheric mass-squared differences with respect to the sterile one and taking ∆m 2 41 ≃ ∆m 2 42 ≃ ∆m 2 43 , the above probabilities can be further approximated as l.p.f. 41 .