Improving CP Measurement with THEIA and Muon Decay at Rest

We explore the possibility of using the recently proposed THEIA detector to measure the $\bar \nu_\mu \rightarrow \bar \nu_e$ oscillation with neutrinos from a muon decay at rest ($\mu$DAR) source to improve the leptonic CP phase measurement. Due to its intrinsic low-energy beam, this $\mu$THEIA configuration ($\mu$DAR neutrinos at THEIA) is only sensitive to the genuine leptonic CP phase $\delta_D$ and not contaminated by the matter effect. With detailed study of neutrino energy reconstruction and backgrounds at the THEIA detector, we find that the combination with the high-energy DUNE can significantly reduce the CP uncertainty, especially around the maximal CP violation cases $\delta_D = \pm 90^\circ$. Both the $\mu$THEIA-25 with 17kt and $\mu$THEIA-100 with 70kt fiducial volumes are considered. For DUNE + $\mu$THEIA-100, the CP uncertainty can be better than $8^\circ$.


I. INTRODUCTION
The charge-parity (CP) symmetry violation is a key to understand the existence of baryon asymmetry in the Universe, namely, why there are more matter than antimatter [1][2][3][4][5]. There are at least two possible sources of CP violation in the Standard Model (SM) of particle physics: the CP phase in the quark mixing matrix [6,7] and the leptonic CP phases in the neutrino mixing matrix [8]. Especially, the leptonic CP phases at low energy play an important role [9,10] in the leptogenesis mechanism [11][12][13]. Both the Dirac and Majorana CP phases can contribute to the leptogenesis mechanism. However, only the Dirac CP phase manifests itself in neutrino oscillation and can be measured by oscillation experiments [14].
The nonzero reactor mixing angle (θ r ≡ θ 13 ) measured by Daya Bay [15] and RENO [16] heralds the precision era of neutrino oscillation experiments. A nonzero θ r allows the Dirac CP phase δ D to have physical effect since these two variables always appear together as sin θ r e ±iδ D in the standard parametrization [8] of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [17,18]. Typically, the neutrino oscillations from the muon flavor to the electron flavor (ν µ → ν e andν µ →ν e ) are used by the longbaseline accelerator experiments to measure δ D [19].
The current long-baseline experiments T2K [20] and NOνA [21] are approaching the discovery threshold. The 2019 T2K result [22] [23] for NO with best-fit value at vanishing CP phase, δ D = 0 • . In 2021, T2K and NOνA updated their results with the best-fit value from T2K remaining the same, δ D = 252 •+40 • −33 • (NO) and 281 •+28 • −31 • (IO) [24] while the NOνA best fit changes to δ D = 148 •+49 • −157 • for NO [25]. Although not significant, there is a tension between the T2K and NOνA data that they exclude each other at 1σ C.L. [26,27]. New physics can explain the tension. Both belonging to accelerator neutrino experiments, T2K and NOνA have very different configurations. While the T2K baseline is 295 km and the peak energy is at 0.6 GeV, the NOνA baseline is 810 km and peak energy at 2 GeV. These differences in baseline and beam energy leave room for new physics. For example, the non-standard interaction (NSI) contributes extra matter potential [28] and hence its effect on oscillation probabilities is energy dependent [29,30] to provide a possible solution [31,32]. The tension is reduced when the data are analyzed in the context of Lorentz invariance violation (LIV). However, it is accompanied by a new mild tension between the best-fit values of sin 2 θ 23 [33]. Besides, the non-unitarity mixing due to heavy neutrinos [34][35][36] allows extra CP phases to fake the genuine CP effect which can also explain the tension [37]. However, a more recent work [38] points out the non-unitarity cannot explain the tension with the bounds on non-unitarity parameters from the combination of short-and long-baseline data. Similar thing happens for the light sterile neutrino scheme [39]. Whether this tension is truly new physics or not needs further investigation at current and future experiments.
Even if the current tension between T2K and NOνA measurements vanishes with more data, correct interpretation of CP measurements still faces intrinsic issues including event rate inefficiency, δ D ↔ π − δ D degeneracy, and large CP uncertainty around the maximal values [40][41][42]. These issues still remain for the next-generation experiments like T2HK [43] and DUNE [44]. Although its wide spectrum can help to reduce the δ D degeneracy, DUNE has much larger matter effect than T2K and NOνA due to higher energy peaking around 2.5 GeV [45]. At long-baseline experiments, the genuine CP effect can be faked by the ubiquitous matter effect, reducing the experimental sensitivity to δ D [46][47][48]. In addition, the uncertainties in the matter effect can also reduce the CP sensitivity at DUNE [45].
In this paper, we propose µTHEIA as combination of the THEIA detector [49][50][51][52] and a µDAR neutrino flux to improve the Dirac leptonic CP phase measurement together with DUNE. Sec. II summarizes the contamination of matter effect in the CP measurement and explains why the µDAR neutrino flux with lower energy can help. Then in Sec. III, we describe the low-energy mode at µTHEIA and the high-energy mode at DUNE, including selection criteria, energy reconstruction, smearing, and backgrounds. A combination of µTHEIA and DUNE can significantly improve the CP sensitivity as illustrated in Sec. IV. Therein, we also give the details of simulation and χ 2 analysis. Our µTHEIA proposal is compared with the existing configurations/proposals in Sec. IV D and summarized in Sec. V.

II. IMPROVING CP MEASUREMENT WITH MULTIPLE BASELINES AND BEAM ENERGIES
The current T2K and NOνA experiment aims for the discovery of leptonic CP violation, namely, excluding δ D = 0 and π. Once a nontrivial δ D is measured, the next step is precision measurement of its value. Several experimental configurations have been proposed to improve the CP measurement after the reactor mixing angle θ r was measured by Daya Bay and RENO. The upgrade from existing experiments includes: (1) Intensity Upgrade: the beam intensity is significantly enhanced, such as T2K-II [53]; (2) Detector Upgrade: T2HK [43] has a much larger detector Hyper-K than Super-K; (3) Spectrum Upgrade: DUNE [44] adopts wider on-axis spectrum than the off-axis one of NOνA; and (4) Baseline Upgrade with longer baseline such as T2HKK [54][55][56] and DUNE. In addition, there are also several new proposals: (5) the accelerator experiments such as P2O [57], ESSνSB [58], and MOMENT [59]; (6) CP measurement with sub-GeV atmospheric neutrino oscillation such as Super-PINGU [60,61], Super-ORCA [62], and even at JUNO [63] or DUNE [64]. Comparison among various experimental configurations can be found in [65][66][67][68].

