Flavor-dependent long-range neutrino interactions in DUNE & T2HK: alone they constrain, together they discover

: Discovering new neutrino interactions would represent evidence of physics beyond the Standard Model. We focus on new ﬂavor-dependent long-range neutrino interactions mediated by ultra-light mediators, with masses below 10 − 10 eV, introduced by new lepton-number gauge symmetries L e − L µ , L e − L τ , and L µ − L τ . Because the interaction range is ultra-long, nearby and distant matter — primarily electrons and neutrons — in the Earth, Moon, Sun, Milky Way, and the local Universe, may source a large matter potential that modiﬁes neutrino oscillation probabilities. The upcoming Deep Underground Neutrino Experiment (DUNE) and the Tokai-to-Hyper-Kamiokande (T2HK) long-baseline neutrino experiments will provide an opportunity to search for these interactions, thanks to their high event rates and well-characterized neutrino beams. We forecast their probing power. Our results reveal novel perspectives. Alone, DUNE and T2HK may strongly constrain long-range interactions, setting new limits on their coupling strength for mediators lighter than 10 − 18 eV. However, if the new interactions are subdominant, then both DUNE and T2HK, together, will be needed to discover them, since their combination lifts parameter degeneracies that weaken their individual sensitivity. DUNE and T2HK, especially when combined, provide a valuable opportunity to explore physics beyond the Standard Model.


Introduction
Discovering a new fundamental interaction would be striking evidence of physics beyond the Standard Model.Yet, because new interactions are likely feeble, they are difficult to detect.And because they may manifest in a variety of ways, they are difficult to search for comprehensively.So far, there is no evidence for them, despite a long history of searches, though there are stringent limits on their strength [11][12][13][14][15][16].
We consider flavor-dependent neutrino-matter interactions, originally introduced in Refs.[46][47][48][49][50], and explored and constrained in earlier literature, e.g., in Refs.[1-3, 7, 8, 38, 51-65].Two reasons motivate our choice.First, if these interactions are long-range, i.e., if they act across long distances, then large collections of nearby and distant matter may source a sizable matter potential that affects neutrino flavor oscillations appreciably.Thus, we concentrate on interactions mediated by new, ultra-light mediators, with masses below 10 −10 eV, that subtend ultra-long interaction ranges.Second, the flavor-dependent interactions we consider are born from gauging, anomaly-free, global symmetries of the Standard Model [66][67][68][69][70][71]: L e − L µ , L e − L τ , and L µ − L τ , where L e , L µ , and L τ are the electron, muon, and tau lepton numbers.This makes them arguably natural and economical extensions of the Standard Model.Gauging each one introduces a single new neutral vector boson that mediates new neutrino interactions with electrons or neutrons (interactions with other particles are suppressed, as we elaborate on later).
Previous works have explored the sensitivity of existing and future long-baseline neutrino experiments to flavor-dependent long-range interactions.However, they either fixed the interaction range, typically to be equal to the Sun-Earth distance (see, e.g., Ref. [56]), or considered mediator masses only as small as about 10 −18 eV (see, e.g., Ref. [1]).We abandon both limitations and explore mediator masses down to 10 −35 eV.Doing so opens up a largely unexplored regime of ultra-long-range interactions.As pointed out in Ref. [7], a mediator this light allows for electrons and neutrons in the Earth, Moon, Sun, Milky Way, and the cosmological distribution of matter to affect neutrino oscillations.To make our forecasts realistic, we base them on detailed simulations of DUNE and T2HK, including their different detection channels, efficiency, backgrounds, and run times.
Figure 1 conveys the novel perspectives revealed by our work.It shows the first half of our main results, concerning constraints: separately or, as in Fig. 1, together, DUNE and T2HK may place the strongest constraints on long-range interactions, especially for mediators lighter than 10 −18 eV .(Future sensitivity from flavor measurements of high-energy astrophysical neutrinos in the IceCube-Gen2 neutrino telescope might be comparable [7], but, for now, they are subject to large uncertainties in the neutrino flux, not pictured in Fig. 1, unlike the constraints from DUNE and T2HK.)The other half of our main results, not contained in Fig. 1, concerns discovery.We find that, separately, DUNE and T2HK will likely be unable to discover subdominant flavordependent long-range neutrino interactions, due to degeneracies between their effect on neutrino oscillations and that of the standard mixing parameters.Yet, together, their complementary capabilities may lift degeneracies and enable the discovery of the new interactions; see Fig. 7. Below, we elaborate on these perspectives.
This paper is organized as follows.Section 2 introduces lepton-number gauge symmetries, long-range interactions, and their effect on neutrino oscillations.Section 3 overviews DUNE and T2HK, and shows oscillation probabilities and event rates in them.Section 4 shows projected constraints and discovery prospects.Section 5 summarizes and concludes.2 Flavor-dependent long-range neutrino interactions

Gauged lepton-number symmetries
In the Standard Model (SM), the baryon number and the lepton numbers, L e , L µ , and L τ , are accidental U (1) global symmetries.Linear combinations of the lepton-number symmetries can be gauged anomaly-free, i.e., without introducing a new fermion or righthanded neutrino (although not simultaneously) [47][48][49][50].To showcase the capabilities of DUNE and T2HK, we explore three such new U (1) gauge symmetries, generated by L e −L µ , L e − L τ , and L µ − L τ , that introduce new flavor-dependent neutrino-matter interactions; later, we show how they affect neutrino oscillations.(Other combinations of baryon and lepton numbers can also be gauged anomaly-free; see Ref. [1].) Figure 2 shows the Feynman diagrams for neutrino-matter interaction that we consider.For a particular lepton-number symmetry, the corresponding effective Lagrangian is (2.1) The first term describes the SM contribution, mediated by the Z boson, i.e., where e/(sin θ W cos θ W ) = 0.723, e is the unit charge, θ W is the Weinberg angle, ν α and l α are a neutrino and charged lepton of flavor α = e, µ, τ , P L is the left-handed projection operator, and u and d are up and down quarks.Because the Z boson is heavy, the interaction that it mediates is short-range; in our work, it matters only inside the Earth.(Equation (2.2), and also Eq. (2.4) below, assumes that matter is electrically neutral, i.e., that it has equal abundance of electrons and protons [54], which is also what we assume later when computing the new matter potential; see Section 2.2.) The second term in Eq. (2.1) describes the interaction between ν α and l α mediated by the new Z ′ αβ boson [46,49,54], i.e., for the L α − L β symmetry, where g ′ αβ is a dimensionless coupling constant.Due to the dearth of naturally occurring muons and tauons with which neutrinos can interact, we neglect this contribution under L µ −L τ and consider it only under L e −L µ and L e −L τ , for which the interaction is sourced by comparatively abundant electrons.
The final term in Eq. (2.1) describes the mixing between Z and Z ′ αβ [54,60,72], which can arise directly or by radiative mixing [73,74].In the physical basis, this term is [72] , where χ is the kinetic mixing angle between the two bosons and ξ is the rotation angle between gauge eigenstates and physical states.This introduces a four-fermion interaction between neutrinos and charged leptons, protons, and neutrons via Z-Z ′ αβ mixing, i.e., where However, the contribution of electrons is nullified by that of protons, leaving only neutrons to source the new interaction via mixing.The term (ξ − sin θ W χ) effectively describes the strength of the Z-Z ′ αβ mixing.Its value is unknown, but there are upper limits on it [54,75,76].We do not consider its value independently, but together with g ′ αβ , as an effective coupling (more on this below).In order to showcase the effect of mixing, we include L mix only under L µ − L τ .
In summary, under L e −L µ and L e −L τ , the new interactions are described by L Z ′ , and are sourced by electrons only, whereas under L µ − L τ , the new interactions are described by L mix , and are sourced by neutrons only.In all cases, in addition, standard neutrinoelectron interactions, described by L SM , are active only inside the Earth.