A. Matter Contamination on CP Measurement
However, the neutrino energy for all these designs is not low enough to avoid contamination from the ubiquitous matter effect [45,[69][70][71][72][73][74]. We will try to describe how the CP measurement is contaminated by the matter effect. Based on this, one can see the possible solutions.
The neutrino propagation through matter is described by the following Hamiltonian [28,75], The first term is the vacuum Hamiltonian that is a product of the PMNS matrix U [17,18] and the diagonal mass matrix with solar ∆m 2 s ≡ ∆m 2 21 and atmospheric ∆m 2 a ≡ ∆m 2 31 mass squared differences. The neutrino energy E ν in the denominator contributes as an overall factor. So it is not a suppression in the vacuum term but actually an enhancement of the matter potential in the second term.
Induced by the SM weak interaction, the matter potential V ≡ √ 2G F n e is proportional to the Fermi constant G F and the electron number density n e , and only appears in the first element of the potential matrix in (1) for the electron flavor [28]. The matter potential V (x) is not just proportional to the matter density ρ(x) but also the number of electrons per nucleon Y e (x). Both the matter density and chemical composition vary with position. Consequently, the matter potential is in general also a function of x. In the Earth crust, the electron fraction is approximately Y e ∼ 0.5 and the average matter density isρ ≡ ρ(x) ≈ 2.845 g/cm 3 [45]. For simplicity, we ignore the density variation and adopt the averaged matter densityρ as constant along the baseline of DUNE.
To make the CP and matter effect explicit, we perform a series expansion of the ν µ → ν e (ν µ →ν e ) oscillation probability in terms of the ratio between the two mass squared differences α (≡ ∆m 2 s /∆m 2 a ≈ 3%) and the reactor mixing angle s r (≡ sin θ r ≈ 0.15) [76][77][78], The sign ± (∓) is for neutrino (the upper one) and anti-neutrino (the lower one), respectively. For convenience, we have used (s a,r , c a,r ) ≡ (sin θ a,r , cos θ a,r ) to denote the sine and cosine functions of the atmospheric (θ a ≡ θ 23 ) and reactor (θ r ) mixing angles while θ s ≡ θ 12 is the solar mixing angle. The matter term A ≡ 2E ν V /∆m 2 a appears in two combinations, sin(A∆ a )/A and sin[(1 ∓ A)∆ a ]/(1 ∓ A). In the limit of tiny matter effect, A → 0, the two combinations reduce to approximately ∆ a and sin ∆ a , respectively. In addition, the CP phase δ D appears as a linear combination with the atmospheric oscillation phase ∆ a ≡ |∆m 2 a |L/4E ν where L is the oscillation baseline [79].
For neutrino CP measurement, the essential observable is the difference between the neutrino and antineutrino oscillation probabilities, P νµ→νe − Pν µ→νe ∝ sin ∆ a sin δ D , that is proportional to sin δ D in the absence of matter potential. However, a realistic measurement has sign difference in not just the Dirac CP phase δ D but also the matter term A. With a typical size of A ≈ (0.05, 0.13, 0.21) estimated with the peak neutrino energies (0.55, 2, 2.5) GeV at T2K/T2HK, NOνA, and DUNE, respectively, the matter effect on neutrino CP measurement cannot be ignored. The matter potential can fake the genuine CP violation and blur the CP measurement.
At a single long-baseline neutrino oscillation experiment, there is only one independent CP observable but two parameters (A and δ D ). To disentangle the matter contamination (A) from the genuine CP effect (δ D ), a combination of two different baselines is a promising choice. For example, the T2HKK with a 295 km baseline to the Kamioka site and a much longer baseline around 1100 km to the Korea site can effectively remove the faked CP by matter potential and achieve a better CP uncertainty [66,67,80,81]. In addition, the atmospheric measurement with neutrinos produced around the Earth intrinsically has multiple baselines [60][61][62][63][64]. Nevertheless, the matter effect at accelerator and atmospheric neutrino oscillation experiments is not negligible in the first place. Even with multiple baselines, intrinsic uncertainty from the matter effect can still reduce the sensitivity of δ D .

B. Improvement with µDAR Neutrinos
A better way is significantly reducing the matter effect with low-energy neutrino beam to make A 1 [82]. One possibility is using the muon decay at rest [83], such as DEAδALUS [84], JUNO supplemented with µDAR sources [85,86], TNT2K/TNT2HK [73,[87][88][89][90], and C-ADS [91]. All these designs share the feature of multiple baselines. Especially, the TNT2K/TNT2HK configuration incorporates both low-and high-energy beams by supplementing the existing T2K/T2HK with µSK/µHK (µDAR source together with the Super-K/Hyper-K detectors) to significantly improve the CP sensitivity. The possibility of detecting µDAR neutrinos at long baseline has also been studied but is unfortunately diluted too much over such a long baseline [92].
As mentioned above, the effect of matter potential on the neutrino oscillation is modulated by the neutrino energy. The higher neutrino energy, the larger contamination on the CP measurement. Comparing with the 30% effect at T2K for E ν ≈ 550 MeV [73], DUNE with peak energy around 2.5 GeV suffers from larger matter effect. So it is more urgent for the DUNE experiment to have a complementary baseline with low-energy beam to further improve the CP measurement.
However, it is impossible to simply add a µDAR source and share the same liquid Argon detectors of DUNE in a similar way as TNT2K. This is because there are no free protons to provide inverse beta decay (IBD) for unique probe of the electron anti-neutrino and hence theν µ →ν e oscillation. Besides, theν e −Ar cross section is too small to detect the µDAR flux at DUNE [44].
A new THEIA detector at the same site of SURF was recently proposed [49][50][51][52]93]. With a new technique of water-based liquid scintillator (WbLS), it is possible to use both scintillation and Cherenkov lights [94][95][96][97][98]. This opens the possibility of detecting the low-energy µDAR neutrino oscillation to supplement the high-energy mode at DUNE. For convenience, we call the combination of µDAR and THEIA as µTHEIA.
The difference between the oscillation probabilities with and without matter, δP µe ≡ P µe (A) − P µe (A = 0), in Fig. 1 shows explicitly the matter effect at the DUNE and µTHEIA configurations. For DUNE (L = 1300 km with blue and yellow lines), the difference can be as large as δP µe ≈ 0.03 for E ν /L ≈ 1.5 MeV/km which is roughly 52% (62%) of the CP-violating oscillation probability P µe = 0.058 (|P ν µe − Pν µe | = 0.048). Even for the neutrino energy peak at E ν /L ≈ 1.9 MeV/km, the size of the matter effect is still as large as 0.025 and 33% (52%) of P µe = 0.075 (|P ν µe − Pν µe | = 0.048). It is interesting to see that the probability difference is almost independent of the Dirac CP phase δ D . This is because the major matter effect, δP µe ≈ 8s 2 r s 2 a A ≈ 0.02, comes from the second term of (3) at the oscillation peak without involving δ D . Although the matter effect also appears through the third term of (3), the effect is further suppressed by a prefactor of ∆ a sin 2θ s × (α/s r ) ≈ 0.3 and hence is a minor effect. Altogether, the matter effect at DUNE is at the same order as the genuine CP effect.
In contrast, the matter effect at µTHEIA (L = 38 km with red and green lines) is negligibly small. Being essentially insensitive to the matter potential, µTHEIA can focus on the genuine CP phase while DUNE probes both. Their combination can significantly improve the CP sensitivity. We will further discuss the details of µTHEIA and its interplay with DUNE in Sec. III for neutrino detection and Sec. IV for CP sensitivity.

III. NEUTRINO DETECTION AT THEIA AND DUNE DETECTORS
As proposed above, the essential feature of the DUNE and µTHEIA complex is a combination of different baselines and neutrino beams. In addition, the DUNE and THEIA detectors are also quite different with liquid Argon and WbLS targets, respectively. Each combination of neutrino beam and detector has its own characteristics in neutrino detection. For clarity, we elaborate separately the details of the "low-energy mode" (LEM) that the µDAR beam is detected by the THEIA detector in Sec. III A as well as the the "high-energy mode" (HEM) that the LBNF beam is detected by both the DUNE and THEIA detectors in Sec. III B. We study in detail the  [30,50] MeV while the DUNE energy window spans the whole range. event reconstruction for both signal and background to obtain their normalized transfer tables. The event rate including the information of flux, running time, cross section, and detector size can be found in the following Sec. IV when estimating the CP sensitivity.

A. The Low-Energy Mode
The µDAR neutrinos are produced by a cyclotron complex. For example, a typical 800 MeV proton beam hits a thick target to first generate pions. Although both π ± can be produced, π − is mostly absorbed by the positively charged nuclei while π + decays at rest via π + → µ + ν µ . The decay product µ + also loses its energy and decays at rest via µ + → e + + ν e +ν µ . During this process, three neutrinos (ν µ , ν e , andν µ ) are produced. Of them,ν µ experiences theν µ →ν e oscillation that is of interest to CP measurement. Since µ + decays at rest,ν µ has a wellunderstood spectrum with maximum energy of 53 MeV [99].

Signal at THEIA Detector
As mentioned in the previous section, the low-energȳ ν µ cannot be detected by the DUNE detector. But the THEIA detector is an ideal equipment. For the major oscillation channelν µ →ν e for CP measurement with µDAR neutrinos, the IBD process (ν e + p → e + + n) is an ideal detection method. With a large volume of WbLS, THEIA has a significant fraction of free protons (hydrogen) to allow IBD. Both positron and neutron in the final state are detectable by THEIA.
The WbLS allows THEIA to detect both scintillation and Cherenkov lights. There are several differences between the Cherenkov and scintillation lights that can be used for separation [94,97]: (1) the arrival time, i.e, Cherenkov light arrives nanoseconds earlier than the delayed scintillation light; (2) the angular topology, i.e, the Cherenkov light emission will cause a local enhancement on top of the isotropic scintillation signal; (3) the wavelength, i.e, the scintillation light has shorter wavelength while the Cherenkov light has relatively longer one.
For particle identification, the Cherenkov ring from e ± and γ has a smeared pattern. In contrast, the heavier muon or charged pion has much sharper Cherenkov ring. In addition, the production rates of Cherenkov and scintillation lights are different for different particles [94]. We will give more details once needed in later discussions.
Both Cherenkov and scintillation lights can be used for energy reconstruction. For low energy reactor and supernovaν e neutrinos, the energy is typically reconstructed with linear form, Eν = E e + m n − m p [100] where m n (m p , m e ) is the neutron (proton, electron) mass and E e the positron energy. Nevertheless, this linear formula cannot reconstruct the neutrino energy very precisely. To see the deviation clearly, we simulate the IBD events with GENIE [101,102] and reconstruct the neutrino energy according to the linear formula. As shown in Fig. 2, the reconstructed energy spectrum for neutrinos with E truē ν = 40 MeV spreads into a trapezoid (green solid line) and the central value shifts leftward by almost 2 MeV. In other words, for the µDAR neutrinos with O(10 MeV) energy, the neutrino energy reconstruction is no longer a linear dependence on the positron energy with a constant shift.
Another possibility is adding a correction term with nonlinear dependence [103], where ∆ ≡ m n − m p is the difference between the neutron (m n ) and proton (m p ) masses. The blue solid line clearly shows that the reconstructed trapezoid spectrum shifts back to center around the true neutrino energy at E truē ν = 40 MeV. Nevertheless, the reconstructed spectrum is still a trapezoid even without smearing due to the detector resolution. With a half width being almost 1.3 MeV, the corresponding energy uncertainty can be as large as 3.25% which is even larger than the statistical uncertainty as elaborated below.
To get a precise reconstruction of the neutrino energy, we adapt the reconstruction formula from long-baseline experiments [104] to incorporate the scattering angle, where p e is the electron momentum while θ e is the angle between the neutrino beam and positron direction. The directional information is a benefit due to the capability of detecting Cherenkov light for a WbLS detector. In the absence of detector resolution, (5) can obtain exactly the true neutrino energy and the reconstructed spectrum becomes a δ-function (red solid line) in Fig. 2. Note that the positron energy E e cannot be directly measured. The positron not just loses its kinetic energy but also annihilates with an environmental electron to produce two 511 keV γ's. The total deposited energy that is visible in the detector receives an extra electron mass in addition to the positron energy, E vis ≡ E e + m e . In other words, the positron energy in the linear as well as non-linear energy reconstruction formula (4) and (5) is reconstructed as E e = E vis − m e .
With systematical uncertainty from the reconstruction formula eliminated, the remaining uncertainty mainly comes from the detector resolutions of the visible energy E vis and scattering angle θ e . Since the light yield is proportional to the deposited energy, the visible energy E vis can be reconstructed from the number of photoelectrons (p.e.) which is roughly 80 (130) p.e./MeV for the Cherenkov (scintillation) light [97]. The relative statistical uncertainties are (σ Evis /E vis ) Ch = 1/ 80 × E vis /MeV for the Cherenkov light and (σ Evis /E vis ) sc = 1/ 130 × E vis /MeV for the scintillation light, respectively. The combined uncertainty is Since the µDAR neutrino energy is typically O(10 MeV), the relative error σ Evis /E vis 2% is actually quite good. Note that the scintillation light yield is more efficient than the Cherenkov one and both can play important roles in energy reconstruction. At E truē ν = 40 MeV, the positron energy resolution leads to about 1.17% uncertainty on the reconstructed neutrino energy.
In addition, the scattering angle θ e in (5) also receives an imperfect detector resolution and blur the reconstructed neutrino energy. Although the actual resolution for the µDAR flux that has higher energy should be better, we take a conservative value σ θ = 10 • [51] estimated for supernova neutrino. The angular uncertainty leads to a 0.44% uncertainty on the neutrino energy reconstruction which is slightly less than the one from positron energy resolution. In Fig. 2, the dashed lines show the combined neutrino energy uncertainty, which is 1.28% for E truē ν = 40 MeV assuming Gaussian smearing. Once produced, the neutron from IBD would experience thermalization and lose energy before being captured. The scintillation light emitted during this process can mix with the light from positron. Consequently, the visible energy receives an extra contribution from the neutron kinetic energy T n , E vis = E e + m e + Q F × T n where Q F is the neutron quenching factor. For the µDAR neutrino at 40 MeV, the neutron kinetic energy T n is typically 3 MeV according to GENIE simulation. With WbLS, it is possible to reconstruct the neutron kinetic energy T n by measuring the positron scattering angle [105]. The neutron quenching factor Q F keeps decreasing with T n and is 15% at T n = 0.25 MeV [105]. Then from the measured visible energy E vis , one can first solve the positron energy E e = E vis − m e − Q F × T n and put it into (5) to reconstruct Eν. In principle, the neutron kinetic energy is not a problem for the IBD energy reconstruction. But the neutron quenching factor Q F still remains to be measured in a WbLS as far as we know. So for simplicity, we omit the neutron kinetic energy in our current phenomenological study.

Background Suppression with WbLS
As illustrated in Sec. III A 1, theν e signal contains not just a positron but also a neutron in the final state. Although a single positron can already allow precise energy reconstruction, it receives various backgrounds. The final-state neutron is extremely important for selecting out the IBD signal with double coincidence.
The neutron is produced at the same time as the positron. But the neutron signal comes out much later with a delay of 250 µs [106]. Since the time resolution required to separate the Cherenkov and scintillation lights (typically 0.1 ns) is much smaller [51], the delayed neutron signal can be well separated from the positron one. When neutron is captured by nuclei, the delayed γs also produce e-like Cherenkov rings and scintillation lights to allow neutron tagging at THEIA [106]. The neutron tagging efficiency can reach almost 90% [107].
The double coincidence requiring both positron and neutron signals can remove those events that contain only one faked positron. But there are still various backgrounds that can survive [87,108] including the beam neutrinos from the complementary high-energy beam neutrinos, the low-energy reactor, solar and supernova neutrinos, the intrinsic beam neutrinos from the same µDAR flux, as well as the atmospheric neutrinos. The invisible muon decay from atmospheric neutrinos is particularly difficult to remove since both positron and neutron can be produced to fake the signal. Fortunately, the new technique of WbLS can be very efficient in suppressing these backgrounds as we elaborate below.
Reactor, Solar, Geo-, Supernova, and DUNE Beam Neutrinos -At low energy, there are four neutrino sources from reactor [109,110], solar [111], geo-radioactivity [112,113], and supernova [114,115]. While the solar neutrinos are not anti-neutrinos and hence cannot fake theν e IBD signal, the other three may be incorrectly identified as the signal. Fortunately, all these neutrinos typically have much lower energy than majority of the µDAR neutrinos. With an energy cut, Eν > 30 MeV, they can be safely removed [87]. On the other side, the LBNF beam can also contributeν e flux. With pulsed beam at much higher energy, a combined cut on the arrival time and energy, Eν < 55 MeV, will eliminate the contamination from the complementary DUNE experiment.
Intrinsic µDAR Beam -During the production of µDAR with cyclotron, various neutrinos can be produced via π ± → µ ± + ν µ (ν µ ) and µ ± → e ± +ν µ (ν µ ) + ν e (ν e ). One needs to examine all possible intrinsic backgrounds. Theν e flux can produce exactly the same IBD signal and becomes an intrinsic background from the same µDAR source. There is no way to remove or suppress this background on the detector side. Fortunately, the negatively charged π − is efficiently absorbed by the positively charged nuclei inside the target before decay. So the µ − production is much smaller than µ + with a suppression as large as 10 −4 [87]. Since theν µ →ν e oscillation probability is typically O(1%) and hence two orders larger, the µ − contamination can be safely neglected.
The second beam background is the electron neutrino ν e from the same µDAR. After oscillation, most ν e neutrinos can survive and scatter with oxygen (ν e + 16 O → e − + 16 F ) or carbon (ν e + 12 C → e − + 12 N ) 1 to produce an electron in the final state. Fortunately, the scattering cross section with oxygen is much smaller than the IBD one in the µDAR energy range [116,117] and hence can be neglected. Even if some events can still be produced, the double coincidence of neutron capture can remove this background.
The muon neutrino and anti-neutrino cannot experi-ence charged-current scattering with not enough energy to produce a µ ± in the final state. However, all neutrinos and anti-neutrinos can elastically scatter with electron to produce an energetic electron in the final state. But these backgrounds can also be removed by requiring neutron capture [94,97].
Atmospheric Neutrinos -The atmospheric neutrino flux [118,119] contains all neutrino flavors (ν e ,ν e , ν µ ,ν µ , ν τ , andν τ ) to contribute as background [87]. With much higher energy, the atmospheric neutrino backgrounds can experience all types of neutrino-nucleus scattering and hence need much more dedicated treatment. Fortunately, the excellent performance of the WbLS with both Cherenkov and scintillation detection can effective suppress these backgrounds. Exceptν e andν µ , the other components can be removed at the THEIA detector as we elaborate below. ν e : The charged-current quasi-elastic scattering (CC-QES) process ν e + 16 O → e − + 16 F cannot produce a neutron. Then the neutron tagging with WbLS can efficiently remove this background. For the atmospheric ν e with higher energies, 30% of the CC-QES scattering process can kick out a neutron and a proton from the nuclei via ν e + 16 O → e − + n + 15 O * + p. Since energy is needed to kick out n and p, the primary electron energy is much smaller than the neutrino energy. Mainly E ν ∈ [80,130] MeV will produce electrons with energy reconstructed in the µDAR range [30,55] MeV. In addition, roughly half of the single neutron events are also accompanied by a monoenergetic photon (∼6.2 MeV) due to the de-excitation of 15 O * , which can also serve as background veto. Both electron (T e ) and proton (T p ) kinetic energies deposit as visible energy E vis . Requiring E vis ∈ [30,55] MeV, only 0.07 events per year at THEIA-25 can survive, which is a negligible amount.
Moreover, those neutrinos with even higher energy can also have resonant (RES) and deep inelastic (DIS) scatterings with at least a pion (π 0 or π + ) starting around E ν ≈ 200 MeV, and possibly a neutron in the final state. In addition, the interaction can also produce protons in the final state or kick off protons from the nucleus [120]. The neutron capture process emits a single 2.2 MeV photon [97] that can be used to separate RES and DIS events. In the energy range E ν ∈ [200, 600] MeV that can contribute to the µDAR IBD energy window, roughly 30% of CC-RES can have a single neutron. The events containing π 0 can be vetoed by energy cut E vis < 55 MeV, since the π 0 decays immediately into a pair of photons each with energy E γ ≥ m π /2. This reduces the remaining atmospheric ν e background by another 50% to only 15% of the total CC-RES. The charged π + first deposits its kinetic energy as scintillation light and then decays at rest to produce a µ + , since its decay length is typically at least one order longer than the radiation length in 30 35 Atmos ν e @ THEIA-25 interactions as a function of the visible energy Evis at THEIA-25. With energy deposits from the primary electron (Te, red), charged π + (+Tπ, blue), µ + (+Tµ, green), the Michel electron (+T e + , yellow), and proton (+Tp, purple) gradually added up, the νe event spectrum moves out of the IBD energy window [30,55] MeV to higher energy. Similarly, theνe spectrum has the same feature with the primary electron (Te, red), charged π − (+Tπ, blue) or with Cherenkov cut (w/ Cher. cut, green), and proton (+Tp, purple). Multiplication factors have been used to make the highly suppressed spectrum visible.
water [121]. Then the energy deposition is composed of five parts: 1) the kinetic energy T e of the primary electron; 2) the proton kinetic energy MeV that is uniquely determined by the pion decay at rest; and 5) the e + energy T e + ranging from 0 to m µ /2. Since the positron always annihilation with an environmental electron, we also include the released 2m e photon energy in T e + for convenience. For a π + with energy E π , the minimal total energy deposit is E π − m π + T µ + 2m e ≈ E π − 134 MeV if both the primary and Michel electrons have vanishing momentum. Those events with E π 189 MeV can be vetoed by E vis ≥ 55 MeV. Requiring the primary and Michel electrons to be produced at rest is actually very stringent. Even with less energetic π + , the atmospheric ν e background can be easily vetoed. As shown in the left figure of Fig. 3, once the π + , µ + , and e + energies are considered, the ν e CC-RES event rate is suppressed to negligible level with just around 8 × 10 −4 events/year for E vis ∈ [30,55] MeV at the THEIA-25 detector. Fig. 3 also shows the CC-DIS case which has only 20% for single neutron events and its contribution is only around 10 −5 event per year. The typical value of the proton kinetic energy is around O(10 MeV) for sub-GeV neutrinos. It can further suppress the background by one order 1 Since the scintillator composition and fraction for THEIA are not decided yet and only a possible fraction 1% ∼ 10% is mentioned [51], we leave the carbon contribution open and focus on the major water target.
for RES and almost down to zero for DIS, shown as the purple curves of Fig. 3.
The electron anti-neutrinoν e that scatters with Hydrogen via CC-QES (or precisely IBD) process is an irreducible background to the IBD signal and there is no experimental solution. It contributes the major background for µTHEIA. Since the direction of the incomingν e is unknown, the scattering angle θ e in (5) cannot be correctly measured but only inferred from the µDAR source direction. This wrong scattering angle effect introduces significant smearing for the reconstructed neutrino energy Eν. With more details elaborated in App. A, Fig. 4 shows that the smearing can be as large as 3.4% ∼ 5.5% at half height for E ν ∈ [30,55] MeV. The wrong scattering angle effect is much larger than the detector resolution in Fig. 2 with known direction. With higher neutrino energy, the wrong scattering angle effect becomes more severe. The atmosphericν e CC-QES background is estimated as 1.1 event per year at THEIA-25.
Theν e CC-RES and CC-DIS events mainly consist of a π 0 or π − in addition to the primary positron. Slightly higher than the ν e mode, the single neutron events contribute around 50% (60%) of the total CC-RES (CC-DIS) events. The π 0 can be vetoed by the energetic photon with E γ ≥ m π /2. Since the pion decay length is much longer than its radiation length in water, π − first loses its kinetic energy and then is absorbed by the positively charged nuclei. The minimal energy deposit is then only its kinetic energy T π = E π − m π and no µ − or subsequent electron from µ − decay. For comparison with the ν e mode, the deposited E vis is shown in the right panel of Fig. 3. As mentioned earlier, 80 (130)  can be produced for each MeV energy deposit [97] as Cherenkov (scintillation) lights. So the required µDAR energy window [30,55] MeV corresponds to [2400, 4400] Chereknov photons and [3900,7150] scintillation ones, respectively. The atmosphericν e CC-RES (CC-DIS) background only contributes 1.2 × 10 −3 (1.5 × 10 −3 ) event per year at THEIA-25. This is negligibly small compared with the atmosphericν e CC-QES background with 1.1 events per year. This background is further suppressed to 1/4 for the RES and almost 0 for the DIS processes if the proton kinetic energy is considered.
ν µ /ν µ : For the atmospheric muon neutrinos ν µ /ν µ , the charged-current scattering produces a µ ∓ in the final state. If the muon energy is above the Cherenkov threshold, the µ-like ring has different pattern to be distinguished from the e-like one. In addition, the ratio between the Cherenkov (N ch ) and scintillation (N sc ) lights is also quite different between electron (N ch /N sc ≈ 0.2) and muon (N ch /N sc 0.08) rings [94].
On the other hand, the Cherenkov light alone cannot see an "invisible muon" below the Cherenkov threshold but only its decay product e ± [87]. Fortunately, with a muon lifetime of 2.2 µs, the time resolution (∼ 0.1 ns) of WbLS is good enough to separate the muon scintillation light from the e ± Cherenkov and scintillation lights [94]. In other words, the WbLS can identify the invisible muon with triple coincidence (µ scintillation light, e ± Cherenkov and scintillation lights, and the 2.2 MeV delayed γ from neutron capture). The triple coincidence can essentially remove all the invisible muon background. Even for a conservative study by requiring the Cherenkov photon and the combined scintillation photon to satisfy N ch /N sc > 0.4, the background rate is suppressed to only 2% [94].
ν τ /ν τ : Although tau neutrinos also exist in the atmospheric flux, the primary τ ± from the CC scattering can decay into either electron or muon with roughly 17% branching ratio each. However, the τ ± lepton is very heavy with mass at 1.78 GeV. The electron and muon from tau decay is then much more energetic than the IBD signal. A simple cut on the visible energy can effectively veto the atmospheric tau neutrino backgrounds.
Neutral Current: In addition to the CC events, the neutral-current (NC) scattering is also a potential background. Since the neutral current process is flavor blind, all flavors can contribute. First, the NC-QES scattering can not contribute as background since there is no positron in the final state. The NC-RES and NC-DIS processes both allow a single π ± and π 0 production in the final state. The interaction can also produce neutron and proton. Since π 0 decays into a pair of photons each with energy larger than m π /2 ≈ 70 MeV, the deposit energy is already beyond the IBD window of µDARν e . Among the remaining π ± , π − is absorbed by nuclei and only π + below Chereknov threshold can decay through π + → µ + → e + to fake the IBD positron. For NC-RES (NC-DIS) events, around 43% (45%) have a single neutron while 16% (22%) have a single charged π + in the final state. The fraction reduces to 11% (18%) if requiring both neutron and π + . It has a probability for E vis being within the energy window of interest as shown in Fig. 5. Requiring E vis ∈ [30,55] MeV, the integrated event number gives around 0.03 (0.06) events per year for NC-RES (NC-DIS) at THEIA-25 which is also a negligible amount. Since there is also proton, the background can be further reduced by at least one order if the proton kinetic energy is also taken into account. Fig. 6 shows the three survival background spectra as discussed above, including the beamν e with 1.83 × 10 25 protons on target (POT) as well as the atmosphericν e and ν µ /ν µ . While the atmospheric invisible muon background dominates at a Cherenkov detector such as Super-K and Hyper-K [87], it is a small minor contribution at the WbLS detector THEIA even with a conservative selection procedure. The intrinsic µDAR beam background is even smaller. THEIA is an ideal detector for the CP measurement with µDAR source.

Event Selection
As elaborated above, the IBD signal is quite distinctive at the THEIA detector. All backgrounds can be suppressed to negligibly small amount. Here we summarize the selection criteria for the IBD events, (1) Only one e-like Cherenkov ring; (2) For the total visible energy E vis , the number of scintillation photons N sc is within the range of [3900, 7150] and Cherenkov photons N ch ∈ [2400, 4400] that correspond to a (30 ∼ 55) MeV positron; (3) The ratio of Cherenkov and scintillation photons N ch /N sc is larger than 0.4; (4) Existence of delayed γs from neutron capture; For the above requirements, the energy window used in the criterion (2) can remove reactor, solar, geo-, supernova and DUNE beam neutrino backgrounds. The cri- teria (2) and (4) as double coincidence is efficient in removing the µDAR flux background except the intrinsic ν e . The combination of all criteria forms the triple coincidence which is even more powerful to reduce the invisible muon background. Finally, the criteria (1) and (2) can remove the atmospheric CC-RES, CC-DIS, and NC backgrounds.
Note that these criteria are already quite conservative. For the WbLS technique to be used by THEIA, more information can help to distinguish signal from background. Especially, the time information and pulse shape can be used to distinguish the Cherenkov and scintillation lights [94]. This could be extremely useful to further suppress the atmospheric invisible muon background with triple coincidence fully implemented.
In Fig. 7, we show the signal and background event rates for the low-energy mode at µTHEIA. The WbLS can highly suppress the background especially for the "invisible muon" that is reduced to a negligible amount. With CP dependence dominating the event rate, we can expect improvement of the CP sensitivity, which is elaborated later in Sec. IV.