Long-range matter potential
The above interactions induce flavor-dependent Yukawa potentials, sourced by electrons and neutrons, that affect the mixing of neutrinos [7,8,46,[48][49][50].Under L e −L β (β = µ, τ ), a neutrino located at a distance d from a collection of N e electrons experiences a potential where m ′ eβ is the mass of the mediating Z ′ eβ boson.Under L µ − L τ , a neutrino located at a distance d from a collection of N n neutrons experiences a potential where m ′ µτ is the mass of the mediating Z ′ µτ boson.In Eqs.(2.5) and (2.6), the effective coupling strength is At distances longer than the interaction range of 1/m ′ αβ , the potential is suppressed due to the mediator mass.Like Ref. [7], we explore light mediators with m ′ αβ = 10 −35 -10 −10 eV, corresponding to interaction ranges from 10 3 Gpc -larger than the observable Universe -to hundreds of meters; see Fig. 1.
We adopt the methods introduced in Ref. [7] to compute the total potential sourced by nearby and faraway electrons and neutrons in the Earth (⊕), Moon ( ), Sun (⊙), Milky Way (MW), and by the cosmological distribution of matter (cos) in the local Universe, i.e., The specific value of m ′ αβ determines the relative size of the contributions of the above sources to the total potential.We do not compute the changing potential along the underground trajectories of the neutrinos from source to detector inside the Earth; see Ref. [1] for such treatment.Instead, like Ref. [7], we compute the average potential experienced by the neutrinos at their point of detection.This approximation is especially valid for mediators lighter than about 10 −14 eV, for which the interaction range is longer than the radius of the Earth (see Fig. 1), and so all of the electrons and neutrons on Earth contribute to the potential experienced by a neutrino regardless of its position along its trajectory.Below 10 −14 eV is also where we place novel projected limits.
We assume that the matter that sources the potential is electrically neutral, so that the number of electrons and protons is the same, and isoscalar, so that the number of electrons and neutrons is the same, except for the Sun [54] and for the cosmological distribution of matter [77][78][79].We treat the Moon (N e, = N n, ∼ 5 • 10 49 ) and the Sun (N e,⊙ ∼ 10 57 , N n,⊙ = N e,⊙ /4) as point sources of electrons and neutrons, and the Earth (N e,⊕ ≈ N n,⊕ ∼ 4 × 10 51 ), the Milky Way (N e,MW ≈ N n,MW ∼ 10 67 ), and the cosmological matter (N e,cos ∼ 10 79 , N n,cos ∼ 10 78 ) as continuous distributions.We defer to Ref. [7] for a detailed calculation of Eq. (2.8), but adopt two differences introduced by Ref. [80].First, unlike Ref. [7], which studied extragalactic neutrinos and so averaged the contribution of cosmological matter over redshift, here we consider V cos αβ to be only the contribution from the local Universe, i.e., we evaluate Eq. (A8) in Ref. [7] at redshift z = 0. Second, unlike Ref. [7], which only computed the potential sourced by electrons under L e −L µ and L e −L τ , here we compute also the potential sourced by neutrons under L µ − L τ .

Neutrino oscillation probabilities under long-range interactions
We consider mixing between the three active neutrinos, ν e , ν µ , and ν τ .Under the L α − L β symmetry, the Hamiltonian that drives neutrino propagation, in the flavor basis, is (2.9) The first two terms on the right-hand side induce standard oscillations, including SM matter effects; the third one, oscillations due to the new interactions.
In vacuum, oscillations are driven by where E is the neutrino energy, ∆m 2 ij ≡ m 2 i − m 2 j are the mass-squared splittings between neutrino mass eigenstates, and U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix, parametrized [81] in terms of the mixing angles θ 12 , θ 23 , and θ 13 , and the CP-violating phase, δ CP .In the main text, we show results assuming normal neutrino mass ordering (NMO), where m 1 < m 2 < m 3 ; Table 1 shows the values of the mixing parameters that we use, taken from Ref. [82].Appendix C contains results obtained under the inverted mass ordering (IMO), where m 3 < m 2 < m 1 .
In Eq. (2.9), the contribution of SM coherent forward scattering on electrons, mediated by the W boson, is where 6Y e [ρ/(10 14 g cm −3 )] eV is the charged-current neutrinoelectron interaction potential, G F is the Fermi constant, n e is the electron number density, Y e ≡ n e /(n p +n n ) is the electron fraction, i.e., its abundance relative to that of protons and neutrons, n p and n n , and ρ is the matter density.In our work, this contribution is relevant only inside Earth, where matter densities are high.We take ρ to be the average density of underground matter along the trajectory from source to detector, calculated using the Preliminary Reference Earth Model [83]: 2.848 g cm −3 for DUNE and 2.8 g cm −3 for T2HK.The potential above is for neutrinos; for antineutrinos, it flips sign, i.e., V mat → −V mat .
Finally, in Eq. (2.9) the contribution from the new matter interaction is where the potential, V αβ , Eq. (2.8), depends on the mediator mass, m ′ αβ , and coupling, G ′ αβ .The potential above is for neutrinos; for antineutrinos, it flips sign, i.e., V αβ → −V αβ .
The ν α → ν β transition probability associated to the Hamiltonian, Eq. (2.9), is where L is the distance traveled by the neutrino from production to detection, ∆ m2 ij ≡ m2 i − m2 j , with m2 i /2E the eigenvalues of the Hamiltonian, modified from those of H vac by matter effects, and U ′ is the unitary matrix that diagonalizes the Hamiltonian.We parametrize U ′ with the same shape as the PMNS matrix, but evaluated at mixing parameters θ m 12 , θ m 23 , θ m 13 , and δ m CP modified by matter effects.In our work, we compute the oscillation probability, Eq. (2.13), exactly and numerically to arbitrary precision; see Refs.[56,58,[84][85][86][87][88][89][90] for approximate analytical solutions.
For the new matter interactions to affect the oscillation probability, the new matter potential must be at least comparable to the standard contributions in Eq. (2.9), i.e., in vacuum, V αβ ≳ (∆m 2  31 /2E) [inside the Earth, this is instead V αβ ≳ max ∆m 2 31 /2E, V CC ].In DUNE and T2HK, where the first oscillation maxima occur at 2.6 GeV and 0.6 GeV, respectively, this implies that they become important for V αβ ≳ 10 −13 eV.This sets the scale of the potential to which our analysis is sensitive.Later, in Section 3.2, we show how the new interactions affect the probabilities in DUNE and T2HK.
Pioneering studies in Refs.[2] and [3] identified the potential of neutrino oscillations to test new long-range interactions, possibly more stringently than gravitational probes.They focused on interactions with a range equal to the Earth-Sun distance [though ignoring the Yukawa suppression in Eq. (2.5)] and sourced by solar electrons.Reference [2] used Super-Kamiokande atmospheric neutrino data to find g ′ eµ < 8.32 × 10 −26 and g ′ eτ < 8.97 × 10 −26 , at 90% confidence level (C.L.) Reference [3] used solar and reactor neutrino data from KamLAND to find g ′ eµ < 2.06 × 10 −26 and g ′ eτ < 1.77 × 10 −26 , at 3σ, assuming θ 13 = 0 • .Reference [54] studied the effect of the L µ − L τ symmetry via kinetic mixing (see Section 2.1) on ν µ in the long-baseline experiment MINOS.By comparing the potential V µτ sourced by a neutron in the Sun to the fifth-force gravitational potential sourced by it, and applying upper limits on the strength of the latter from torsion-balance experiments [75,76], Ref. [54] set an upper limit on the mixing strength of (ξ − sin θ W χ) < 5 × 10 −24 at 95% C.L. for a long-range interaction with range equal to the Earth-Sun distance.(This is the limit that we saturate when computing the V µτ potential, Eq. (2.6); see Section 2.1.)This translates into an upper limit of g ′ µτ ≤ 2.51 × 10 −26 .For an interaction with a range of the size of the Earth, the upper limit degrades to g ′ µτ ≤ 10 −24 .Reference [8] showed that upper limits on the coefficients that parametrize the strength of non-standard neutrino interactions (NSI) can be translated into upper limits on the coupling strength of flavor-dependent long-range interactions.Figure 1 shows the resulting limits, based on the NSI limits from Refs.[4][5][6].
Recently, Ref. [1] performed a global oscillation analysis of new U (1) symmetries, including L e − L µ and L e − L τ , by using the same experimental data sets used in Nu-FIT 5.0 [92,93].Unlike our analysis, Ref. [1] computed the changing long-range matter potential due to underground matter in the Earth along the trajectory of the neutrinos.Their procedure is more detailed than ours for mediators lighter than 10 −14 eV.However, they explore masses only as low as 10 −18 eV, i.e., an interaction range of 1 A.U.
Reference [7] first showed that the flavor composition of TeV-PeV astrophysical neutrinos, i.e., the relative number of ν e , ν µ , and ν τ , can be used to probe long-range interactions under L e − L µ and L e − L τ sourced by the same collections of nearby and distant electrons that we consider here.Reference [80] refined the statistical methods and included also L µ −L τ .The main effect is that, if the potential sourced by electrons or neutrons were to be  1: Values of the standard mixing parameters used in our analysis.We assume normal neutrino mass ordering (NMO) in the main text.The benchmark values are the best-fit values from Ref. [82].For each parameter over which we minimize our test statistic (see Section 4.1), the minimum is searched for within the range shown, which is the 3σ allowed range from Ref. [82].
We assume no correlation between the parameters.Table C1 shows the parameter ranges that we use in Appendix C to obtain results under the inverted mass ordering (IMO) instead.
dominant, oscillations would turn off, and the flavor composition emitted by the astrophysical sources and received at Earth would be the same; see also Ref. [94]. Figure 1 shows the projected upper limits obtained in Ref. [7] based on estimates of flavor-composition measurements in the envisioned IceCube-Gen2 neutrino telescope [95].
Finally, following Ref.[8], Fig. 1 includes two indirect limits.First, Ref. [9] excluded three mediator mass windows ("Black-hole superradiance") by considering the superradiant growth rate of a gravitationally bound accumulation of light vector bosons around selected stellar-mass and supermassive black holes.Second, Ref. [10] placed a tentative lower limit on the coupling ("Weak gravity conjecture") by studying low-energy effective theories that contain gravity and U (1) gauge fields where at least one particle charged under U (1) is essential for gravity to be the weakest force.
In Fig. 1, we show existing limits as they were published in their original references.Hence, they do not extend to mediators lighter than 10 −14 -10 −20 eV, depending on the limit (except for the proof-of-principle sensitivity based on projected IceCube-Gen2 measurements of the flavor composition [7]).These limits could be recomputed and extended to span lighter mediators, using the same long-range matter potential that we have used, Eq. (2.8), though doing so lies beyond the scope of this work.
3 Long-range interactions in DUNE and T2HK
In our forecasts, we focus on DUNE and T2HK.Below, we overview their features.For each one, we compute appearance and disappearance oscillation probabilities and event rates in neutrino (ν) and antineutrino (ν) beam modes: Appearance, ν mode: This is sensitive mainly to ν µ → ν e transitions.The beam works in neutrino mode, and the detector targets ν e -initiated events.
Appearance, ν mode: This is sensitive mainly to νµ → νe transitions.The beam works in antineutrino mode, and the detector targets νe -initiated events.
Disappearance, ν mode: This is sensitive mainly to ν µ → ν µ survival.The beam works in neutrino mode, and the detector targets ν µ -initiated events.
Disappearance, ν mode: This is sensitive mainly to νµ → νµ survival.The beam works in antineutrino mode, and the detector targets νµ -initiated events.
DUNE will also detect ν τ with energies larger than 3.4 GeV via their charged-current interactions, which allows for interesting physics opportunities [111][112][113][114][115].However, in our analysis, we focus on ν e appearance only and treat ν τ appearance as background; see below.Table 1 shows the values and allowed ranges of the mixing parameters that we use in our analysis, taken from the global oscillation fit of Ref. [82].In the main text, we show results assuming that the true neutrino mass ordering is normal, since there is currently weak preference for it [82,92,93,116].However, as part of our statistical analysis in Section 4, we report sensitivity after minimizing over the mass ordering.Appendix C contains results assuming instead that the true mass ordering is inverted.

DUNE
DUNE will consist of a near detector, about 600 m downstream of the neutrino production point on the Fermilab site, and a far detector, 1285 km away and about 1.5 km underground, in the Sanford Underground Research Facility in South Dakota [117].The near detector will monitor and characterize the neutrino beam (though it has physics capabilities itself, too [118,119]).We focus on the far detector since it offers prime sensitivity to neutrino oscillations.It is a state-of-the-art liquid-argon time projection chamber with a net volume of 40 kton; to generate our results, we consider single-phase detection only [107].Neutrino detection is via charged-current neutrino-argon interaction.Detector deployment will be phased [108], but in our simulations we consider only the final, total detector volume.
DUNE will use the Long Baseline Neutrino Facility (LBNF) neutrino beam produced at Fermilab.There, the Main Injector of the LBNF fires a 1.2-MW beam of protons of 120 GeV onto a graphite target, producing charged mesons that decay in flight to neutrinos.The resulting neutrino flux is wide-band, ranges from a few hundreds of MeV to a few tens of GeV, and is expected to peak at 2.5 GeV, with most neutrinos in the 1-5 GeV range.By changing the polarity of the focusing horns [120,121], the experiment can run in neutrino or antineutrino mode.Following the DUNE Technical Design Report [107], we adopt a run time of 5 years in neutrino mode and 5 years in antineutrino mode.This amounts to 1.1 × 10 21 protons-on-target per year and a net exposure of 480 kton MW year.To produce our results, we use the DUNE simulation configuration from Ref. [17].    1. See Section 3.2 for details.

T2HK
T2HK will consist of near detectors, about 280 m downstream from the neutrino production point at the Japan Proton Accelerator Research Complex (JPARC), and a far detector, 295 km away and about 1.7 km underground, in the Tochibora mines of Japan, 8 km from Super-Kamiokande [45].Like in DUNE, the near detectors will monitor and characterize the neutrino beam, and we focus on the far detector.It will be a tank filled with purified water, with a net volume of 187 kton, whose internal wall is lined with photomultipliers (PMTs).Neutrino detection is via quasielastic charged-current scattering (QECC), i.e., ν l + n → p + l − and νl + p → n + l + (l = e, µ, τ ), and via charged-current deep inelastic scattering (DIS), i.e., ν l + N → l − + X and νl + N → l + + X (l = e, µ), where X represents final-state hadrons (ν τ DIS is suppressed due to the large tauon mass).Electrons emit gamma rays by bremsstrahlung and e + e − annihilation, which register as a fuzzy ring on the PMTs.Muons emit Cherenkov light, which registers as a sharply defined ring.Like its predecessor, T2K (Tokai-to-Kamioka) [102], T2HK will use the 2.5 • -off-axis JPARC neutrino beam [122].To produce it, JPARC fires a 1.3-MW beam of protons of 30 GeV onto a graphite target.The resulting neutrino flux is narrow-band, ranges from a few MeV to a few GeV, and is expected to peak at 600 MeV, with most neutrinos in the 100-3000 MeV range.As in DUNE, by changing the polarity of the focusing horns, T2HK can run in neutrino or antineutrino mode [123].Following Ref. [18], we adopt a run time of 2.5 years in neutrino mode and 7.5 years in antineutrino mode, in accordance with the default 1:3 ratio planned for them.This amounts to 2.7 × 10 22 protons-on-target per year and a net exposure of 2431 kton MW year.To produce our results, we match the binned event spectra that we generate under standard oscillations with those of Ref. [18].

Oscillation probabilities
Appendix A, especially Fig. A1 therein, shows in detail the effects of long-range interactions on the modified mixing parameters θ m 12 , θ m 23 , and θ m 13 ; here, we summarize them.Differences in their behavior under the different symmetries stem from differences in the flavor structure of the new matter potential, V αβ in Eq. (2.9).
The solar angle in matter rapidly approximates its maximum value of θ m 12 = 90 • already at a few GeV, for all symmetries.(We use this later, in Section 4.1, to justify why we neglect the effect on our forecasts of the uncertainty in its value in vacuum, θ 12 .)For DUNE and T2HK, the mixing angles that drive the probabilities are the atmospheric angle, θ m 23 , and the reactor angle, θ m 13 .Assuming θ m 12 = 90 • , Ref. [58] showed that the transition probabilities for ν µ → ν e and νµ → νe are ∝ sin 2 θ m 23 sin 2 θ m 13 and the survival probabilities for ν µ → ν µ and νµ → νµ are ∝ sin 2 2θ m 23 , with a more nuanced dependence on θ m 13 .The deviation of θ m 23 relative to θ 23 grows with energy, though there are differences depending on which symmetry is active: θ m 23 grows under L e − L τ and L µ − L τ , and shrinks under L e − L µ .The reactor angle in matter, θ m 13 , grows appreciably with energy under L e − L µ and L e − L τ , and falls to about 0 • under L µ − L τ .
Figure 3 shows the oscillation probabilities, Eq. (2.13), computed under the three symmetries in each of the four detection channels listed in Section 3.1, for DUNE and T2HK.To illustrate the effects of long-range interactions, we pick a relatively high value of the potential, V αβ = 1.3 × 10 −13 eV; later, when producing our results, we vary this value.Via Eq. ( of long-range interactions with underground matter on the probabilities in the former are more prominent than in the latter [90].The effects are more clearly visible in the transition probabilities: the oscillation maxima shift to lower energies, due to a change in the effective mass-splitting ∆m 2 31,m , in agreement with Ref. [56], and the oscillation amplitudes grow, especially after the first maximum.The effects are more prominent under L e − L τ because θ m 23 and θ m 13 are enhanced, whereas under L e − L µ and L µ − L τ only one of them is; see Appendix A for details.Naturally, for weaker potentials, the above effects are lessened.

Event rates
We compute event rates in DUNE and T2HK using GLoBES [124,125], extended with the snu matrix-diagonalization library [126,127], by modeling their technical design specifications [17,19] of efficiency, operation times, and backgrounds.Because we are interested in assessing the mean sensitivity of the experiments (Section 4.1), we compute only mean event rates and do not generate event spectra that include fluctuations from the mean rates.We bin event rates in reconstructed energy, E rec , built from the detected secondaries born in neutrino interactions.In both experiments, because the far detectors cannot distinguish between neutrinos and antineutrinos, there is irreducible contamination from "wrong-sign" events; we add it to the signal.
DUNE.-We consider events with E rec in the range 0-110 GeV, with 64 bins within 0-8 GeV, each 0.125 GeV wide, and 16 bins within 8-110 GeV, of varying widths.In the appearance channel, the signal is due to the charged-current (CC) interactions of ν e , in neutrino mode, and of νe , in antineutrino mode.The background consists of (i) the CC interactions of "intrinsic" ν e and νe , i.e., those created as such that survive the flight to the detector (from ν e → ν e and νe → νe ); (ii) the CC interactions of ν µ and νµ whose final-state muons are misidentified as electrons (from ν µ → ν µ and νµ → νµ ); (iii) the   1.See Section 3.3 for details.
CC interactions of ν τ and ντ (from ν µ → ν τ and νµ → ντ ); and (iv) the neutral-current (NC) interactions of neutrinos of all flavors.In the disappearance channel, the signal is due to the CC interactions of ν µ , in neutrino mode, and of νµ , in antineutrino mode.The background consists of (i) the CC interactions of ν τ and ντ (from ν µ → ν τ and νµ → ντ ); and (ii) the NC interactions of neutrinos of all flavors.
T2HK.-We consider events with E rec in the range 0.1-3 GeV, with 29 bins, each 0.1 GeV wide.In the appearance channel, the signal is due to the CC interactions of ν e , in neutrino mode, and of νe , in antineutrino mode.The background consists of (i) the CC interactions of intrinsic ν e and νe (from ν e → ν e and νe → νe ); (ii) the CC interactions of ν µ and νµ whose final-state muons produce fuzzy Cherenkov rings; and (iii) the NC interactions of neutrinos of all flavors.In the disappearance channel, the signal is due to the CC interactions of ν µ , in neutrino mode, and of νµ , in antineutrino mode.The background consists of (i) the CC interactions of intrinsic ν e and νe (from ν e → ν e and νe → νe ); and (ii) the NC interactions of neutrinos of all flavors.
Figure 4 shows the mean event-rate spectra under long-range interactions for each detection channel in DUNE and T2HK, including all the above backgrounds, and computed using the same illustrative value of the long-range potential V αβ as in Fig. 3.The event rates in T2HK are higher than in DUNE due to its larger size.The shapes of the event spectra in Fig. 4 reflect those of the oscillation probabilities in Fig. 3. Long-range interactions affect each detection channel differently, but there are common features among them.Broadly stated, in the appearance channels, they enhance the event rates relative to the standardoscillations rates (with the exception of L e − L µ in neutrino mode for DUNE).For DUNE, additionally, they slightly shift the event rates to lower energies, reflecting the shift in the oscillation maxima.In the disappearance channels, the effect of long-range interactions is more nuanced; the event rate is enhanced or reduced depending on the symmetry and the energy.The above features in the event spectra hold for other values of the potential, though, naturally, their prominence varies depending on the value.
Table 2 shows the mean expected number of signal and background events for each detection channel, assuming standard oscillations.In all channels, the signal is dominant.In T2HK, unlike DUNE, neutrino and antineutrino event rates are comparable, due to the 1:3 ratio between run times in neutrino and antineutrino modes that compensates for the smaller antineutrino cross sections.In DUNE, neutrino event rates are higher than T2HK due to its longer run time in neutrino mode.These general features of the event rates hold also in the presence of long-range interactions.
Below, we show how the above features grant DUNE and T2HK the capability to probe long-range interactions, and how they organically complement each other.
4 Projected constraints and discovery potential

Statistical methods
We forecast the capability of DUNE and T2HK to probe long-range interactions that stems from the modification of the oscillation probabilities (Section 3.2), based on the detailed computation of event rates outlined above (Section 3.3).Our forecasts are two-fold: we forecast constraints on long-range interactions -on the long-range matter potential and ultimately on the mediator mass and coupling -assuming that no evidence for them is found, and we forecast prospects of discovering them and measuring their parameter values.
We study each symmetry, L e − L µ , L e − L τ , and L µ − L τ , separately.For a given symmetry, we generate two sets of event spectra, including signal plus backgrounds, for each of the four detection channels of T2HK and DUNE (Section 3.3): a "true" spectrum, which we take to be the observed spectrum, and a set of "test" spectra, generated for test values of the parameters, that we compare against it.When forecasting constraints, in Section 4.2, we compute the true spectrum fixing the true value of the potential to be V true αβ = 0, which corresponds to standard oscillations.When forecasting discovery prospects, in Section 4.3, we compute the true spectrum fixing V true αβ to a specific nonzero choice.We expand on this below.
To compare true and test event spectra, we follow Refs.[128][129][130] and adopt a Poissonian χ 2 function.For each experiment e = {T2HK, DUNE}, and for each detection channel c = {app ν, app ν, disapp ν, disapp ν}, this is where N true e,c,i and N test e,c,i are the true and test event rates in the i-th bin of E rec , N e is the number of bins of E rec (Section 3.
where N s,true e,c,i and N b,true e,c,i are, respectively, the number of signal (s) and background (b) events, summed over all channels, computed using the true values of the mixing parameters, mass ordering, and potential.The test number of events is 3) where π s e,c and π b e,c,k are normalization errors on the signal and background rates, which lie between 2% and 10%; see Table D1 in Appendix D for their values, taken from Refs.[17,19].The background rates do not vary significantly upon changing the mass ordering.
For T2HK or DUNE, separately or together, we compute the total χ 2 by adding the contributions of all the detection channels, i.e., and we take the contributions of different channels to be uncorrelated.We report sensitivity by comparing the minimum value of the χ 2 function, χ 2 e,min , reached when it is evaluated at the true values of the parameters, V true αβ , θ true , and o true , against the value of χ 2 evaluated at test values of the parameters.We treat θ and o as nuance parameters and profile over them (more on this below).For instance, for DUNE, and similarly for T2HK and DUNE+T2HK.In the main text, we fix θ true to its presentday best-fit value under NMO, from Ref. [82] (Table 1) and o true also to NMO; we fix them to inverted mass ordering in Appendix C. When reporting constraints on long-range interactions, we set V true αβ = 0 and extract from ∆χ 2 DUNE , ∆χ 2 T2HK , and ∆χ 2 DUNE+T2HK the upper limits on the inferred value of V αβ , for 1 degree of freedom (d.o.f).When reporting discovery potential, we fix V true αβ to a nonzero illustrative value and report the inferred range of values of V αβ , again for 1 d.o.f.When reporting discovery, we also study the correlation between V αβ and δ CP or sin 2 θ 23 .In those cases, we use instead, respectively, ∆χ When profiling, we minimize the ∆χ 2 functions above with respect to sin 2 θ 23 , δ CP , and |∆m 2  31 | by varying them within their 3σ allowed ranges from Ref. [82]; see Table 1.We do not include correlations between them, since these are expected to weaken in coming years (see, e.g., Ref. [131]), nor do we include pull terms on the mixing parameters in the test-statistic, in order to be conservative.We keep θ 13 and θ 12 fixed at their present-day best-fit values [82].For θ 13 , the precision that Daya Bay has achieved, of 2.8% [132], is not expected to be improved upon by upcoming experiments.For θ 12 , whose presentday uncertainty is of 4.5%, we expect only weak sensitivity in the oscillation probabilities (see Section 3.2), so fixing its value is a justified approximation.For the neutrino mass ordering, we adopt a simplified approach where switching from NMO to IMO amounts only to flipping the sign of ∆m 2  31 to make it negative.This approach is motivated by the fact that the present-day 3σ allowed ranges of the mixing parameters are similar in the NMO V eτ , and V µτ , using DUNE, T2HK, and their combination.The ∆χ 2 function is Eq.(4.7) and similar ones, assuming V true αβ = 0 for the true value of the potentials.We profile over the values of the most relevant standard mixing parameters and over the neutrino mass ordering; see Table 1.See Table 3 for the resulting upper limits on the potentials and Fig. 6 for the corresponding constraints on the mass and coupling of the new mediator.Combining DUNE and T2HK not only provides sensitivity to lower values of the potential, but also removes degeneracies in the test-statistics that would otherwise weaken the sensitivity.See Sections 4.1 and 4.2 for details.
and IMO [82,92,93], except for δ CP and ∆m 2  31 .(Further, in the next decade, DUNE is expected to determine the mass ordering, though, admittedly, non-standard oscillations like those induced by long-range interactions may confound this [56].) Below, we compute the test-statistics by varying the long-range matter potential in the range 10 −15 ≤ V αβ /eV ≤ 3 × 10 −13 , where its effects are potentially visible in DUNE and T2HK.This range is wide enough to comfortably place constraints or make measurements with significant statistical confidence.

Projected constraints on long-range interactions
Figure 5 shows how the test-statistics for constraints, e.g., Eq. (4.7), vary with the potential, for the three symmetries and for DUNE, T2HK, and their combination.As expected, they are smallest closer to the true value of the potential, V true αβ = 0, and grow as they move away from it.At high values of V αβ , the test-statistics for DUNE and T2HK dip, reflecting a loss of sensitivity due to V αβ being degenerate with θ 23 and δ CP .Appendix B expands on this.Combining DUNE and T2HK removes the dips: T2HK lifts the degeneracies due to θ 23 and δ CP , while DUNE fixes the mass ordering, i.e., the sign of ∆m 2  31 .Thus, our results reveal novel insight: the interplay of DUNE and T2HK facilitates degeneracy-free constraints on flavor-dependent long-range neutrino interactions.
Table 3 shows the resulting upper limits on the potential.They are strongest for V µτ , followed by V eτ and then V eµ .For DUNE, the limits are driven predominantly by the runs in neutrino mode, which contribute most of the total event rate; see Table 2 and Fig. 4. For T2HK, the runs in neutrino and antineutrino modes contribute comparably.
For L µ − L τ , the limits on V µτ are strongest because long-range interactions affect mainly the disappearance probabilities, ν µ → ν µ and νµ → νµ (see Fig. 3), whose associated disappearance detection channels have high event rates (see Fig. 4), making deviations from standard oscillations easier to spot.For L e − L τ , long-range interactions enhance instead the appearance probabilities, ν µ → ν e and νµ → νe , but the appearance detection channels have lower rates, so the limits on V eτ are weaker.For L e −L µ , long-range interactions affect both the appearance probabilities -though less so than under the other two symmetries -and the disappearance probabilities -though less so than under the L µ −L τ symmetry; as a result, the limits on V eµ are the weakest.Figure 6 (also Fig. 1) shows the corresponding upper limits on G ′ αβ (refer to Eq. (2.7)) for varying m ′ αβ , translated from the upper limits on V αβ in Table 3 via the definition of the potential, Eq. (2.8).Each curve in Fig. 6 is an isocontour of potential that saturates each of the upper limits in Table 3.The curves show step-like transitions at various values of m ′ αβ : as explained in Ref. [7] (see also Section 2.2), each transition reflects the interaction range becoming long enough for a new source of electrons or neutrons to contribute to the total potential, Eq. (2.8).For m ′ αβ ∼ 10 −18 -10 −10 eV, the Earth and the Moon dominate the upper limits; for m ′ αβ ≲ 10 −18 eV, the Sun dominates; for m ′ αβ ≲ 10 −27 eV, the Milky Way dominates; and, for m ′ αβ ≲ 10 −33 eV, cosmological electrons and neutrinos dominate.
Down to m ′ αβ ∼ 10 −18 eV, where direct limits on flavor-dependent long-range neutrino interactions exist, our projected limits improve on existing ones that use atmospheric neutrinos [2], are comparable to limits that use solar, reactor [3], and accelerator neutrinos [54], and to limits culled from non-standard interactions [4][5][6], but are weaker than limits from a global fit to oscillation data [1].
Below m ′ αβ ∼ 10 −18 eV, our projected limits tread into a largely unexplored range.To our knowledge, the only constraints that exist there, other than the indirect, tentative one from the weak gravity conjecture [10], are from measurements of the flavor composition of high-energy astrophysical neutrinos in the IceCube neutrino telescope, from Ref. [7], which, however, were only produced at the 1σ level as a proof of principle of the sensitivity.A recent recalculation [80] found comparable results at higher statistical significance, using more solid statistical methods.Ideally, the sensitivity that could be reaped from highenergy astrophysical neutrinos is unmatched due to them having energies in the TeV-PeV range -which enhances the possible contribution of long-range interactions relative to

Detector
Upper limits (2σ) on potential [10 Table 3: Projected upper limits (2σ) on the long-range matter potentials V eµ , V eτ , and V µτ , using DUNE, T2HK, and their combination.See Fig. 5 for the test-statistics from whence they originate and Fig. 6 for constraints on the mass and coupling of the associated mediator.See Sections 4.1 and 4.2 for details.αβ , with mass m ′ αβ , that mediates flavor-dependent long-range neutrino interactions, using DUNE, T2HK, and their combination.Same as Fig. 1, but now showing also limits using DUNE and T2HK separately.Left: For neutrino-electron interactions under the L e − L µ symmetry.Center: For neutrino-electron interactions under L e − L τ .Right: For neutrinoneutron interactions under L µ − L τ ; existing limits from accelerator neutrinos in MINOS (95% C.L.) are from Ref. [54], and for non-standard interactions (NSI) we follow Appendix C in Ref. [80].See Section 4.2 for details. standard oscillations -and to the fact that neutrinos of all flavors are detected.However, it is presently downplayed by large astrophysical uncertainties, limited event rates, and the difficulty in measuring the flavor composition in neutrino telescopes.These issues will likely be surmounted in the future [133][134][135].For now, Fig. 6 shows the projected proof-of-principle sensitivity of the envisioned IceCube-Gen2 upgrade [95], from Ref. [7].Our limits from DUNE and T2HK improve on it significantly due to high event rates and well-characterized neutrino beams.

Discovering subdominant long-range interactions
Figure 7 shows the inferred allowed ranges of the coupling, for varying mediator mass, for three illustrative choices of the true value of the long-range potential, one for each symmetry: V true eµ = 2.12 × 10 −14 eV, V true eτ = 1.6 × 10 −14 eV, and V true µτ = 1.12 × 10 −14 eV.They represent long-range interactions that are subdominant to standard oscillations; see Section 2.3.To generate Fig. 7, first we compute the allowed ranges of V αβ using the discovery test-statistics, e.g., Eq. (4.7) with V true αβ fixed to the above illustrative choices, and then we use Eq.(2.8) to translate those into allowed ranges of G ′ αβ for different values of m ′ αβ .We present results for a modest discovery significance of 3σ.
Figure 7 shows that DUNE and T2HK, by themselves, can only place upper limits on G ′ αβ .Thus, our discovery forecasts also reveal novel insight: DUNE and T2HK, a relative measurement uncertainty of 90%-100%.For larger values of the true potential, discovery claims should be stronger and the uncertainty in its measurement should shrink.
Figures 8 and 9 reveal that the reason behind the difficulty of T2HK and DUNE to discover subdominant long-range interactions by themselves are the uncertainties in δ CP , θ 23 , and the neutrino mass ordering.Appendix B shows this in detail; here we summarize.On the one hand, in T2HK, the shorter baseline provides less contamination from fake CP violation induced by SM matter effects and, therefore, higher precision in measuring δ CP [136][137][138][139][140], while the high event rates in the disappearance channels provide high precision on θ 23 .However, at the same time, the shorter baseline reduces the sensitivity to the mass ordering and provides a shorter neutrino travel time during which long-range interactions may act.On the other hand, in DUNE, the longer baseline helps to pin down the mass ordering, but introduces more contamination from fake CP violation, which degrades the sensitivity to δ CP compared to T2HK.Also, in the presence of long-range interactions, DUNE can measure θ 23 significantly less precisely than T2HK.
Thus, combining DUNE and T2HK improves the sensitivity with which δ CP , θ 23 , and the mass ordering can be measured, weakens the degeneracies between them and the longrange potential, and allows for its measurement at a high statistical significance.

Summary and outlook
Neutrinos are powerful probes of new physics.Extant uncertainties in their properties leave room for the conceivable possibility that they experience interactions with matter beyond weak ones.Discovering them would not only further our view of neutrino physics, but also represent striking evidence of physics beyond the Standard Model.In the 2030s, the nextgeneration long-baseline neutrino experiments DUNE and T2HK may provide us with an opportunity to look for new neutrino interactions more incisively than ever before, thanks to high event rates and well-characterized neutrino beams.We have forecast their reach.Because we use detailed simulations of the detectors, including efficiencies, run times, and backgrounds, our predictions are realistic.
Our forecasts are geared at new flavor-dependent neutrino interactions that are introduced by gauging three different accidental global lepton-number symmetries of the Standard Model, generated by L e − L µ , L e − L τ , and L µ − L τ , that have received prior attention in other experimental settings [1-3, 7, 8, 54, 56, 58].We focus on them because they can be gauged anomaly-free, so the only new particle introduced is a neutral vector boson that mediates the interaction.Its mass and coupling strength are a priori undetermined.Gauging L e − L µ and L e − L τ introduces new neutrino-electron interactions.Gauging L µ − L τ introduces new neutrino-neutron interactions.
Under these interactions, electrons and neutrons source a flavor-dependent potential that may affect neutrino oscillations.We concentrate on ultra-light mediators, with masses below 10 −10 eV.They induce interactions whose range is ultra-long -ranging from hundreds of meters to Gpc, depending on the mass -so that neutrinos may experience the potential sourced by a large number of nearby and distant electrons and neutrons in the Earth, the Moon, the Sun, the Milky Way, and the cosmological matter distribution [7].Yet, because the coupling strength may be tiny, their effects on the oscillation probabilities may be subtle and, therefore, testable with future experiments, like DUNE and T2HK.
Our forecasts reveal two novel, promising perspectives.First, while DUNE and T2HK, individually, should be able to improve on present-day upper limits on the coupling strength of the new interaction, their individual sensitivities are hampered by degeneracies due to uncertainties in the mixing angle θ 23 , the CP-violating phase, δ CP , and the neutrino mass ordering.Yet, DUNE and T2HK have complementary capabilities: while T2HK is especially well-suited to measure θ 23 and δ CP , DUNE is especially well-suited to measure the neutrino mass ordering.Thus, combining DUNE and T2HK removes parameter degeneracies, which tightens the upper limits on long-range neutrino interactions.Second, and more importantly, neither DUNE nor T2HK, by itself, may discover subdominant long-range interactions, owing to parameter degeneracies, but their combination may .Thus, our forecasts stress the need for combining measurements in DUNE and T2HK to probe long-range interactions.
More broadly, our results illustrate the known need for complementarity in the future long-baseline neutrino program, not only to measure the standard mixing parameters, but to search for new physics.For flavor-dependent long-range neutrino interactions, a future global fit to oscillation data is poised to deliver substantially improved limits or transformative discovery.

A Effects of long-range interactions on neutrino oscillation parameters
Figures A1 and A2 show the variation with neutrino energy of, respectively, the mixing angles and mass splittings under the new matter interactions for the three symmetries, and geared at DUNE for illustration.Their behavior when geared at T2HK is similar.
Regarding the mixing angles, their behavior in Fig. A1 backs the explanation of the behavior of the oscillation probabilities in Section 3.2, in agreement with Refs.[56,58,90].
Regarding the mass splittings, the energy at which the first oscillation maximum occurs in the probabilities (Fig. 3) is determined mostly by ∆m 2  31,m .Figure A2 shows that, under L e −L µ and L e −L τ , ∆m 2  31,m evolves similarly with energy; this explains why the oscillation maxima for these two symmetries occur at approximately the same energy.The change with energy of ∆m 2  21,m affects the ν µ → ν µ and ν µ → ν τ probabilities, but only for baselines of around 10000 km, as shown in Ref. [89].Table B1: Effect of profiling over the different parameters on the test-statistic for constraints.The test-statistic is computed using Eq.(4.7), assuming V true αβ = 0.In this table, we report the values of the test-statistic calculated at an illustrative choice of the potential, V αβ = 3 × 10 −14 eV.We show effects separately for DUNE (D), T2HK (T), and their combination (D+T).Our main results are obtained by profiling over δ CP , sin 2 θ 23 , |∆m 2  31 |, and the neutrino mass ordering, o, i.e., the sign of ∆m 2  31 , in our prescription (see Section 4.1).Results for partial profiling are shown only for the purpose of singling out the effect of different parameters.The values of the parameters that are not profiled over are fixed to the benchmark values in Table 1.

B Effects of oscillation parameter uncertainties on constraints
Table B1 illustrates how the uncertainty on different oscillation parameters affect the constraint the test-statistic, ∆χ 2 , i.e., Eq. (4.7) and similar ones, via profiling over them.The table reports the test-statistic computed assuming that the true value of the potential is V true αβ = 0, and evaluated at an illustrative test value of V αβ = 3 × 10 −14 eV.The observations we make regarding Table B1 hold also for other test values of V αβ .Our main results are computed by profiling the test-statistic over δ CP , sin 2 θ 23 , |∆m 2  31 |, and o.Table B1 shows what the effect is of profiling only over one or two of them at a time.
Because DUNE can measure |∆m 2  31 | and the mass ordering precisely, profiling over their values has little effect compared to keeping them fixed.Similarly, because T2HK can measure θ 23 and δ CP precisely, profiling over their values has little effect compared to keeping them fixed.Combining DUNE and T2HK affords both.
Profiling over sin 2 θ 23 and δ CP has the largest effect on the test-statistic.The shrinking of ∆χ 2 by profiling over sin 2 θ 23 , compared to keeping it fixed, comes via the ν µ → ν µ and νµ → νµ disappearance probabilities, which are ∝ sin 2 θ 23 (see, e.g., Eq. ( 33) in Ref. [87]).The shrinking of ∆χ 2 by profiling over δ CP comes instead via the ν µ → ν e and νµ → νe appearance probabilities, which depend on δ CP .Thus, under L e − L µ , the test-statistic is driven by the uncertainty in sin 2 θ 23 and δ CP -via profiling over them; under L e − L τ , it is driven mostly by δ CP ; and, under L µ − L τ , it is driven mostly by sin 2 θ 23 .

C Constraints assuming inverted mass ordering
In the main text, we showed results obtained assuming that the true neutrino mass ordering is normal.Here we show constraints assuming instead that the true ordering is inverted.They are broadly similar to the ones obtained under normal ordering (Section 4.2).
Table C1 shows the values and allowed ranges of the mixing parameters that we use when assuming that the inverted mass ordering is true, taken from the global oscillation fit of Ref. [82].The values are similar as for normal ordering (Table 3), with the important difference that for inverted ordering the benchmark value of θ 23 lies in the higher octant, which has consequences for the sensitivity, as we point out below.
Figure 5 shows how the test-statistics for constraints, e.g., Eq. (4.7), vary with the longrange potential, for the three symmetries and for DUNE, T2HK, and their combination, assuming the true mass ordering is inverted.Their behavior is broadly similar to that in Fig. 6 for normal ordering.However, in Fig. C1 the weakening of the sensitivity due to parameter degeneracies is milder than in Fig. 5 because, assuming inverted mass ordering, θ 23 lies in the higher octant, which lessens the influence of δ CP [136,140].
Table C2 shows the resulting upper limits on the long-range potentials.Like for normal ordering (Table 3), they are strongest for V µτ .Compared to normal ordering, the limits on V eµ are stronger by a factor of 2 or more, and the limits on V eτ are similarly weaker.-26 - V eµ , V eτ , and V µτ , using DUNE, T2HK, and their combination, assuming that the true neutrino mass ordering is inverted.Same as Fig. 5, made assuming that the true mass ordering is normal, but instead assuming that it is inverted.See Table C2 for the resulting upper limits on the potentials.See Sections 4.1 and 4.2 in the main text, and Appendix C for details.

Detector
Upper limit (2σ) on potential [10  C2: Projected upper limits (2σ) on the long-range matter potentials V eµ , V eτ , and V µτ , using DUNE, T2HK, and their combination, assuming that the true neutrino mass ordering is inverted.Same as Table 3, made assuming that the true mass ordering is normal, but instead assuming that it is inverted.See Fig. C1 for the test-statistics from whence they originate and Fig. Table D1: Normalization errors of the event rates associated to the signal and background detection channels in DUNE and T2HK.They are shown separately for the neutrino (ν) and antineutrino (ν) modes, and for the appearance ("App.")and disappearance ("Disapp.")channels.Normalization errors are used to compute test event spectra, Eq. (4.3).The values are taken from Ref. [17,18].See Section 4.1 for details.

1 Figure 2 :
Feynman diagrams for neutrino-matter interactions.Each diagram corresponds to a term in the Lagrangian, Eq. (2.1): (a) SM contribution mediated by the Z boson, (b) new contribution from the gauge symmetry U (1) Lα−L β , mediated by the new boson Z ′ αβ , and (c) mixing between Z and Z ′ αβ .In our analysis, (a) is significant only for neutrinos inside the Earth.For L e −L β symmetries, (b) is the only additional contribution sourced by electrons.For L µ − L τ , (c) is instead the only additional contribution sourced by neutrons.See Section 2.1 for details.

Figure 3 :
Figure 3: Neutrino oscillation probabilities for T2HK (left column) and DUNE (right column) under flavor-dependent long-range neutrino interactions.The interactions are induced by the lepton-number symmetry L e − L µ , L e − L τ , or L µ − L τ .For this figure, we fix the long-range potential to V αβ = 1.3 × 10 −13 eV as illustration, and the standard mixing parameters to their benchmark values from Table1.See Section 3.2 for details.

Figure 4 :
Figure 4: Expected mean number of detected events in T2HK (left column) and DUNE (right column) under flavor-dependent long-range neutrino interactions.The interactions are induced by the symmetry L e − L µ , L e − L τ , or L µ − L τ .For T2HK, we use 2.5 years in ν mode and 7.5 years in ν mode.For DUNE, we use 5 years in ν mode and 5 years in ν mode.For this figure, we fix the long-range potential to V αβ = 1.3 × 10 −13 eV as illustration, and the standard mixing parameters to their benchmark values from Table1.See Section 3.3 for details.
3), θ ≡ {sin 2 θ 23 , δ CP , |∆m 2 31 |} are the test values of the most relevant mixing parameters (more on this later), o = {NMO, IMO} is the test mass ordering, and ξ s and ξ b,c,k are, respectively, pull terms for the systematic uncertainties on the signal and the k-th background contribution to detection channel c, from the list of contributions in Section 3.3.The pull terms have the same values in neutrino and antineutrino mode, and are uncorrelated with one another.The true number of events is N true e,c,i = N s,true e,c,i + N b,true e,c,i ,

Figure 5 :
Figure5: Projected test-statistic used to constrain the long-range matter potentials V eµ , V eτ , and V µτ , using DUNE, T2HK, and their combination.The ∆χ 2 function is Eq.(4.7) and similar ones, assuming V true αβ = 0 for the true value of the potentials.We profile over the values of the most relevant standard mixing parameters and over the neutrino mass ordering; see Table1.See Table3for the resulting upper limits on the potentials and Fig.6for the corresponding constraints on the mass and coupling of the new mediator.Combining DUNE and T2HK not only provides sensitivity to lower values of the potential, but also removes degeneracies in the test-statistics that would otherwise weaken the sensitivity.See Sections 4.1 and 4.2 for details.

10 − 35 Figure 6 :
Figure 6: Projected upper limits on the effective coupling, G ′ αβ (Eq.(2.7)), of the new boson, Z ′αβ , with mass m ′ αβ , that mediates flavor-dependent long-range neutrino interactions, using DUNE, T2HK, and their combination.Same as Fig.1, but now showing also limits using DUNE and T2HK separately.Left: For neutrino-electron interactions under the L e − L µ symmetry.Center: For neutrino-electron interactions under L e − L τ .Right: For neutrinoneutron interactions under L µ − L τ ; existing limits from accelerator neutrinos in MINOS (95% C.L.) are from Ref.[54], and for non-standard interactions (NSI) we follow Appendix C in Ref.[80].See Section 4.2 for details.

Figure A1 :
FigureA1: Variation of the effective neutrino mixing angles with neutrino energy, for the three lepton-number symmetries, L α − L β .For this plot, we adopt the baseline, average matter potential, and approximate energy range of DUNE; see Section 3.1.1.For all the symmetries, we adopt an illustrative value of the new matter potential of V αβ = 1.3 × 10 −13 eV.For comparison, we include results using only standard matter effects and in a vacuum.The values of the mixing parameters in vacuum are from Table1, except with sin 2 θ 23 = 0.5.See Appendix A for details.

2 ]Figure A2 :
Figure A2: Variation of the effective neutrino mass splittings with neutrino energy, for the three lepton-number symmetries, L α − L β .Same as Fig. A1, but for the mass splittings.See Appendix A for details.

Figure C1 :
FigureC1: Projected test-statistic used to constrain the long-range matter potentials V eµ , V eτ , and V µτ , using DUNE, T2HK, and their combination, assuming that the true neutrino mass ordering is inverted.Same as Fig.5, made assuming that the true mass ordering is normal, but instead assuming that it is inverted.See TableC2for the resulting upper limits on the potentials.See Sections 4.1 and 4.2 in the main text, and Appendix C for details.

Table 2 :
[82], multiple combinations of m ′ αβ and G ′ αβ can yield this value of the potential, or any other.Because the baseline for DUNE is longer than for T2HK, the effects Mean number of signal and background events, summed over all background channels, expected in DUNE and T2HK after their full run times.For this table, whose aim is illustrative only, we assume standard oscillations and normal mass ordering (NMO).DUNE runs for 5 years in ν mode and 5 years in ν mode.T2HK runs for 2.5 years in ν mode and 7.5 years in ν mode.To compute the rates in this table, we fix the values of the standard mixing parameters to their benchmark values from Ref.[82]; see Table1.In the main text, to produce results, we also compute event rates in the presence of the new long-range neutrino interactions (not shown in this table).In those cases, the relative sizes of the event rates in the different detection channels are roughly as in this table.See Section 3.3 for details.
this plot, we adopt the baseline, average matter potential, and approximate energy range of DUNE; see Section 3.1.1.For all the symmetries, we adopt an illustrative value of the new matter potential of V αβ = 1.3 × 10 −13 eV.For comparison, we include results using only standard matter effects and in a vacuum.The values of the mixing parameters in vacuum are from Table1, except with sin 2 θ 23 = 0.5.See Appendix A for details.
the Department of Science and Technology (DST), Govt. of India, and the Research Grant (sanction order no.SB/SJF/2020-21/21) provided by the Science and Engineering Research Board (SERB), Govt. of India, under the Swarnajayanti Fellowship project.S.K.A. would like to thank the United States-India Educational Foundation (USIEF) for providing the financial support through the Fulbright-Nehru Academic and Professional Excellence Fellowship (Award no.2710/F-N APE/2021).M.S. acknowledges the financial support from the DST, Govt. of India (DST/INSPIRE Fellowship/2018/IF180059). M.B. is supported by the Villum Fonden under the project no.29388.The numerical simulations are carried out using the "SAMKHYA: High-Performance Computing Facility" at the Institute of Physics, Bhubaneswar, India.

Table C1 :
[82]es of the standard mixing parameters used in our analysis, assuming that the true neutrino mass ordering is inverted.Same as Table1, made assuming that the true mass ordering is normal, but instead assuming that it is inverted (IMO).The benchmark values are the best-fit values from Ref.[82].