B. The High-Energy Mode
The high-energy LBNF neutrinos from Fermilab to SURF is a wide beam spanning from 0.5 GeV to 5 GeV with peak at 2.5 GeV [44] for both neutrino and antineutrino modes. Each mode contains four different flavor components: ν e ,ν e , ν µ andν µ . Since the THEIA detector is at the same SURF experimental site as the DUNE far detectors, it can also probe the LBNF neutrinos. While the event reconstruction of high energy neu-trinos at the DUNE liquid Argon detectors has already been studied carefully [122], we focus on the detection at THEIA.
The key element for the neutrino CP measurement is neutrino flavor reconstruction for ν µ → ν e (ν µ →ν e ). It is interesting to see that the low-energy µDAR neutrinos outside the energy window [0.5 GeV, 5 GeV] cannot contribute as background. Although atmospheric neutrinos can overlap in energy, the pulse shape of the LBNF beam provides an efficient way to suppress the atmospheric backgrounds. Both signal and the remaining background actually come from the same LBNF beam.
With broad energy range, several types of CC scatterings with a target nuclei N can happen. In addition to a charged lepton α , the final state is either a single nuclei N for the quasi-elastic (ν α + N → α + N ), nuclei plus mesons for the resonant (ν α + N → α + N + meson), or nuclei plus hadrons for the deep-inelastic (ν α +N → α + N + hadrons) CC scatterings [117]. Typically CC-QES dominates below 1 GeV and CC-RES between 1.2 GeV and 5 ∼ 7 GeV, while CC-DIS takes over above 7 GeV. With the LBNF neutrino beam being below 5 GeV, most of the interactions are CC-QES and CC-RES. In order to make better neutrino reconstruction, it is desirable to distinguish these different CC scattering events.
The following discussions focus on the signal of appearance channels ν µ → ν e andν µ →ν e . Similar procedures shall also apply for the disappearance channels ν µ → ν µ andν µ →ν µ . Although the disappearance channels do not contribute significantly to the leptonic CP measurement, they can serve as supplementary probe of the other oscillation parameters and the neutrino flux. Both appearance and disappearance channels are taken into account in our GLoBES simulation in Sec. IV.

The CC-QES Category
Signals -The CC-QES process has a two-body final state with a primary lepton and a nuclei. A combination of Cherenkov and scintillation lights can achieve outstanding lepton identification as discussed in Sec. III A. Similar to the IBD case (5), the neutrino energy can be reconstructed from the charged lepton energy E (or momentum |p |) and scattering angle θ [104], where m f is the final-state nucleon mass. On the other hand, the initial nucleon mass m i always appears together with the binding energy E b for a nucleon inside 16 O nuclei as m i ≡ m i − E b . The binding energy E b = 42 MeV (19 MeV) for the neutrino (anti-neutrino) mode is adopted to make the reconstructed energy peak at the true value for E ν = 2.5 GeV as shown in Fig. 8. This energy reconstruction formula is similar to the IBD one (5) with the only difference that the initial free proton (Hydrogen) mass m p is replaced by m i to account for the binding energy. Experimentally, there is no efficient way to distinguish whether it is the Oxygen or Hydrogen nuclei that is scattered. Consequently, the neutrino energy reconstruction (7) for Oxygen target is used universally since the mass fraction of Oxygen is eight times larger than Hydrogen in the water target. Similar to the low energy mode, the same Gaussian smearing 7%/ E vis /MeV is used for simulating the detector resolution of the deposit visible energy. In addition, we take the angular resolutions 1.48 • for e ± and 1.00 • for µ ± from SK-IV [123]. Fig. 8 shows the reconstructed spectrum for a 2.5 GeV ν e (ν e ) at the THEIA detector. For the CC-QES events, the reconstructed neutrino energy spectrum is the narrowest whose width at half height is roughly 80 MeV. The high-energy CC-RES or CC-DIS event always has at least one π ± /π 0 in the final state which can be used to distinguish from the CC-QES signal. That means the CC-QES events can be separately from the CC-RES and CC-DIS counterparts. For illustration, Fig. 8 also shows the individual features of CC-RES and CC-DIS events separately. The fact that CC-RES and CC-DIS partially overlap with each other will be elaborated in later discussions.
Backgrounds -The LBNF beam can contribute intrinsic ν e (ν e ) background ν e → ν e (ν e →ν e ) which is irreducible. For the neutrino (anti-neutrino) mode, the ν e (ν e ) flux is two orders smaller than the dominant ν µ (ν µ ) [122]. Considering the ν µ → ν e (ν µ →ν e ) oscillation probability that is typically 5%, the intrinsic background can be comparable or roughly one order smaller than the signal. In other words, the intrinsic beam background is sizable and hence cannot be neglected. Moreover, the beam also has intrinsicν e (ν e ) fluxes which are smaller by another order. Although the positron annihilation at the THEIA detector can in principle be used to distinguish the electric charge of the final-state leptons, we assume theν e (ν e ) also contributes as background to be conservative. The intrinsic beam backgrounds are shown as the pink filled region in Fig. 9.
For the disappearance channel of ν µ (ν µ ), the misidentification of muon as electron is also a potential background. Nevertheless, the misidentification rate is 0.05% for a muon misidentified as an electron and 0.02% for the converse one [123]. This contribution is much smaller than the intrinsic ν e /ν e background and hence can be neglected for simplicity.
As for the neutral-current events, the NC-QES does not have charged lepton in the final state which can serve as an effective veto. However, the NC-RES and NC-DIS scatterings have both single and multi-pion final states. The multi-pion background can be removed by simply imposing a single-ring cut. Single pion production has π ± or π 0 in the final state. The π ± above the Cherenkov threshold can produce a µ-like ring, since both are heavy particles with similar masses. For π + below the Cherenkov threshold, it can experience a chain decay π + → µ + → e + with both π + and µ + decaying at rest. Although the final e + can produce an e−like ring, its maximal energy is only m µ /2 = 53 MeV and hence the reconstructed energy is far below the signal energy window [0.5, 5] GeV. Moreover, the invisible π + and µ + can actually been seen by the THEIA detector with scintillation lights before the e + Cherenkov ring, which provides an extra way to veto the invisible π + background. Furthermore, the invisible π − is absorbed in the detector and hence cannot produce an electron to fake signal. Neither π + nor π − can become background. However, the photon pair from π 0 decay can fake an electron. This happens if (1) the two photons have a small opening angle θ γγ ≤ 17 • and overlap with each other or (2) one photon is soft enough to be invisible [124]. Roughly 3.3% (16.2%) of NC-RES (NC-DIS) single π 0 events at the peak energy E ν = 2.5 GeV have overlapping photons. For the soft photon in case (2), we take a conservative E γ < 30 MeV Cherenkov detection threshold [124]. The soft photon events with θ γγ > 17 • contribute 7.6% (6.8%) of the NC-RES (NC-DIS) scatterings. In total, 10.9% (23.0%) of the NC-RES (NC-DIS) events have a π 0 to fake the electron. The π 0 direction and energy are then used as the electron information to reconstruct the neutrino energy via (7) for small opening angle case. Since the soft photon direction is difficult to reconstruct, the direction of the hard photon is used instead.
Event Selection -As elaborated above, the CC-QES category has only one primary lepton in the final state. Being different from the IBD signal of the low-energy mode, there is no neutral for double coincidence to significantly suppress the background. The selection criteria for the CC-QES signal events are summarized below, The event spectra of the CC-QES category after selection and the corresponding backgrounds are shown as functions of the reconstructed neutrino energy E rec ν in Fig. 9. For both neutrino and anti-neutrino modes, the signal lines clearly show the CP dependence. Taking δ D = 90 • for illustration, the corresponding blue line is at the bottom for the neutrino mode while it is at the top for the anti-neutrino one. The opposite feature happens for the other maximal CP phase δ D = 270 • . This is a reflection of the fact that the CP-violating term sin δ D in (3) differs by a sign. The NC background is contributed by all flavors and the CC one by the intrinsic ν e /ν e beam background in addition to the misidentified ν µ /ν µ comes from the disappearance channels. While the signal and CC backgrounds can extend to 5 GeV, the NC background mainly contributes below 3.5 GeV since typically more particles are produced to split the energy. The two NC background peaks here are produced due to the soft photon contribution (lower energy peak) and the small opening angle contribution (higher energy peak), respectively. For CC-QES, the CC background spectrum is quite flat and the NC backgrounds dominate below 2.5 GeV. Around the peak energy, E rec ν = 2.5 GeV, the signals dominate to provide a good CP sensitivity. Between neutrino and anti-neutrino modes, the typical event rate differs by a factor around 2.5.

The CC-RES Category
Signals -The CC-RES process for the appearance signal ν µ → ν e (ν µ →ν e ) has 1 electron (positron) plus 1 pion (π ± or π 0 ). Therefore, the CC-RES event can be distinguished from the CC-QES one by ring counting if the final-state pion is a charged one π ± . For π 0 , the major contribution to its decay final state is two resolved photons. Nevertheless, both CC-RES and CC-DIS can produce pion in the final state to have some overlap with each other. But it is still possible to partially distinguish CC-RES CC-DIS events by the number of pions. Our simulation shows that for the neutrino (anti-neutrino) mode at the peak energy of 2.5 GeV, only 15% (11%) of the CC-RES events have multiple pions, while CC-DIS reaches 64% (48%). Hence the CC single-pion event can categorized as CC-RES while the multi-pion one as CC-DIS.
The energy reconstruction formula (7) for CC-QES can no longer apply due to different particles in the final state. Instead of having nucleons in the initial and final states, CC-RES has heavy ∆ baryons as intermediate resonance. Most of the pions are produced from ∆ decays ∆ ++ → p + π + , ∆ + → p + π 0 , and ∆ − → n + π − [120]. Usually the proton is not energetic enough to produce a Cherenkov ring and hence cannot be uniquely identified but only leaves some scintillation lights. Observationally, the ν e /ν e CC-RES has a primary lepton and a pion (π 0 or π ± ). A reasonable neutrino energy reconstruction for CC-RES is [104], Comparing with (7), the final-state nucleon mass is replaced by the ∆ baryon mass, m ∆ = 1.232 GeV. As demonstrated with blue lines in Fig. 8, the CC-RES energy reconstruction formula (8) gives a correct peak position. However, for the high-energy scattering process, not only ∆ is produced but also other resonant particles like N (1440) which can also produce pion particles. So (8) cannot describe all the resonant processes exactly and gives a wider distribution than CC-QES.
Backgrounds -The beam background for the CC-RES category has three major contributions. First, the beam electron-flavor neutrinos (anti-neutrinos) of the disappearance channel ν e → ν e (ν e →ν e ) via either CC-RES or CC-DIS scatterings contributes as irreducible background with the same 1e1π final state. Being depicted as the pink regions in Fig. 10, this type of background for neutrino (left) and anti-neutrino (right) modes contribute roughly 13.9% and 14.9% of the total detected events for δ D = −90 • . The second beam background comes from the muon neutrinos ν µ interacting via CC-RES or CC-DIS to produce µ − and π 0 in the final state. The π 0 can fake electron if one of the decay photons is soft (E γ < 30 MeV) or the two photons are almost collinear (θ γγ < 17 • ). For the peak energy E ν = 2.5 GeV, only 7.6% (13%) of such π 0 can fake an electron. In addition, the muon needs to be misidentified as a charged pion. With both misidentification, the beam ν µ can fake the CC-RES signal. Although the new reconstruction algorithm of fitQun [123] shows some capability of separating pion from muon at T2K using the hadronic kinks in the pion propagation [125], we assume µ − and π + can not be separated to be conservative. Not to say, µ − has an electron but π + has a positron in their decay products, respectively, and the e + e − annihilation photons from π + can also provide a distinguishable feature. Note that there is no beamν µ CC background for the anti-neutrino mode. This is because the π − in theν e CC-RES signal is absorbed in the detector and does not produce a delayed Michel electron which is different from the beamν µ that can decay to µ + and finally a Michel e + . For both neutrino and anti-neutrino modes, muon misidentified as electron/positron might happen together with the same π ± to fake the CC-RES category. But this background can be neglected since the misidentification rate is negligibly small (< 0.05%) as mentioned in Sec. III B 1. The beam ν µ CC background is shown as green region in the left panel of Fig. 10 only for the neutrino mode. The third part of the background comes from NC-RES and NC-DIS with two pions (2π 0 or π 0 + π ± ) in the final state. In such case, one π 0 decays into unresolved photons to fake the primary electron signal while the other pion coincides with the one in the signal. The estimated NC background is shown as the orange region of Fig. 10. We can see that the reconstructed NC background mostly occur at E rec ν < 3 GeV.
Event Selection -In summary, the selection criteria for CC-RES signal is, (1) Only one primary e-like Cherenkov ring with |p e | > 100 MeV.
(3) Single pion particle in the final state. To be conservative, the events with a charged π ± or a neutral π 0 are not mixed into a single CC-RES category. More detailed study is necessary for maximizing the CP sensitivity. Fig. 10 depicts the ν e (left) andν e (right) CC-RES signal and the corresponding backgrounds. Comparing the two panels, we can see that the neutrino event rates are typically 4 ∼ 5 times larger than its anti-neutrino counterparts. This is a direct consequence of the relatively larger cross section for neutrino than anti-neutrino. For backgrounds, the neutrino mode has four components while the anti-neutrino mode has only three. This is because the disappearance channel ν µ →ν µ has µ + and Michel e + to be distinguished from π − of theν e CC-RES signal. For neutrino mode, all backgrounds contribute roughly the same size while for anti-neutrino the NC and intrinsic ν e background dominate. Around the major peak, E rec ν ≈ 2 GeV, the signal is slightly larger than the background for neutrino and anti-neutrino modes.

The CC-DIS Category
Signals -In order to enhance the statistics, we also include CC-DIS events as a separate category. As mentioned in the previous CC-RES section, the events identified as CC-DIS are those with an e-like ring accompanied with more than one pion in the final state. The majority of events contains two pions. The 2-pion events include any of the following combination (π ± , π ± ), (π ± , π 0 ) and (π 0 , π 0 ). Similar to the CC-RES category, only a pair of resolved photons are identified as π 0 . For events with one or more π 0 and consequently 3 or more e-like rings in the final state, we take the most energetic e-like ring as the primary electron/positron. Then the neutrino energy can be reconstructed using the same (8). Nevertheless, the scattering process of CC-DIS differs a lot from CC-RES and the reconstructed neutrino energy is almost flat as shown with green curves in Fig. 8.
Backgrounds -The background of the CC-DIS process also has three major contributions. The first is the irreducible one from the electron-flavor neutrinos via the CC-RES or CC-DIS with the 1eN π final state with N > 1. Being depicted as the pink regions in Fig. 11, this background is the largest contribution for E rec ν > 1.5 GeV and can be as large as 45% of the total events at 2.5 GeV for δ D = −90 • .  The second component comes from muon neutrinos ν µ interacting via CC-RES or CC-DIS to produce a µ − , a π 0 and one or two other pions of any type. While µ − is misidentified as a pion, the π 0 can be misidentified as an electron if it has a soft photon E γ < 30 MeV or the two photons are almost collinear, θ γγ < 17 • . The multi-pion production of two or three pions occurs for 4% of the total ν µ CC events. In addition, the misidentification of the π 0 as an electron occurs in 17% (9%) of the events for 2-pion (3-pion) final state. This ν µ background is shown as green curve in the left panel of Fig. 11 only for the neutrino mode. Again, there is no ν µ CC background due to the delayed Michel electron veto.
The third component is the NC multi-pion background of 1π 0 1π x 1π y with x, y = 0, ±. The triple-pion final state is reconstructed as background if one π 0 is misidentified as electron. This occurs for only 0.1% of the total NC events at the 2.5 GeV peak energy, which is represented as the orange region in Fig. 11. It is most relevant at low energies, E rec ν < 1.5 GeV where the CP value has small impact.
Event Selection -We can see that the CC-DIS category has much similarity as the CC-RES one with singlepion final state replaced by the multiple one. Although the CC-DIS signal would not contribute much as shown in Fig. 8, this separation reveals the feature of various channels. Below we summarize the selection criteria for the CC-DIS signal, (1) Only one primary e−like Cherenkov ring with |p e | > 100 MeV. In the presence of multiple e-like rings, the most energetic one is identified as the primary lepton.
(3) At least two pions in the final state with energy π(|p π |) > 200 MeV. Similar to the CC-RES category, all CC-DIS events are put into a single category without division according to the final-state pions. This conservative treatment can be further improved with more careful studies.
The CC-DIS signals and their backgrounds are shown in Fig. 11. Different from the CC-QES in Fig. 9 and CC-RES in Fig. 10, the CC-DIS signal can be smaller than the background, especially for the neutrino mode. Around the peak, E res ν ≈ 2.5 GeV, the dominant background comes from the CC-RES multi pion background which can reach more than 50% of the total events at E rec ν ≈ 2 GeV. The anti-neutrino case is slightly better with signal still dominating around the peak energy. In addition, the three background components for theν e mode all have sizable contributions. The NC background background dominates for E rec ν < 2 GeV, while the CC-RES dominates at intermediate energies 2 GeV < E rec ν < 3.5 GeV and the intrinsic ν e dominates for E rec ν > 3.5 GeV. Nevertheless, ν e events are 7 times larger thanν e .

IV. CP SENSITIVITY WITH µTHEIA AND DUNE
As discussed in Sec. II, the combination of µTHEIA and DUNE is expected to improve the CP measurement. Below we give a quantitative estimation. We first summarize the experimental setups in Sec. IV A for complete-ness and then establish the χ 2 formalism in Sec. IV B. In order to maximize the CP sensitivity, Sec. IV C explores the optimal baseline between the µDAR source and the THEIA detector for µTHEIA. Based on these, the CP sensitivity and the influence of matter effect are studied in Sec. IV D.

A. Experimental Setups
The µDAR neutrinos are typically produced by cyclotrons [126]. We adopt the configuration that the µDAR flux is generated by a 9 mA proton beam with 800 MeV protons hitting the high-Z target to deliver 1.83 × 10 25 POT in 10 years [87]. Since the µDAR cyclotron is not a pulsed beam, it is possible to achieve almost full duty factor. Of the charged pions produced by proton hitting the target, π − is mainly absorbed by the positively charged nuclei while π + first loses energy in the thick material and then decays at rest to produce µ + . Similar process of decay at rest also happens for µ + . Consequently, the µDAR neutrino spectra are well defined and predicted by the SM interactions.
For DUNE, the LBNF neutrinos are produced by 120 GeV protons with 1.2 MW beam power and the flux can reach 1.1 × 10 21 POT/year. We take 6.5 years running for each neutrino and anti-neutrino mode [44]. With completely different energy windows [30,55] MeV for µDAR and [0.5, 5] GeV for LBNF, there is no energy overlap between the low-energy µDAR flux and the highenergy LBNF one. So the two fluxes can run simultaneously.
For the THEIA detector we consider two possibilities, THEIA-25 with a fiducial mass of 17 kt and THEIA-100 with a fiducial mass of 70 kt [51]. The DUNE detector has four modules, each with 10 kt fiducial mass [122]. Both THEIA and DUNE detectors are at the same SURF site and 1289 km away from the LBNF source. To install the THEIA detector, one of the DUNE modules needs to be replaced. As mentioned in Sec. II, the average matter density 2.85 g/cm 3 is used for both DUNE and µTHEIA.
In the following, we consider three combinations: (1) the full DUNE configuration with 4 far detector modules, (2) the reduced DUNE with 3 far detector modules and µTHEIA-25, and (3) the reduced DUNE and µTHEIA-100. The high-energy mode (HEM) can be detected by both the DUNE and THEIA detectors while the lowenergy mode (LEM) only applies for the THEIA detector as discussed in Sec. III. The energy window, fiducial mass, baseline, and running time of each experimental configurations are summarized in Tab. I.

B. Simulation and χ 2 Analysis
We then use GLoBES [127,128] to simulate the event rates and evaluate the CP sensitivities. A quantitative evaluation is realized by minimizing the χ 2 function that contains three contributions, for statistical (χ 2 stat ) and systematical (χ 2 sys ) uncertainties in addition to the prior constraint on the oscillation parameters (χ 2 para ). For event rate > 10 in a single bin, a Gaussian χ 2 is much more convenient with analytical fit [100,129,130]. However, the event rates considered in this paper are not large enough. To make the sensitivity evaluation exact, the statistical part in the first term takes the Poisson form, where N data i , N sig i and N bkg i are the pseudo data, signal and background event numbers in the i-th bin, respectively. The coefficients a sig and a bkg are nuisance parameters for the signal and background normalizations, respectively.
The χ 2 sys term contains the uncorrelated Gaussian priors of signal and background normalizations. We take σ sig = 5% and σ bkg = 10% for each channel of both the low-and high-energy modes that are observed at the THEIA detector. For the low-energy mode, the normalization uncertainties are the same as the configuration given in [87] while the uncertainties are more conservative than the values (2% for ν e and 5% forν e ) used in [51]. The systematics of the LBNF beam detection by the DUNE detector are described by the official configuration files [122].
Among these oscillation parameters, the solar mass squared difference ∆m 2 s and the solar mixing angle θ s are kept fixed throughout our analysis. On one hand, the contribution of the solar mass squared difference ∆m 2 s enters via a coefficient parameter α ≡ ∆m 2 s /∆m 2 a as shown in (3). The parameter α has a small value (≈ 0.03) and the error in ∆m 2 s has an even smaller contribution. Hence it can be neglected comparing with the O(1) uncertainty on the CP phase. On the other hand, the solar mixing angle appears in the first and the third term on the righthand side of (3). Since the current prior on the solar mixing angle is roughly 3%, it can also be neglected for the study of the large CP phase uncertainty. Moreover, the next generation reactor neutrino experiment like JUNO [63] will provide a sub-percent uncertainty on θ s that is negligibly small.
For the other mixing parameters, the reactor mixing angle θ r , the atmospheric mixing angle θ a , and the atmospheric mass-square difference ∆m 2 a are treated as free parameters. We use the marginalized one-dimensional χ 2 curves [131] as our priors, χ 2 para ≡ χ 2 θr +χ 2 θa +χ 2 ∆m 2 a . Due to the existing tension between the T2K and NOνA results, which is discussed in Sec. I, we do not include any prior on δ D .
Since the matter effect is a natural source of fake CP, we take its uncertainty by including a parameter q to scale its average value defined in Sec. II, A → qA [132].
The average matter density corresponds to q = 1 while the vacuum case takes q = 0. Note that q is undetermined and hence treated as the fifth free parameter. In the following discussions, we consider two different scenarios for its uncertainty: fixed q = 1 (no uncertainty) and a conservative 10% Gaussian uncertainty. For comparison, previous studies have used 1% ∼ 2% uncertainty [45,[133][134][135].
The CP uncertainty ∆δ D is defined as the half-width of the ∆χ 2 = χ 2 (δ D ) − χ 2 min = 1 band where χ 2 min corresponds to the best-fit value δ BF D of the CP phase. Since our simulation uses pseudo-data, the best-fit value is the same as the true value, δ true D . Note that it is not necessary for the χ 2 (δ D ) function to be symmetric around the minimum. For this case, the previous definition, ∆χ 2 = 1, gives two boundaries below (δ − D ) and above (δ − D ) the bestfit value δ BF D . Then we take the average deviation as the CP uncertainty, ∆δ D ≡ (δ + D − δ − D )/2.

C. Baseline Options of µTHEIA
The baseline between the µDAR source and the THEIA detector can significantly affect the CP uncertainty ∆δ D . Fig. 12 shows ∆δ D as a function of the µTHEIA baseline L in the range from 10 km to 80 km. Since the true value of the CP phase δ true To see the impact of matter effect on the optimal baseline, the left panel of Fig. 12 is obtained by fixing q = 1 while the right one takes q = 1 ± 0.1. For both cases, the result shows two local minima in the CP uncertainty for maximal CP violation. One is around L = 38 km and the other around L = 55 km. The longer one, L = 55 km, is a local optimal option for δ D = 270 • that is preferred by the T2K measurement. Although the true local minimum for vanishing CP violation cases δ D = 0 • and 180 • actually happens with L 65 km, the difference in χ 2 is not significant while the maximal CP violation cases δ D = 90 • and 270 • (or equivalently −90 • ) become much worse. Since the data-driven δ D = −90 • is of larger interest, L = 55 km is preferred than the longer 65 km. For the shorter one, the choice is more difficult. The global minimum around L ≈ 38 km for δ D = ±90 • is very close to the global maximum for the vanishing CP violation cases. So choosing L = 38 km needs to pay too much price and we take L = 30 km to balance among various CP values. Our simulations takes these three baselines L = 30 km, 38 km, and 55 km as possible options. The final choice is up to the on-going T2K and NOνA experiments. The comparison between the left (fixed q = 1) and right (10% uncertainty around q = 1) panels of Fig. 12 shows that although the uncertain matter effect contaminates the CP sensitivity, increasing ∆δ D by 3 or 4 degrees to be exact, it does not affect the optimal baselines. The optimal baseline options L = (30,38,55) km for DUNE + µTHEIA are all different from the TNT2HK one L = 23 km [87]. Not just the baseline length is different, but also TNT2K/TNT2HK obtains only a single local minimum. The key difference is the atmospheric invisible muon background. Since both SK and HK are water Cherenkov detectors, the atmospheric invisible muon dominates the background. As L increases, the beam flux decreases with 1/L 2 , but the atmospheric background remains the same and eventually dominates the statistics. So the relatively large amount of the atmospheric background limits the optimal baseline length. For comparison, THEIA with WbLS reduces the invisible muons to negligible amount as shown in Fig. 7. Therefore, µTHEIA can have longer baseline than TNT2K/TNT2HK to optimize the CP sensitivity.
While the options L = 30 km and 38 km are not so far from the 23 km of TNT2K/TNT2HK and hence easier to understand, the longer baseline L = 55 km also achieving comparable CP sensitivities seems counter-intuitive. From 30 km to 55 km, the flux decreases quadratically with distance and is suppressed by a factor (30/55) 2 ≈ 0.3. However, this flux reduction is compensated by the increasing oscillation amplitude. To make this feature explicit, we show in the left panel of Fig. 13 the oscillation probability difference, ∆Pν µ→νe (≡ P fit νµ→νe −P truē νµ→νe ), between the true (δ true D ) and fit (δ fit D ) CP values, ∆Pν µ→νe = 16∆ s s r c s s s c a s a sin ∆ a Note that this formula is obtained from (3) in the vac-uum limit (A → 0). The oscillation and CP phases are defined as ∆ a,s ≡ ∆m 2 a,s L/4E and ∆ δ D ≡ δ fit D − δ true D . To further illustrate the CP sensitivity, the parameter (∆Pν µ →νe ) 2 /P truē νµ→νe that has a similar form as the χ 2 calculation is also calculated in the right panel of Fig. 13. For convenience, we name (∆P ) 2 /P as pseudoχ 2 at the oscillation probability level. Both variables are shown as a function of the atmospheric oscillation phase ∆ a (≡ ∆m 2 a L/4E ν ) for five different true CP values, δ true D = 0 • , 90 • , 180 • , 270 • , and 150 • . The difference between the fitting and true CP values are assigned to have ∆ δ D = 10 • for illustration.
The left panel shows that the second peak of the oscillation probability difference ∆P is much larger than the first one and the relative size between the pseudo-χ 2 peaks in the right panel is also significantly enhanced. Take the δ true D = 0 • curve in the right panel as an example, the second peak is 3 times of the first one which can roughly compensate the flux suppression (∼ 0.3). Similar feature applies also for the other true CP values. This explains why the local minimum around L = 55 km has roughly the same value at L = 30 km.
The three filled regions correspond to different baselines L = 30 km, 38 km, and 55 km, respectively, while the µDAR neutrino energy E ν spans a wide range of [30,55] MeV. With ∆ a inversely proportional to E ν , the left boundary of each region corresponds to the upper energy limit 55 MeV and the right one to the lower limit 30 MeV. Since theν µ spectrum from µDAR source peaks at the upper limit [87], the left sides of the filled regions give the largest contribution. This important feature can explain the location of those local minimums in Fig. 12. For example, the left side of the pink region for L = 30 km covering the first oscillation peak/valley of δ true Fig. 13 explains why this baseline corresponds to the best sensitivity of this true CP value in Fig. 12. The same thing happens for L = 38 km with the left side of the purple region covering the δ true D = ±90 • peak/valley in Fig. 13 to justify the local minimum in Fig. 12. Not to say the left side of the light green region of L = 55 km covers the second peak/valley of δ true D = ±90 • in Fig. 13.

D. CP Sensitivity and Matter Effect
As emphasized in earlier discussions, the CP sensitivity suffers from matter effect contamination. The DUNE + µTHEIA configuration we propose in this paper can overcome this issue to provide a clean measurement of the Dirac CP phase δ D . Fig. 14 shows the CP uncertainty ∆δ D as a function of the true value δ true   FIG. 14. The CP phase uncertainty ∆δD as a function of the true CP value δ true D . For illustration, the three µTHEIA baseline options L = 30 km (top), 38 km (middle), and 55 km (bottom) are shown separately. Each panel takes three different experimental setups, DUNE with 4 modules (blue), DUNE with 3 modules and µTHEIA-25 (red) or µTHEIA-100 (green). As for the matter effect, both the fixed case (q = 1, solid) and 10% uncertainty (q = 1 ± 0.1, dashed) are implemented to show how it affects the CP sensitivity. In all panels, the DUNE configuration contains the full 40 kt detector and the blue curves do not change with varying L since it does not depend on µTHEIA baseline. For fixed matter effect (q = 1), the CP uncertainty peaks around δ true D ≈ 70 • and 250 • which are consistent with [65,133,135,136] where either fixed matter effect or just 2% uncertainty is adopted. Slightly shifted peaks at δ true D ≈ 90 • and 270 • are obtained and attributed to different treatment of systematics in [137].
From DUNE alone to DUNE + µTHEIA-25, the CP uncertainty ∆δ D significantly reduces by roughly 1/4. This is especially true around the maximal CP values, δ true D = ±90 • since the µDAR spectrum is especially wide to provide both cos δ D and sin δ D terms. With only sin δ D term, the CP uncertainty around the maximal CP phase is intrinsically large, ∆δ D ∝ 1/ cos δ D . But a wide spectrum can also introduce a large enough cos δ D term to make the CP uncertainty decrease. Previous study shows that this feature is expected to appear when µDAR flux is added to supplement the narrow beam accelerator experiments, such as TNT2K/TNT2HK [40][41][42]87]. Nevertheless, the result turns out that this is also true for the addition of µTHEIA to DUNE, although the LBNF flux spectrum is already quite wide. Adding a larger µTHEIA-100 can even further reduce the CP uncertainty to almost only 5 • for the best case.
As expected, the matter effect can fake the CP violation and hence its uncertainty can significantly modify the CP uncertainty. In addition to the fixed q = 1 scheme (solid lines), Fig. 14 also shows the results obtained with 10% uncertainty in the matter effect or equivalently q = 1 ± 0.1 (dashed lines). With q relaxed, the CP uncertainty becomes much worse especially around the vanishing CP violation, δ true D ≈ 0 • or 180 • . The most significantly affected point is around δ true D ≈ 150 • or 330 • . With the uncertainty of matter effect taken into account, the improvement brought by µTHEIA is even more significant. This is exactly because of the fact that matter effect plays more important role at DUNE with much higher energy than the low-energy µTHEIA. With µTHEIA added, even switching on matter effect uncertainty would not make the situation much worse. It also is interesting to see that the CP uncertainty around the maximal CP violation is almost not affected by switching on/off the matter effect uncertainty. The µTHEIA improvement is quite stable against matter effect.
To understand the interplay between the Dirac CP phase and the matter effect qualitatively, Fig. 15 shows the partial derivative ∂P/∂δ D of oscillation probability (3) with respect to the Dirac CP phase δ D and ∂P/∂q to the matter effect parameter q at q = 1 as a function of neutrino energy E ν . The matter effect mimics the CP effect quite well at δ true D = 150 • with the two derivative curves having similar shapes and peak positions. For comparison, the two derivatives have very different features at 70 • .
Note that the green curves with DUNE and µTHEIA-100 are much more flat than the original DUNE alone. For most of the parameter space, the CP uncertainty is better than 8 • . Among the three panels of Fig. 14, the L = 38 km one has the most flat CP uncertainty curves for the DUNE + µTHEIA-100 configuration. Especially, the CP certainty is always better than 8 • no matter what is the value of δ true D . In this sense, L = 38 km is probably the optimal baseline for µTHEIA.
To further illustrate the advantages of the DUNE + µTHEIA combination, we compare with other existing experiments or designs in Fig. 16. The first two rows show the latest measurements from the NOνA (δ D = 148 •+49 • −157 • ) [25] and T2K (δ D = −108 •+40 • −33 • ) [24]. The T2K experiment has two major upgrades: T2HK [43] with a much larger Hyper-K detector and T2HKK [56] with another detector at the second oscillation peak. For both of them, a matter effect uncertainty of 6% is considered [56]. Since T2HK is in construction and T2HKK still being planned, there is no real data yet. We take two typical values δ D = −90 • and 150 • for illustration. At these two true values, the CP uncertainty at T2HK (T2HKK) can reach 22 • (13 • ) and 10 • (7 • ), respectively [43,56].
The next group is the accelerator-based DUNE [44] and MOMENT [59]. The CP uncertainty at DUNE is  derived from our own simulation with a conservative 10% uncertainty in the matter potential, as shown in Fig. 14. Different from the µDAR flux, the neutrino flux from muon decay in flight is adopted by the MOMENT experiment to 15 • [138].
Finally, TNT2HK and the two DUNE+µTHEIA-25/100 configurations use µDAR neutrinos to supplement the accelerator measurements. Since the Hyper-K detector is going to be built, TNT2HK will always dominate over TNT2K once a µDAR source is added around the Kamioka site. So we only show TNT2HK in Fig. 16. Among these three options, TNT2HK has the advantage of using the full Hyper-K detector including Super-K for the detection of both the accelerator and µDAR neutrinos. Nevertheless, the atmospheric invisible muon background limits its CP uncertainty to 12 • at regions around maximum CP phase. [87]. In addition, the J-PARC beam energy and flux [20] are also lower than the LBNF ones [44]. Although the µTHEIA-25 can only use a fiducial volume of 17 kt, the CP uncertainty 11 • (10 • ) at DUNE + µTHEIA-25 is already slightly better than TNT2HK. With 70 kt fiducial volume at µTHEIA-100, the CP uncertainty further reduces to only 7 • ∼ 8 • . This clearly shows the advantages of supplementing DUNE with µTHEIA-25 or µTHEIA-100.

V. CONCLUSION AND OUTLOOK
The leptonic CP phase measurement at acceleratorbased neutrino oscillation experiments suffers from the contamination of matter effect. The higher neutrino energy, the more severe contamination. In this paper, we put forward a possible combination of intrinsically low-energy µDAR neutrinos and the recently proposed THEIA detector to overcome this problem.
Our simulation shows that the THEIA detector using WbLS has very good capability of particle identification. This is especially useful for suppressing the atmospheric invisible muon background, which was the major background at the TNT2K/TNT2HK configuration, to negligible amount. Then the µDARν µ →ν e oscillation leaves a very clear IBD signal in the detector. In addition, the high-energy LBNF flux can also be measured at the THEIA detector in addition to DUNE.
With essentially a background free measurement, the enhancement on CP sensitivity from µTHEIA is significant. The CP uncertainty around the maximal CP violation δ D = ±90 • reduces up to 20% (40%) when compared to the standard DUNE configuration. Especially, the CP uncertainty is controlled to be below 8 • and the best case can be as good as 6 • for the baseline L = 38 km. In addition, the dependence of CP uncertainty on the true CP phase value is largely mitigated. If realized, either the DUNE + µTHEIA-25 or DUNE + µTHEIA-100 configuration can bring the CP measurement into a precision era. During the scattering process, the incoming neutrino with momentum pν (long red arrow) becomes an electron/positron with momentum pe (blue arrow). The neutrino momentum (−pν ) is parameterized by θz and φz in the lab frame (green) while the electron/positron momentum (pe) by θe and φe in the neutrino frame. Without knowing the incoming neutrino direction, the genuine scattering angle θe is wrongly reconstructed as θw. Fig. 17 shows the geometry of the atmospheric IBD background. With O(10) km baseline, the µDAR neutrinos essentially travels horizontally, providing a natural definition of the x-axis while the other horizontal direction is the y-axis of the lab frame (cyan). Then, the zenith angle θ z is measured from the vertical z-axis. Given θ z , the incoming atmospheric neutrino direction is parametrized by the azimuth angle φ z , measured from the x-axis. The final-state lepton typically has a nonzero scattering angle θ e from the neutrino direction. For convenience, a neutrino frame (blue) is established around the neutrino momentum p ν with x , y , and z -axes.
The wrong scattering angle is defined as the one between the direction of e − /e + (blue arrow) and the direction of the µDAR flux (cyan x-axis). For comparison, the true scattering angle θ e and the corresponding azimuth angle φ e are also plotted in the figure. To calculate the wrong angle as a function of θ e , θ z , φ e , and φ z , one need to first establish the connection between the horizontal and neutrino frames. The z direction can be easily read out from the figure, z = − sin θ z cos φ z x − sin θ z sin φ z y − cos θ z z. But the xand y -axes can be randomly set, as long as they satisfy n x · n z = n y · n z = n x · n y = 0, | n x | = | n y | = 1 and n x × n y = n z . Accordingly, the following transfor-mation matrix between two coordinate systems is chosen, The e − /e + direction in the neutrino frame is (sin θ e cos φ e , sin θ e sin φ e , cos θ e ). Using the frame transformation (A1), the e − /e + direction in the lab frame is