Exploring Flavor-Dependent Long-Range Forces in Long-Baseline Neutrino Oscillation Experiments

The Standard Model gauge group can be extended with minimal matter content by introducing anomaly free U(1) symmetry, such as $L_e-L_{\mu}$ or $L_e-L_{\tau}$. If the neutral gauge boson corresponding to this abelian symmetry is ultra-light, then it will give rise to flavor-dependent long-range leptonic force, which can have significant impact on neutrino oscillations. For an instance, the electrons inside the Sun can generate a flavor-dependent long-range potential at the Earth surface, which can suppress the $\nu_{\mu} \to \nu_e$ appearance probability in terrestrial experiments. The sign of this potential is opposite for anti-neutrinos, and affects the oscillations of (anti-)neutrinos in different fashion. This feature invokes fake CP-asymmetry like the SM matter effect and can severely affect the leptonic CP-violation searches in long-baseline experiments. In this paper, we study in detail the possible impacts of these long-range flavor-diagonal neutral current interactions due to $L_e-L_{\mu}$ symmetry, when (anti-)neutrinos travel from Fermilab to Homestake (1300 km) and CERN to Pyh\"asalmi (2290 km) in the context of future high-precision superbeam facilities, DUNE and LBNO respectively. If there is no signal of long-range force, DUNE (LBNO) can place stringent constraint on the effective gauge coupling $\alpha_{e\mu}<1.9 \times 10^{-53}~(7.8 \times 10^{-54})$ at 90% C.L., which is almost 30 (70) times better than the existing bound from the Super-Kamiokande experiment. We also observe that if $\alpha_{e\mu} \geq 2 \times 10^{-52}$, the CP-violation discovery reach of these future facilities vanishes completely. The mass hierarchy measurement remains robust in DUNE (LBNO) if $\alpha_{e\mu}<5 \times 10^{-52}~(10^{-52})$.


Introduction and Motivation
Over the past two decades or so, active attempts have been made both in experimental and theoretical fronts to improve our knowledge about neutrinos [1][2][3]. The three most important fundamental issues that have taken center stage in neutrino physics as a part of these activities are the following. First issue: how tiny is the neutrino mass? Second issue: can a neutrino turn into its own anti-particle? Third issue: do different neutrino flavors oscillate into one another? To shed light on the first issue, recently, the Planck Collaboration has reported an upper limit on the sum of all the neutrino mass eigenvalues of m i < 0.23 eV at 95% C.L. in combination with the WMAP polarization and baryon acoustic oscillation (BAO) measurements [4]. Here, the sum runs over all the neutrino mass eigenstates which are in thermal equilibrium in the early Universe. As far as the second issue is concerned, the hunt for neutrinoless double beta decay process is still on which violates the total lepton number and requires Majorana neutrinos [5][6][7][8]. The third question has been answered positively only in recent years [9][10][11] and now, we have compelling evidence in favor of neutrino flavor oscillation [12,13], suggesting that leptonic flavors are not symmetries of Nature. This entails that neutrinos are massive and they mix with each other, providing an evidence for physics beyond the Standard Model (SM) [14,15]. This milestone has been achieved due to several fantastic world-class oscillation experiments involving neutrinos from the Sun [16][17][18][19][20][21][22], the Earth's atmosphere [23][24][25], nuclear reactors [26][27][28][29][30][31][32][33][34][35][36], and accelerators [37][38][39][40][41][42][43][44][45][46] which have enabled us to obtain a clear understanding of the lepton mixing pattern in three-flavor scenario. Using the standard Particle Data Group convention [14], we parametrize the vacuum Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [12,47,48] in terms of the three mixing angles: θ 12 , θ 23 , θ 13 , and one Dirac-type CP phase δ CP (ignoring Majorana phases). In a three-flavor framework, the transition probability also depends on two independent mass-squared differences: ∆m 2 21 ≡ m 2 2 − m 2 1 in the solar sector and ∆m 2 32 ≡ m 2 3 − m 2 2 in the atmospheric sector where m 3 corresponds to the neutrino mass eigenstate with the smallest electron component. The smallest lepton mixing angle θ 13 connects these two sectors and determines the impact of the sub-leading three-flavor effects [49][50][51]. All the neutrino oscillation data available till date have been explained quite successfully in terms of these mass-mixing parameters [52][53][54], excluding few anomalous results obtained at very-short-baseline experiments [55].
Neutrinos acquire additional phases while travelling in matter, the so-called 'MSW effect' [56][57][58][59] which determined the ordering of the 1-2 mass splitting using solar neutrinos. In light of large θ 13 , matter effects are also going to play an important role in presently running and future long-baseline [49,51,60] neutrino oscillation experiments to settle the remaining unsolved issues, namely, the neutrino mass hierarchy 1 , possibility of leptonic CP-violation if δ CP differs from both 0 • and 180 • , and the octant degeneracy of θ 23 [61] if θ 23 turns out to be non-maximal. The combined data from the current off-axis ν e appearance experiments, T2K [62,63] and NOνA [64][65][66], hold promise of providing a first hint for these missing links for only favorable ranges of oscillation parameters [67][68][69][70][71][72][73]. Hence, future facilities consisting of intense, high power, on-axis wide-band beams and large smart detectors have been proposed to cover the entire parameter space at unprecedented confidence level [60]. Future superbeam long-baseline facilities with liquid argon detectors, Deep Underground Neutrino Experiment (DUNE) [74][75][76][77][78] in the United States with a baseline of 1300 km from Fermilab to Homestake mine in South Dakota and Long-Baseline Neutrino Oscillation Experiment (LBNO) in Europe involving a path length of 2290 km between CERN and Pyhäsalmi mine in Finland [79][80][81][82] are the two major candidates in this roadmap which are capable enough to claim the discovery for the above mentioned issues [50,51]. For both the DUNE and LBNO baselines, the matter effects are quite significant which break the eight-fold degeneracies [83,84] among the various oscillation parameters and improve the physics reach by considerable amount.
Apart from the SM W -exchange interaction in matter, there may well be flavordependent, vector-like, leptonic long-range force (LRF), like those mediated by the L e −L µ,τ gauge boson which is very light and neutral, leading to new non-universal flavor-diagonal neutral current (FDNC) interactions of the neutrinos which can give rise to non-trivial three-neutrino mixing effects in terrestrial experiments, and could affect the neutrino propagation through matter [85][86][87][88][89]. Can we constrain/discover these long-range FDNC interactions in upcoming long-baseline neutrino experiments? If this LRF exists in Nature, can it become fatal in our attempts to resolve the remaining unknowns in neutrino oscillation? In this paper, we attempt to address these pressing issues in the context of future high-precision superbeam facilities, DUNE and LBNO.
This paper is organized as follows. We start section 2 with a discussion on flavordependent LRF and how it arises from abelian gauged L e − L µ,τ symmetry. Then, we discuss the strength of the long-range potential V eµ/eτ at the Earth surface generated by the electrons inside the Sun. This is followed by a brief discussion on the current constraints that we have on the effective gauge couplings α eµ, eτ of the L e −L µ,τ symmetry from various neutrino oscillation experiments. In section 3, we study in detail the three-flavor oscillation picture in presence of long-range potential. We derive the compact analytical expressions for the effective oscillation parameters in case of the L e − L µ symmetry, and also present a simple expression for the resonance energy, where 1-3 mixing angle in matter becomes 45 • in the presence of long-range potential. Next, we demonstrate the accuracy of our approximate analytical probability expressions by comparing it with the exact numerical results. We also discuss some salient features of the appearance and disappearance probabilities (calculated numerically) for the Fermilab-Homestake and CERN-Pyhäsalmi baselines in the presence of long-range potential towards the end of this section. We give the similar plots for the anti-neutrino case in appendix A. At the beginning of section 4, we give a brief description of the main experimental features of the DUNE and LBNO set-ups. Then, we study the impact of the long-range potential due to L e − L µ symmetry on the expected event spectra and total event rates for the DUNE and LBNO experiments. In section 5, we describe our simulation strategy. Next, we derive the expected constraints on α eµ from DUNE and LBNO in section 6.1, and estimate the discovery reach for α eµ in section 6.2. In section 6.3, we study how the long-range potential affects the CP-violation search of these future facilities. Section 6.4 is devoted to assess the impact of LRF on mass hierarchy measurements. Finally, we summarize and draw our conclusions in section 7.

Flavor-Dependent Long-Range Forces from Gauged U (1) Symmetries
The SM gauge group SU (3) C × SU (2) L × U (1) Y can be extended with minimal matter content by introducing anomaly free U (1) symmetries under which the SM remains invariant and renormalizable [90]. There are three lepton flavor combinations: L e − L µ , L e − L τ , and L µ − L τ which can be gauged in an anomaly free way with the particle content of the SM [91][92][93]. These U (1) gauge symmetries have to be broken in Nature in order to allow the different neutrino species to mix among each other giving rise to neutrino oscillation as demanded by the recent data [14,[52][53][54]. This can be understood from the following example. If we assume that neutrino masses are generated by effective five-dimensional operator following say, seesaw mechanism, then this operator would remain invariant under these U (1) symmetries if they are exact. This will ultimately give us an effective neutrino mass matrix containing a Dirac and a Majorana neutrino which remain unmixed and there will be no neutrino oscillation. For more discussions on these issues, see references [85,88]. Now, there are two possibilities for the extra gauge boson associated with this abelian symmetry 2 : either it can be very heavy or it is very light but in both the cases, it couples to matter very feebly to escape direct detection. If the electrically neutral gauge boson (Z ) corresponding to a flavor-diagonal generator of this new gauge group is relatively heavy 3 , then the phenomenological consequences can be quite interesting [90][91][92][93]100]. On the other hand, if the mass of the gauge boson is very light, then it can introduce LRF with terrestrial range (greater than or equal to the Sun-Earth distance) without call upon extremely low mass scales [85,86]. This LRF is composition dependent (depends on the leptonic content and the mass of an object) and violate the universality of free fall which could be tested in the classic lunar ranging [101,102] and Eötvos type gravity experiments [103,104]. In fact, this idea was given long back by Lee and Yang [105] and later, using this idea, Okun [106,107] gave a 2σ bound of α < 3.4 × 10 −49 (α denotes the strength of the long-range potential) for a range of the Sun-Earth distance or more. Now, the light gauge boson Z should have a mass m Z ∼ g v where g is the gauge coupling of the new U (1) symmetry and it should be 0.6 × 10 −24 since g ≈ √ α and v is the vacuum expectation value of the symmetry breaking scale. If the range of the LRF is comparable to the Sun-Earth distance (≈ 1.5 × 10 13 cm) then m Z ∼ 1/(1.5 × 10 13 cm) = 1.3 × 10 −18 eV which means v 2 MeV.

Abelian Gauged
If the extra U (1) corresponds to L e −L µ or L e −L τ flavor combination 4 , the coupling of the solar electron to the L e −L µ,τ gauge boson generates a flavor-dependent long-range potential for neutrinos [108][109][110] and can have significant impact on neutrino oscillations [85][86][87][88][89] inspite of such stringent constraints on α. This is caused due to the fact that the (L e −L µ,τ )charge of the electron neutrino is opposite to that of muon or tau flavor, leading to new non-universal FDNC interactions of the neutrinos in matter on top of the SM W -exchange interactions between the matter electrons and the propagating electron neutrinos. These new flavor-dependent FDNC interactions alter the 'running' of the oscillation parameters in matter by considerable amount [111]. Another important point is that the large number of electrons inside the Sun and the long-range nature of interaction balance the smallness of coupling and generate noticeable potential. As an example, the electrons inside the Sun 2 Models with these symmetries necessarily possess a Higgs sector which also discriminates among different lepton families [94], but we will only focus here on the effect of the extra gauge boson. 3 For an example, the Z in gauged Le − Lµ,τ can be produced in e + e − collisions and subsequently decay into e + e − or µ + µ − or τ + τ − pairs and can be constrained using the data from LEP/LEP2 [95][96][97][98][99]. 4 Due to the absence of muons or tau leptons inside the Sun or Earth, we do not consider gauged Lµ − Lτ symmetry here.
give rise to a potential V eµ/eτ at the Earth surface which is given by [85,88] where α eµ/eτ = g 2 eµ/eτ 4π , and g eµ/eτ is the gauge coupling of the L e − L µ,τ symmetry. In Eq. (2.1), N e (≈ 10 57 ) is the total number of electrons inside the Sun [112] and R SE is the Sun-Earth distance ≈ 1.5 × 10 13 cm = 7.6 × 10 26 GeV −1 . Here α eµ/eτ can be identified as the 'fine-structure constant' of the U (1) Le−Lµ,τ symmetry and its value is positive 5 . The corresponding potential due to the electrons inside the Earth with a long-range force having the Earth-radius range (R E ∼ 6400 km) is roughly one order of magnitude smaller compared to the solar long-range potential and can be safely neglected 6 [85,88]. The long-range potential V eµ/eτ in Eq. (2.1) appears with a negative sign for anti-neutrinos and can affect the neutrino and anti-neutrino oscillation probabilities in different fashion. This feature can invoke fake CP-asymmetry like the SM matter effect and can influence the CP-violation search in long-baseline experiments. This is one of the important findings of our paper and we will discuss this issue in detail in the later section. Now, it would be quite interesting to compare the strength of the potential given in Eq. (2.1) with the quantity ∆m 2 /2E which governs the neutrino oscillation probability. For long-baseline neutrinos, ∆m 2 /2E ∼ 10 −12 eV (assuming ∆m 2 ∼ 2×10 −3 eV 2 and E ∼ 1 GeV) which is comparable to V eµ/eτ even for α eµ/eτ ∼ 10 −51 and can affect the long-baseline experiments significantly which we are going to explore in this paper in the context of upcoming facilities DUNE and LBNO.

Existing Phenomenological Constraints on L e − L µ,τ Parameters
There are phenomenological bounds on the effective gauge coupling α eµ/eτ of the L e − L µ,τ abelian symmetry 7 using data from various neutrino oscillation experiments. It was shown in [85] that L e − L µ,τ potential at the Earth due to the huge number of electrons inside the Sun suppresses the atmospheric neutrino ν µ → ν τ oscillations which enabled them to place tight constraints on α eµ/eτ using the oscillation of multi-GeV neutrinos observed at the Super-Kamiokande (SK) experiment. They obtained an upper bound of α eµ < 5.5 × 10 −52 and α eτ < 6.4 × 10 −52 at 90% C.L. [85]. In [88], the authors performed a global fit to the solar neutrino and KamLAND data including the flavor-dependent LRF. They quoted an upper bound of α eµ < 3.4 × 10 −53 and α eτ < 2.5 × 10 −53 at 3σ C.L. assuming θ 13 = 0 • [88]. A similar analysis was performed in [87] to place the constraints on LRF mediated by vector and non-vector (scalar or tensor) neutral bosons where the authors assumed one mass scale dominance. The proposed 50 kt magnetized iron calorimeter (ICAL) detector at the 5 In our work, we consider the case of a light vector boson exchange which makes sure that α eµ/eτ is positive. It means that for an example, the force between an isolated electron and νe is repulsive. 6 The possibility of the local screening of the leptonic force (generated due to the solar electrons) by the cosmic anti-neutrinos is also negligible over the Sun-Earth distance [85]. 7 Flavor-dependent long-range leptonic forces can also be generated via the unavoidable mixing of light Z boson of the Lµ − Lτ symmetry with the Z boson of the SM. This issue was discussed in the context of the MINOS long-baseline experiment in [113,114].
India-based Neutrino Observatory (INO) can also probe the existence of LRF by observing the atmospheric neutrinos and anti-neutrinos separately over a wide range of energies and baselines [89]. With an exposure of one Mton.yr and using the muon momentum as observable, ICAL would be able to constrain α eµ/eτ 1.65 × 10 −53 at 3σ C.L. [89].
3 Three-Flavor Oscillation Picture in Presence of Long-Range Potential In this section, we focus our attention to study the impact of flavor-dependent long-range leptonic potential (due to the electrons inside the Sun) at the Earth surface when neutrinos travel through the Earth matter. In a three-flavor framework, the long-range potential of Eq. (2.1) due to L e − L µ symmetry modifies the effective Hamiltonian for neutrino propagation in Earth matter in the flavor basis to where U is the vacuum PMNS matrix [12,47,48] which can be parametrized as In the above equation, E is the neutrino energy and V CC is the Earth matter potential which appears in the form where G F is the Fermi coupling constant, N e is the electron number density inside the Earth, ρ is the matter density, and Y e = Ne Np+Nn is the relative electron number density. N p , N n are the proton and neutron densities in Earth matter respectively. In an electrically neutral, isoscalar medium, we have N e = N p = N n and Y e comes out to be 0.5. In Eq. (3.1), V eµ is the long-range potential due to L e − L µ symmetry and its strength is given by Eq. (2.1). In case of L e − L τ symmetry, the contribution due to long-range potential in Eq. (3.1) takes the form Diag (V eτ , 0, −V eτ ). Note that the strength of the long-range potential V eµ/eτ does not depend on the Earth matter density 8 and takes the same value for any baseline on the Earth. In case of anti-neutrino propagation, we have to reverse the sign of V CC , V eµ (or V eτ ), and the CP phase δ CP .
To judge the impact of the long-range potential V eµ in long-baseline experiments with multi-GeV neutrinos, we need to compare its strength with the two other main terms in Eq. (3.1) which are ∆m 2 31 2E and the Earth matter potential V CC . In Table 1, we compare the strengths of these three quantities which control the 'running' of the oscillation parameters in matter. The first oscillation maximum for the DUNE (LBNO) experiment occurs at 2.56 GeV (4.54 GeV) assuming ∆m 2 31 = 2.44 × 10 −3 eV 2 (see Table 2). In the fourth

2.40
Not marginalized Table 2. The second column shows the current best-fit values and 1σ uncertainties on the threeflavor oscillation parameters assuming normal hierarchy in the fit [52]. The third column shows the true values of the oscillation parameters used to simulate the 'observed' data set. The fourth column depicts the range over which sin 2 θ 23 and δ CP are varied while minimizing the χ 2 to obtain the final results. In our calculations, ∆m 2 µµ serves as an input parameter and then we estimate the value of ∆m 2 31 using the relation given in Eq. (3.4) (see text for details). In the third column, we take δ CP = 0 • to calculate the value of ∆m 2 31 from ∆m 2 µµ .
column of Table 1, the values of V CC have been estimated using Eq. (3.3) where we take the line-averaged constant Earth matter densities 9 for both the baselines: ρ = 2.87 g/cm 3 for the DUNE baseline and ρ = 3.54 g/cm 3 for the LBNO baseline. Table 1 shows that around first oscillation maximum, the strengths of the terms comparable for both the set-ups under consideration even if α eµ is as small as 10 −52 . It means that these three terms can interfere with each other having significant impact on the oscillation probability which we are going to study in this section with the help of analytical expressions. Before we start deriving our approximate analytical expressions for the effective mass-squared differences and mixing angles in matter in the presence of longrange potential, let us take a look at the current status of the oscillation parameters. The second column of Table 2 shows the present best-fit values and 1σ errors on the three-flavor oscillation parameters assuming normal hierarchy in the fit [52]. We use the benchmark values of the various oscillation parameters as given in the third column of Table 2 to draw the oscillation probability plots in this section and to generate the 'observed' data set while estimating the physics reach of the experimental set-ups. The ranges over which sin 2 θ 23 and δ CP are marginalized while minimizing the χ 2 are given in the fourth column which we will discuss in detail while describing the simulation method in the later section. In Table 2, ∆m 2 µµ is the effective mass-squared difference measured by the accelerator experiments using ν µ → ν µ disappearance channel [41,116] and it is a linear combination of ∆m 2 31 and ∆m 2 21 . The value of ∆m 2 31 is estimated from ∆m 2 µµ using the relation [117,118] ∆m 2 31 = ∆m 2 µµ + ∆m 2 21 (cos 2 θ 12 − cos δ CP sin θ 13 sin 2θ 12 tan θ 23 ) . The value of ∆m 2 31 is calculated separately for NH and IH using the above equation assuming ∆m 2 µµ = ± 2.4 × 10 −3 eV 2 where positive (negative) sign is for NH (IH). Note that through out this paper, whenever we vary δ CP or θ 23 or both, we calculate a new value for ∆m 2 31 using Eq. (3.4).

Analytical Expressions for the Effective Oscillation Parameters
In a CP-conserving scenario (δ CP = 0 • ), the effective Hamiltonian in the flavor basis given in Eq. (3.1) takes the form where H 0 = Diag (0, ∆ 21 , ∆ 31 ) with ∆ 21 ≡ ∆m 2 21 /2E and ∆ 31 ≡ ∆m 2 31 /2E. In the above equation, V = Diag (V CC + V eµ , −V eµ , 0) for L e − L µ symmetry. We can rewrite H f in Eq. (3.5) as where a 11 = A + W + sin 2 θ 13 + α cos 2 θ 13 sin 2 θ 12 , (3.7) a 12 = 1 √ 2 cos θ 13 (α cos θ 12 sin θ 12 + sin θ 13 − α sin 2 θ 12 sin θ 13 ) , (3.8) cos θ 13 (−α cos θ 12 sin θ 12 + sin θ 13 − α sin 2 θ 12 sin θ 13 ) , (3.9) a 22 = 1 2 α cos 2 θ 12 + cos 2 θ 13 − 2α cos θ 12 sin θ 12 sin θ 13 + α sin 2 θ 12 sin 2 θ 13 − 2W , (3.10) a 23 = 1 2 cos 2 θ 13 − α cos 2 θ 12 + α sin 2 θ 12 sin 2 θ 13 , (3.11) a 33 = 1 2 cos 2 θ 13 + α cos 2 θ 12 + α sin θ 13 (sin 2θ 12 + sin 2 θ 12 sin θ 13 ) . (3.12) In the above equations, we introduce the terms A, W , and α which are defined as and we assume that the vacuum value of θ 23 is 45 • . Note that we have kept the terms of all orders in sin θ 13 and α which are quite important in light of the large value of 1-3 mixing angle as was measured recently by the modern reactor experiments. Now, we need to diagonalize the effective Hamiltonian H f in Eq. (3.6) to find the effective mass-squared differences and mixing angles in the presence of the Earth matter potential (V CC ) and the long-range potential (V eµ ). We can almost diagonalize H f with the help of a unitary matrix where off-diagonal terms are quite small and can be safely neglected. The lower right 2 × 2 block in Eq. (3.6) gives us the angle θ m 23 which has the form tan 2θ m 23 = cos 2 θ 13 − α cos 2 θ 12 + α sin 2 θ 12 sin 2 θ 13 W + α sin 2θ 12 sin θ 13 . (3.16) The mixing angles θ m 13 and θ m 12 can be obtained by subsequent diagonalizations of the (1,3) and (1,2) blocks respectively and we get the following expressions (3.17) and where λ 3 = 1 2 cos 2 θ 13 + α cos 2 θ 12 + α sin 2 θ 12 sin 2 θ 13 − W + W + α sin 2θ 12 sin θ 13 cos 2θ m and The eigenvalues m 2 i,m /2E (i = 1, 2, 3) are given by the expressions and (3.24)   Table 2 and we consider δ CP = 0 • .  Table 2. We do the same for the effective mass-squared differences (see Eqs. (3.22), (3.23), and (3.24)) in Fig. 2. The extreme right panel of Fig. 1 shows that θ m 12 approaches to 90 • very rapidly with increasing E in the presence of V CC and this behavior does not change much when we introduce V eµ . This is not the case for θ m 23 and θ m 13 . The long-range potential V eµ affects the 'running' of θ m 23 (see extreme left panel of Fig. 1) significantly. As we go to higher energies, θ m 23 deviates from 45 • and its value decreases very sharply depending on the strength of α eµ . In case of θ m 13 , the effect of V eµ is quite opposite as compared to θ m 23 as can be seen from the middle panel of Fig. 1. Assuming NH, as we increase E, θ m 13 quickly reaches to 45 • (resonance point) in the presence of V eµ and then finally approaches toward 90 • as we further increase E. These two opposite behaviors of θ m 23 and θ m 13 alter the amplitudes and the locations of oscillation maxima in the transition probability substantially for non-zero α eµ which we discuss in the later part of this section. For α eµ = 10 −52 (10 −51 ), the resonance occurs around 4 GeV (0.6 GeV) for 1300 km baseline. We can obtain an analytical expression for the resonance energy demanding θ m 13 = 45 • in Eq. (3.17). In one mass scale dominance approximation where ∆m 2 21 can be neglected i .e. assuming α = 0, the condition for the resonance energy (E res ) takes the form: Now, putting α = 0 in Eqs. (3.19) and (3.16), we get a simplified expression for λ 3 which has the following form since at resonance energies, the term W 2 is small compared to cos 4 θ 13 and can be safely neglected. Now, comparing Eq. (3.25) and Eq. (3.26), we get a simple and compact expression for the resonance energy In the absence of long-range potential V eµ , Eq. (3.27) gives us the standard expression for the resonance energy in the SM framework. Eq. (3.27) suggests that for a given baseline, in the presence of V eµ , the resonance occurs at lower energy as compared to the SM case and this is exactly what we observe in the middle panel of Fig. 1. The right panel of Fig. 2 demonstrates that the effective solar mass-squared difference ∆m 2 21,m increases with energy and can be comparable to the vacuum value of ∆m 2 31 at higher energies in the SM case. In the presence of V eµ , ∆m 2 21,m increases with energy even faster and can become quite large at higher energies depending on the strength of α eµ . In the SM framework, the effective atmospheric mass-squared difference ∆m 2 31,m does not run with energy for this choice of baseline (see left panel of Fig. 2). But, in the presence of V eµ , the value of ∆m 2 31,m enhances a lot depending on the strength of α eµ as energy is increased. Next, to check the accuracy of our approximate analytical results, we compare the oscillation probabilities calculated with our approximate effective running mixing angles and mass-squared differences with those calculated numerically for the same baseline and line-averaged constant matter density along it.

Demonstration of the Accuracy of the Approximation
The neutrino oscillation probabilities in the presence of V CC and V eµ can be obtained by replacing the vacuum expressions of the elements of the mixing matrix U and the masssquare differences ∆m 2 ij with their effective 'running' values Incorporating the modifications due to V CC and V eµ , the new transition probability in a CP-conserving scenario can be written as Using Eq. (3.29), we obtain the following expressions for the appearance and disappearance channels [119] P (ν µ → ν e ) = 4Ũ 2  Figure 3. ν µ → ν e transition probability as a function of neutrino energy E in GeV for 1300 km (2290 km) baseline in left (right) panels. The upper panels are for the SM case without longrange potential. The lower panels correspond to α eµ = 10 −52 . In all the panels, we compare our analytical expressions (dashed curves) to the exact numerical results (solid curves) for NH and IH.
The vacuum values of the oscillation parameters are taken from the third column of Table 2 and we take δ  . ν µ → ν µ transition probability as a function of neutrino energy E in GeV for 1300 km (2290 km) baseline in left (right) panels. The upper panels are for the SM case without longrange potential. The lower panels correspond to α eµ = 10 −52 . In all the panels, we compare our analytical expressions (dashed curves) to the exact numerical results (solid curves) for NH and IH.
The vacuum values of the oscillation parameters are taken from the third column of Table 2 and we take δ CP = 0 • .
In Fig. 3, we present our approximate ν µ → ν e oscillation probabilities (dashed curves) as a function of the neutrino energy against the exact numerical results (solid curves) considering L = 1300 km (left panels) and 2290 km (right panels). We give the plots for both NH and IH in all the panels considering line-averaged constant Earth matter densities for both the baselines. The upper panels are drawn for the SM case (α eµ = 0) where our approximate results match exactly with the numerically obtained probabilities.
In the lower panels, we give the probabilities considering α eµ = 10 −52 and find that our approximate expressions work quite well in the presence of long-range potential and can predict almost accurate L/E patterns of the oscillation probability. In Fig. 4, we study the same for ν µ → ν µ oscillation channel and find that our approximate expressions match quite nicely with the numerical results. Here, we present our analytical results for NH. We can obtain the same for IH by changing ∆m 2 31 → −∆m 2 31 . Following the same procedure and reversing the sign of V in Eq. (3.5), we can derive the analytical expressions for antineutrino as well. Note that in this paper, we limit our investigation to L e − L µ symmetry, though similar procedure can be adopted for L e − L τ symmetry.

Discussion at the Probability Level -Neutrino Case
In this section, we discuss in detail how the long-range potential affects the full three-flavor neutrino oscillation probabilities in matter considering non-zero values of δ CP . In Fig. 5, we show the exact numerical transition probability P µe as a function of the neutrino energy using the line-averaged constant Earth matter densities for 1300 km (left panels) and 2290 km (right panels) baselines. We vary δ CP within the range −180 • to 180 • and the resultant probability is shown as a band, with the thickness of the band reflecting the effect of δ CP on P µe . Inside each band, the probability for δ CP = 0 • case is shown explicitly by the black dashed line. The left panels (right panels) are for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential.
In the upper (lower) panels, we take α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with longrange potential. We study the same for the disappearance (ν µ → ν µ ) channel in Fig. 6. We give the similar plots for the anti-neutrino case in appendix A.
To explain the behavior of the oscillation probabilities in Fig. 5 and Fig. 6 at the qualitative level, we can simplify the analytical expressions given in Eqs. (3.30) and (3.31) in the following fashion. The extreme right panel of Fig. 1 suggests that sin θ m 12 → 1 and cos θ m 12 → 0 very quickly as we increase E in the SM case or with non-zero α eµ . So, we set sin θ m 12 ≈ 1 and cos θ m 12 ≈ 0 in Eqs. (3.30) and (3.31) and obtain the following simple expressions: and In Fig. 5, we can see that for α eµ = 10 −52 case (upper panels), the locations of the first oscillation maxima have been shifted toward lower energies for both the baselines and also the amplitudes of the first oscillation maxima have been enhanced when we assume NH. This can be understood from the 'running' of θ m 23 , θ m 13 (extreme left and middle panels of The transition probability P µe as a function of neutrino energy. The band reflects the effect of unknown δ CP . Inside each band, the probability for δ CP = 0 • case is shown by the black dashed line. The left panels (right panels) are for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential. In the upper (lower) panels, we take α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with long-range potential.
increases and θ m 23 decreases and there is a trade-off between the terms sin 2 θ m 23 and sin 2 2θ m 13 in Eq. (3.32). Also, the value of ∆m 2 32,m (∆m 2 31,m -∆m 2 21,m ) decreases with energy as ∆m 2 21,m increases by substantial amount compared to ∆m 2 31,m , which shifts the location of the first oscillation maxima toward lower energies. For IH, the value of θ m 13 decreases fast with non-zero α eµ compared to the SM case, causing a depletion in the probabilities over a wide range of energies. In case of α eµ = 10 −51 (lower panels), there is a huge suppression in the probabilities at both the baselines over a wide range of energies above 1 GeV assuming NH. The main reason behind this large damping in the probabilities is that θ m 13 approaches very quickly to 90 • around 1 GeV or so for α eµ = 10 −51 (see middle panel of Fig. 1)  Figure 6. The transition probability P µµ as a function of neutrino energy. The band reflects the effect of unknown δ CP . Inside each band, the probability for δ CP = 0 • case is shown by the black dashed line. The left panels (right panels) are for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential. In the upper (lower) panels, we take α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with long-range potential. and therefore, sin 2 2θ m 13 → 0, vanishing the probability amplitude for ν µ → ν e oscillation channel. Below 1 GeV, θ m 13 runs toward 45 • and therefore, sin 2 2θ m 13 → 1, causing the enhancement in the probabilities. When we take IH, θ m 13 quickly advances to zero, causing a huge damping in the probabilities at all the energies. These 'running' behaviors of θ m 23 , θ m 13 , and the mass-squared differences in the presence of long-range potential as discussed above also affect ν µ → ν µ oscillation channel (see Fig. 6) which can be explained with the help of Eq. (3.33). Next, we discuss how the long-range potential due to L e − L µ symmetry modifies the expected event spectra and total event rates of the DUNE and LBNO experiments.

Impact of Long-Range Potential at the Event Level
We start this section with a brief description of the main experimental features of the DUNE and LBNO set-ups that we use in our simulation.

Key Features of DUNE and LBNO Set-ups
The proposed DUNE experiment [74][75][76][77][78]120] in the United States with a baseline of 1300 km from Fermilab to Homestake mine in South Dakota is planning to build a massive 35 kt liquid argon time projection chamber (LArTPC) as the far detector choice. This LArTPC will have excellent kinematic reconstruction capability for all the observed particles, rejecting almost all of the large neutral current background. We use the detector properties which are given in Table 1 of Ref. [79]. As far as the neutrinos are concerned, this facility will have a new, high intensity, on-axis neutrino beam, which in its initial phase, will operate at a proton beam power of 708 kW, with proton energy of 120 GeV, delivering 6 × 10 20 protons on target in 230 days per calendar year. In this work, we have used the latest fluxes being considered by the collaboration [121]. We have assumed five years of neutrino run and five years of anti-neutrino run to estimate the physics capabilities of this set-up.
In Europe, the proposed LBNO experiment [79][80][81][82] offers an interesting possibility to address the fundamental unsolved issues in neutrino oscillation physics using a baseline of 2290 km between CERN and Pyhäsalmi mine in Finland which enables us to cover a wide range of L/E choices, mandatory to resolve parameter degeneracies. The Pyhäsalmi mine will house a giant 70 kt LArTPC as a far detector which will observe the neutrinos produced in a conventional wide-band beam facility at CERN. The fluxes that we use in our simulation have been computed assuming an exposure of 1.5 × 10 20 protons on target in 200 days per calendar year from the SPS accelerator at 400 GeV with a beam power of 750 kW [122]. For LBNO also, we assume five years of neutrino run and five years of anti-neutrino run. We consider the same detector properties as that of DUNE.

Event Spectrum and Rates
In this section, we present the expected event spectra and total event rates for both the setups under consideration in the presence of long-range potential. We calculate the number of expected electron events 10 in the i-th energy bin in the detector using the following expression where φ(E) is the neutrino flux, T is the total running time, n n is the number of target nucleons in the detector, is the detector efficiency, and R(E, E A ) is the Gaußian energy resolution function of the detector. σ νe is the neutrino interaction cross-section which has been taken from Refs. [123,124], where the authors estimated the cross-section for water and isoscalar targets. In order to have LAr cross-sections, we have scaled the inclusive charged current cross-sections of water by a factor of 1.06 for neutrino and 0.94 for antineutrino [125,126]. The quantities E and E A are the true and reconstructed (anti-)neutrino energies respectively, and L is the path length.  In both the panels, the solid grey (blue) vertical lines display the locations of the first (second) oscillation maxima. We assume δ CP = 0 • and NH. For other oscillation parameters, the values are taken from the third column of Table 2.
In our study, we consider the ν e andν e appearance channels, where the backgrounds mainly stem from the the intrinsic ν e /ν e contamination of the beam, the number of muon events which will be misidentified as electron events, and the neutral current events. In Fig. 7, we show the expected signal and background event spectra as a function of reconstructed neutrino energy including the efficiency and background rejection capabilities. The left panel shows the results for the DUNE set-up with 35 kt far detector mass. The right panel displays the same for the LBNO set-up with 70 kt far detector. In both the panels, the thick lines correspond to the SM case, whereas the thin lines are drawn assuming α eµ = 10 −52 . In both the set-ups, we can clearly see a systematic downward bias in the reconstructed energy for the neutral current background events due to the final state neutrino included via the migration matrices. The solid grey (blue) vertical lines display the locations of the first (second) oscillation maxima. The red solid histogram shows the signal event spectrum. Note that in the presence of long-range potential, both the signal and background (intrinsic ν e contamination and misidentified muons) event spectra get modified substantially. For both the baselines, we have considerable number of signal events around the second oscillation maximum. But, these event samples are highly contaminated with the neutral current and other backgrounds at lower energies, limiting their impact.   Table 4. Comparison of the total signal and background event rates in the ν µ /ν µ disappearance channel for DUNE (35 kt) and LBNO (70 kt) set-ups. For the cases denoted by 'SM+LRF', we take α eµ = 10 −52 . The results are shown for both NH and IH assuming δ CP = 0 • . For both the set-ups, we consider five years of neutrino run and five years of anti-neutrino run.
In Table 3, we show a comparison between the total signal and background event rates in the ν e /ν e appearance channel for DUNE (35 kt) and LBNO (70 kt) set-ups. For both the set-ups, we assume five years of neutrino run and five years of anti-neutrino run. For the cases denoted by 'SM+LRF', we take α eµ = 10 −52 . The results are shown for both NH and IH assuming δ CP = 0 • . The Earth matter effects play an important role for both the baselines which is evident from the fact that in the neutrino channel, the number of expected events is quite large compared to the IH case and in the anti-neutrino channel, the situation is totally opposite where we have larger event rates for IH than for NH. The relative difference between the number of events for NH and IH is larger for the CERN-Pyhäsalmi baseline than the FNAL-Homestake baseline, since the impact of matter effects is more significant at the 2290 km baseline compared to the 1300 km baseline. In the presence of long-range potential with a benchmark choice of α eµ = 10 −52 , qualitatively, the trend remains the same as mentioned above. Table 3 clearly shows that in all the cases, the most dominant contribution to the background comes from the intrinsic ν e /ν e beam contamination. Note that though the total signal event rate for the DUNE set-up in the neutrino mode with NH, does not change much due to long-range potential with α eµ = 10 −52 , but, the shape of the signal event spectrum gets affected by considerable amount as can be seen from Fig. 7, which enables us to place tight constraints on α eµ as we discuss in the results section. In our simulation, we also include the information coming from the ν µ /ν µ disappearance channels. For these type of channels, neutral current events are the main source of background. Table 4 shows the total signal and background event rates in the ν µ /ν µ disappearance channels for both the set-ups, considering five years of neutrino run and five years of anti-neutrino run. For the cases marked by 'SM+LRF', we take α eµ = 10 −52 and we present results for both NH and IH assuming δ CP = 0 • . Interestingly, the ν µ /ν µ disappearance channels are also quite sensitive to the long-range potential and in all the cases, we see a significant change in the total signal event rates with α eµ = 10 −52 as compared to the SM case. Also, the rates are different for NH and IH in the presence of long-range potential. The ν µ /ν µ disappearance channels also play an important role to constrain the atmospheric oscillation parameters in the fit.

Bi-events Plot
In this section, we make an attempt to unravel the impact of long-range potential with the help of bi-events plot. In Fig. 8, we have plotted ν e vs.ν e appearance events, for DUNE (left panel) and LBNO (right panel), considering both NH and IH and with and without long-range potential. Since δ CP is not known, events are generated for the full range [−180 • , 180 • ], leading to the ellipses. For the cases labelled by 'SM+LRF', we take α eµ = 10 −52 . We generate these plots with sin 2 θ 23 = 0.5 as mentioned in Table 2. The ellipses in Fig. 8 suggest that both the set-ups can discriminate between NH and IH at high confidence level, irrespective of the choice of δ CP , and the presence of long-range potential with α eµ = 10 −52 does not spoil this picture. We can see from the left panel that for the DUNE set-up, the anti-neutrino (neutrino) event rates get reduced for NH (IH) with LRF as compared to the SM case. But, for CERN-Pyhäsalmi baseline (see right panel) with more matter effect, both the neutrino and anti-neutrino event rates get diminished for NH and IH in the presence of long-range potential as compared to the SM case. Now, let us make an attempt to understand this behavior. In the presence of long-range potential, the locations of the first oscillation maxima shift towards lower energies (see the upper panels of Fig. 5 and Fig. 16), where both the fluxes and the interaction cross-sections are small. On the contrary, at higher energies, we see a suppression in the probability where we have most of the neutrino fluxes and the cross-sections are also high at these energies. These opposite bahaviors are responsible for the large depletions in the event rates. Fig. 8 also portrays that the asymmetries between the neutrino and anti-neutrino appearance events are largest for the combinations: (NH, δ CP = −90 • ) and (IH, δ CP = 90 • ). One striking feature emerging from both the panels is that all the ellipses get shrunk in the presence of long-range potential, reducing the differences in the number of events due to the CP-conserving and CP-violating phases. It ultimately affects the CP-coverage for the leptonic CP-violation searches, where we study the choices of the CP phase, δ CP which can be distinguished from both 0 • and 180 • at a given confidence level. The right panel shows that this effect is more prominent for the LBNO set-up, severely limiting its discovery reach for CP-violation which we discuss in detail in the results section.

Simulation Method
In this section, we give a brief description of the numerical technique and analysis procedure which we adopt to estimate the physics reach of the experimental set-ups. We have made suitable changes in the GLoBES software [127,128] to obtain our results. The entire numerical analysis is performed using the full three-flavor oscillation probabilities. Unless stated otherwise, we generate our simulated data considering the true values of the oscillation parameters given in the third column of Table 2. These choices of the oscillation parameters are well within their 1σ allowed ranges which are obtained in recent global fit analysis [52]. In the fit, we marginalize over test sin 2 θ 23 and test δ CP in their allowed ranges which are given in the fourth column of Table 2, without assuming any prior on these parameters. We also marginalize over both the hierarchy choices in the fit for all the analyses, except for the mass hierarchy discovery studies where our goal is to exclude the wrong hierarchy in the fit. We keep θ 13 fixed in the fit as the Daya Bay experiment is expected to achieve a relative 1σ precision of ∼ 3% by the end of 2017 [129], and needless to say that the global oscillation data will severely constrain θ 13 beyond the Daya Bay limit before these future experiments will come online. For the atmospheric mass-squared splitting, we take the true value of ∆m 2 µµ = ± 2.4 × 10 −3 eV 2 where positive (negative) sign is for NH (IH), and we do not marginalize over this parameter in the fit since the projected combined data from the currently running T2K and NOνA experiments will be able to improve the precision in |∆m 2 µµ | to sub-percent level for maximal θ 23 [130]. On top of the standard three-flavor oscillation parameters, we have also the LRF parameter α eµ which enters into the oscillation probability. As far as the true and test values of α eµ are concerned, we vary our choices in the range 10 −54 to 10 −51 , where the lower limit 11 corresponds to the cases where the oscillation probabilities almost overlap with the SM cases for the set-ups that we consider in this work. We take the upper limit of α eµ as 10 −51 which covers all the existing bounds on this parameter available from the oscillation experiments as discussed in section 2.2. Based on the techniques discussed in Refs. [131,132], we use the following χ 2 function in our statistical analysis: where n is the total number of bins and Above, N th i ({ω, α eµ }) is the predicted number of signal events in the i-th energy bin for a set of oscillation parameters ω and a particular value of α eµ . N b i ({ω, α eµ }) is the number of background events in the i-th bin where the charged current backgrounds are dependent on ω and α eµ , and the neutral current backgrounds do not depend on the oscillation parameters and α eµ . The quantities π s and π b in Eq. (5.2) are the systematic errors on the signal and background respectively. For both the set-ups, we consider π s = 5% and π b = 5% in the form of normalization error for both the appearance and disappearance channels. The quantities ξ s and ξ b are the "pulls" due to the systematic error on signal and background respectively. We incorporate the data in Eq. (5.1) through the variable i is the number of observed charged current signal events in the i-th energy bin and N b i is the background as mentioned earlier. To estimate the total χ 2 , we add the χ 2 contributions coming from all the relevant channels in a given experiment in the following way  where we assume that all these channels are completely uncorrelated, all the energy bins in a given channel are fully correlated, and the systematic errors on signal and background are fully uncorrelated. Finally, χ 2 total is marginalized in the fit over the oscillation parameters, both the hierarchy choices, the LRF parameter α eµ (as needed), and the systematic parameters as mentioned above to obtain ∆χ 2 min .

Results
In this section, we report our main findings. First, we present the expected constraints on α eµ from the proposed DUNE and LBNO experiments. Next, we quantify the discovery reach for α eµ of these future facilities. Then, we address how the flavor-dependent LRF, mediated by the extremely light L e − L µ gauge boson can affect the CP-violation searches and the mass hierarchy measurements at these upcoming facilities.

Expected Constraints on the Effective Gauge Coupling α eµ
In this section, we estimate the upper bounds on α eµ from the proposed DUNE and LBNO experiments if there is no signal of LRF in the data. This performance indicator corresponds to the new upper limit on α eµ if the experiment does not see a signal of LRF in oscillations. We simulate this situation in our analysis by generating the data at α eµ (true) = 0 and fitting it with some non-zero value of α eµ by means of the χ 2 technique as outlined in section 5. The corresponding ∆χ 2 Bound obtained after marginalizing over oscillation parameters (θ 23 , δ CP , and mass hierarchy) and systematic parameters in the fit, is plotted in Fig. 9 as a function of α eµ (test), which gives a measure of the sensitivity reach of the DUNE or LBNO In all the cases, we assume NH as true hierarchy. Fig. 9 clearly shows that the LBNO set-up with 70 kt detector mass can place better limits on α eµ as compared to the DUNE set-up with 35 kt detector and the limits are not very sensitive to the choice of unknown δ CP (true) for both the set-ups. Table 5 lists the precise upper limits on α eµ which are expected from these future facilities if there is no trace of LRF in the data. We present the bounds at 90% (1.64σ) and 3σ confidence levels 12 for four different choices of true values of δ CP : 0 • , 180 • , 90 • , and −90 • . For each δ CP (true) value, we give the results for both NH and IH as true hierarchy choice. For an example, if δ CP (true) = −90 • and true hierarchy is NH, then the 90% C.L. limit from the DUNE (LBNO) experiment is α eµ < 1.9 × 10 −53 (7.8 × 10 −54 ), suggesting that the constraint from the LBNO experiment is ∼ 2.4 times better than the DUNE set-up 13 . This future limit from the DUNE (LBNO) experiment is ∼ 30 (70) times better than the existing limit 14 from the SK experiment [85] which is also mainly sensitive to the atmospheric mass scale like long-baseline experiments. At 3σ, we see the same relative improvement in the LBNO experiment in constraining α eµ compared to the DUNE set-up. Table 5 also suggests that the limits are not highly dependent on the true choices of mass hierarchy. The variation in the upper limits on α eµ is also not significant as we vary δ CP (true) in the range -180 • to 180 • as can be seen from Fig. 10, where we give the expected bounds at 3σ and 5σ confidence levels from both the set-ups assuming NH as true hierarchy. Next, we discuss the discovery reach for α eµ if we find a positive signal of LRF in the expected event spectra at DUNE and LBNO.

Discovery Reach for α eµ
How good are our chances of observing a positive signal for LRF and hence α eµ in these proposed facilities? We answer this question in terms of the parameter indicator which we call the "discovery reach" of the experiment for α eµ . We define this performance indicator as the expected lower limit on true values of α eµ above which the projected data at DUNE or LBNO would give us a signal for LRF at a certain confidence level. To find these limiting values, we simulate the data for various true values of α eµ and fit it with a predicted event spectrum corresponding to α eµ = 0. We marginalize over θ 23 , δ CP , mass hierarchy, and systematic parameters in the fit to estimate the resultant ∆χ 2 Discovery which is plotted in Fig. 11 for DUNE (35 kt) and LBNO (70 kt) set-ups, considering four different choices of true values of δ CP . In the left panel, we take the CP-conserving choices: δ CP (true) = 0 • (solid lines), 180 • (dashed lines). In the right panel, we consider the maximal CP-violating choices: δ CP (true) = 90 • (solid lines), −90 • (dashed lines). In all the cases, we take NH as true hierarchy. The nature of the curves in Fig. 11 are quite similar to the curves which are shown in Fig. 9, and LBNO with 70 kt detector has better discovery reach for α eµ as   compared to DUNE with 35 kt, like in the case of constraints on α eµ . In Table 6, we give the precise lower limits on true values of α eµ which can be separated from α eµ = 0 in the fit at 90% and 3σ confidence levels. The results are given for both the set-ups and for four different choices of true values of δ CP : 0 • , 180 • , 90 • , and −90 • . For each δ CP (true) value, we show the results for both the choices of true hierarchy: NH and IH. If we compare the entries in Table 6 and Table 5, then we can see that the values of the discovery reach for α eµ are slightly different than the constraints on α eµ at a given confidence level and for the same choices of true oscillation parameters. Also, we can see that the values of discovery reach are marginally dependent on the choices of true δ CP and mass hierarchy as we have seen for the constraints in the previous section.  Table 7. Fraction of δ CP (true) for which a discovery is possible for CP-violation from DUNE (35 kt) and LBNO (70 kt) set-ups at 2σ and 3σ confidence levels. We show the coverage in δ CP (true) for both the choices of true hierarchy: NH and IH. For the 'SM' cases, we consider α eµ = 0 in the data and also in the fit. We also give the results generating the data with α eµ (true) = 10 −52 and in the fit, we marginalize over test values of α eµ in its allowed range. The rest of the simulation details are exactly similar to the 'SM' case (see text for details).

How Robust are CP-violation Searches in Presence of LRF?
This section is devoted to study how the long-range potential due to L e − L µ symmetry affects the CP-violation search which is the prime goal of these future facilities. Can we reject both the CP-conserving values of 0 • , 180 • at a given confidence level? The performance indicator "discovery reach of leptonic CP-violation" addresses this question and obviously, this measurement becomes extremely difficult for the δ CP values which are close to 0 • and 180 • . In Fig. 12, we present the CP-violation discovery reach of DUNE (left panel) and LBNO (right panel) as a function of true value of δ CP assuming NH as true hierarchy. In this plot, we generate our predicted event spectrum (data) considering the true value of δ CP as shown in the x-axis, along with the other true values of the oscillation parameters given in the third column of Table 2. Then, we estimate the various theoretical event spectra assuming the test δ CP to be the CP-conserving values 0 • and 180 • , and by varying simultaneously θ 23 in its 3σ allowed range and both the choices of mass hierarchy. We calculate the ∆χ 2 between each set of predicted and theoretical event spectra using the procedure described in section 5. The smallest of all such ∆χ 2 values: ∆χ 2 CPV is plotted in Fig. 12 as a function of δ CP (true) in the range -180 • to 180 • . In both the panels, the solid red lines depict the 'SM' case where α eµ = 0 in the data and also in the fit. For each δ CP (true), we also give the results generating the data with three different true values of α eµ which are mentioned in the figure legends. In all these three cases, in the fit, we also marginalize over test values of α eµ in the range 10 −54 to 10 −51 along with the other three-flavor oscillation parameters as discussed before. Fig. 12 clearly shows that the CP-violation discovery reach can be altered by substantial amount as compared to the 'SM' case depending on the true choice of α eµ . In case of α eµ (true) = 6 × 10 −53 , we see a large suppression in the CP-violation discovery reach of DUNE (left panel) in the range 45 • ≤ δ CP (true) ≤ 135 • . We have checked that this mainly happens due to the marginalization over θ 23 in the fit where we vary sin 2 θ 23 over a wide range (0.38 to 0.64) without imposing any prior on it. In case α eµ (true) = 10 −52 , the LBNO set-up (right panel) suffers a large depletion in the CP-violation discovery reach which can be easily explained with the help of bi-events plot (Fig. 8) shown in section 4.3. In the right panel of Fig. 8, we have seen a large reduction in the ν andν event rates for CERN-Pyhäsalmi baseline with α eµ = 10 −52 and this is true for both NH and IH. The differences in the number of events for the CP-conserving and CP-violating phases get reduced as the ellipses in Fig. 8 get shrunk in the presence of LRF, severely deteriorating the CP-violation discovery reach of LBNO as can be seen from the right panel of Fig. 12. Table 7 also validates this result, where we compare the precise fraction of δ CP (true) for which a discovery is possible for CP-violation from LBNO (70 kt) and DUNE (35 kt) at 2σ and 3σ confidence levels. For LBNO set-up with true NH and α eµ (true) = 10 −52 , the coverage in δ CP (true) at 3σ C.L. reduces to 30% from 55% as we have in the 'SM' case. In case of true IH, the impact of long-range potential is even more dramatic for these future facilities. At 3σ with α eµ (true) = 10 −52 , their CP-violation reach is quite minimal: only 12% for LBNO and 37% for DUNE while in the 'SM' framework, the coverage is 60% for LBNO and 53% for DUNE. Since, the sign of the long-range potential V eµ is opposite for neutrino and anti-neutrino, it affects the neutrino and anti-neutrino oscillation probabilities in different fashion. This feature introduces fake CP-asymmetry like the SM matter effect and severely limits the CP-violation search in these long-baseline facilities which can be clearly seen from Table 7.
Finally, to see the complete picture, the fraction of δ CP (true) for which a discovery is possible for CP-violation is shown in Fig. 13 as a function of true value of α eµ assuming NH as true hierarchy. In each panel, we compare the performances of DUNE (35 kt) and LBNO (70 kt) which are shown by red and blue lines respectively. We give the results at 2σ (left panel) and 3σ (right panel) confidence levels. In both the panels, the solid horizontal lines depict the 'SM' case where α eµ = 0 in the data and also in the fit. For the 'SM+LRF' case (dashed lines), the data is generated with the true value of α eµ as shown in the x-axis and in the fit, we marginalize over test values of α eµ in its allowed range. The rest of the simulation details are exactly similar to the 'SM' case as discussed before. For the values close to α eµ (true) = 10 −54 , the event spectra in the data is almost similar to the 'SM' case, but since we allow α eµ to vary in the fit in the range 10 −54 to 10 −51 along with the other three-flavor oscillation parameters as discussed before, we see a small suppression in the fraction of δ CP (true). In both the panels, around α eµ (true) = 6 × 10 −53 , the CP-violation discovery reach of DUNE deteriorates substantially, which we also observe in Fig. 12, and the marginalization over θ 23 is mainly responsible for this as we have already discussed. Once α eµ (true) approaches toward 10 −52 , the coverages in δ CP (true) for which CP-violation can be observed, shrink very rapidly for both the set-ups, and ultimately around α eµ (true) = 2 × 10 −52 , the coverages almost become zero. We can understand this feature from our discussions in section 3.1, where we have seen that in the presence of V eµ , as we increase E, θ m 13 quickly approaches toward 45 • (see middle panel of Fig. 1), and the resonance occurs at much lower energies as compared to the SM case. Finally, θ m 13 reaches to 90 • as we further increase E, and the ν µ → ν e oscillation probabilities vanish for most of the energies where we have significant amount of neutrino flux. It causes a huge suppression in the event rates and as a result, the sensitivity goes to zero. Next, we turn our attention to the mass hierarchy discovery potential of DUNE and LBNO. The rest of the simulation details are exactly similar to the 'SM' case (see text for details). Also, note that the ranges in the x-axis and y-axis are different in some of the panels.

Impact of LRF on Mass Hierarchy Measurements
The large Earth matter effects at both the DUNE and LBNO baselines enhance the separation between the oscillation spectra of NH and IH, and hence, we have large differences in the event rates for NH and IH, leading to unprecedented sensitivity toward neutrino mass hierarchy. Now, it would be quite interesting to see how robust are these measurements in the presence of LRF? A 'discovery' of the mass hierarchy is a discrete measurement and is defined as the ability to exclude any degenerate solution for the wrong (fit) hierarchy at a given confidence level. For hierarchy sensitivity, we first assume NH to be the true hierarchy and we choose a true value of δ CP and α eµ . We compute the NH event spectrum for these assumptions and the other true values of the oscillation parameters (see the third column of Table 2) and label it to be data. Then, we estimate the various theoretical event spectra assuming IH and a test value of α eµ as shown in the x-axis of Fig. 14, and by varying simultaneously test θ 23 in its 3σ allowed range and test δ CP in the full allowed range (-180 • to 180 • ). Next, we compute the ∆χ 2 between each set of predicted and theoretical event spectra using the numerical technique described in section 5. The smallest of all  such ∆χ 2 values: ∆χ 2 MH is shown in Fig. 14 as a function of α eµ (test) for given choices of α eµ (true) and δ CP (true). As mentioned above, we always assume NH in the data and IH in the fit while generating the curves in Fig. 14. The upper panels portray the performance of DUNE (35 kt), while the lower panels are for LBNO (70 kt). In each panel, the results are given for four different choices of δ CP (true) and the solid horizontal lines depict the 'SM' case where α eµ is zero in the data and also in the fit. For the 'SM+LRF' case (dashed lines), the data is generated with the true value of α eµ as mentioned in the top part of each panel, and in the fit, we vary the test values of α eµ while marginalizing over θ 23 and δ CP . The rest of the simulation details are exactly similar to the 'SM' case as mentioned above. Though in Fig. 14, we have shown the results for three benchmark values of α eµ (true) = 10 −53 (left panels), 10 −52 (middle panels), and 10 −51 (right panels), but, we have checked that for DUNE, the mass hierarchy sensitivity always stays above the standard expectations irrespective of δ CP (true) provided the true value of α eµ < 5 × 10 −52 . There is a large suppression in the appearance event rates when we generate the data with a true value of α eµ around 5 × 10 −52 , and if we further increase the value of α eµ (true), the statistical strength of the data reduces very rapidly, and the sensitivity goes below the standard expectation. The upper right panel in Fig. 14 clearly shows this behavior. In case of LBNO, the mass hierarchy discovery reach never goes below the 'SM' value irrespective of δ CP (true) if the true choice of α eµ is smaller than 10 −52 . Once we consider the true value of α eµ ≥ 10 −52 , the appearance event rates get reduced in data by considerable amount, causing a significant drop in the sensitivity which can be clearly seen from the lower middle and right panels in Fig. 14. In Fig. 15, we generate the data with IH and fit it with NH. We see almost similar behavior in all the panels of Fig. 15 as we have noticed in Fig. 14.

Summary and Conclusions
Flavor-dependent long-range leptonic forces mediated by the extremely light and neutral bosons associated with gauged L e − L µ or L e − L τ symmetries, constitute a minimal extension of the SM preserving its renormalizability and can lead to interesting phenomenological consequences. For an example, the electrons inside the Sun can generate a flavor-dependent long-range potential V eµ/eτ at the Earth surface which can give rise to non-trivial three neutrino mixing affects in terrestrial experiments, and could influence the neutrino propagation through matter. The sign of this potential is opposite for anti-neutrinos, and affects the neutrino and anti-neutrino oscillation probabilities in different fashion. This feature invokes fake CP-asymmetry like the SM matter effect and can severely affect the leptonic CP-violation searches in long-baseline experiments. In this paper for the first time, we have investigated in detail the possible impacts of these long-range flavor-diagonal neutral current interactions in the oscillations of neutrinos and anti-neutrinos in the context of future high-precision superbeam facilities, DUNE and LBNO. The key point here is that for long-baseline neutrinos, ∆m 2 /2E ∼ 2.5 × 10 −13 eV (assuming ∆m 2 ∼ 2.5 × 10 −3 eV 2 and E ∼ 5 GeV) which is comparable to V eµ even for α eµ ∼ 10 −52 , and can influence the long-baseline experiments significantly. For the Fermilab-Homestake (1300 km) and CERN-Pyhäsalmi (2290 km) baselines, the Earth matter potentials are also around 10 −13 eV (see Table 1), suggesting that V CC can also interfere with V eµ and ∆m 2 31 /2E, having substantial impact on the oscillation probability. We have explored these interesting possibilities in detail in this work. In this paper, we have presented all the results considering the L e − L µ symmetry. Similar analysis can be performed for the L e − L τ symmetry which we will present elsewhere.
We have derived approximate analytical expressions for the effective neutrino oscillation parameters to study how they 'run' as functions of the neutrino energy in the presence of both long-range and Earth matter potentials. We have also obtained a compact and simple expression for the resonance energy, where θ 13 becomes 45 • in the presence of both V CC and V eµ . We have observed that in the presence of V eµ , as we increase the neutrino energy, θ 13 in matter quickly approaches toward 45 • , and the resonance occurs at much lower energies as compared to the SM case. Finally, θ 13 in matter reaches to 90 • as we further increase the energy, causing a large suppression in the appearance probability for most of the energies where we have significant amount of neutrino flux for both the set-ups.
As a result, the event rates get reduced which can be clearly seen from the bi-events plot in Fig. 8.
As the long-range potential due to gauged L e − L µ symmetry can change the standard oscillation picture of these future facilities significantly, we can expect to place strong constraints on α eµ if these experiments do not observe a signal of LRF in oscillations. For an example, if δ CP (true) is −90 • and true hierarchy is NH, then the expected bound from the DUNE (35 kt) set-up at 90% C.L. is α eµ < 1.9 × 10 −53 . The same from the LBNO (70 kt) experiment is α eµ < 7.8 × 10 −54 , suggesting that the constraint from LBNO is 2.4 times better than DUNE. This future limit from the DUNE (LBNO) experiment is almost 30 (70) times better than the existing bound from the SK experiment [85]. We have noticed that these future limits on α eµ from DUNE and LBNO are not very sensitive to the true choice of δ CP and mass hierarchy. We have also estimated the discovery reach for α eµ if we find a positive signal of LRF in the expected event spectra at DUNE and LBNO. We have found that the spectral information on the signal and background events is quite crucial to constrain/discover this new long-range force.
We have also studied in detail the CP-violation discovery reach of DUNE (35 kt) and LBNO (70 kt) in the presence of LRF. We have seen that the CP-violation measurements can be deteriorated by considerable amount as compared to the standard expectation depending on the true value of α eµ . At 3σ with α eµ (true) = 10 −52 and true NH, the coverage in δ CP (true) for which a discovery is possible for CP-violation is 41% (30%) for DUNE (LBNO) while in the standard case, the coverage is 48% for DUNE and 55% for LBNO. In case of true IH, the impact of long-range potential is even more striking for these future facilities. As an example, if α eµ (true) = 10 −52 , their chances of establishing CP-violation are quite minimal: only for 37% (12%) values of δ CP (true), DUNE (LBNO) can reject both the CP-conserving values 0 • and 180 • in the fit at 3σ, while in the 'SM' framework, DUNE (LBNO) can do so for 53% (60%) values of true δ CP . As the true value of α eµ approaches toward 10 −52 , the coverages in δ CP (true) for which CP-violation can be observed, diminish very quickly for both the set-ups, and ultimately around α eµ (true) = 2 × 10 −52 , the coverages almost become zero.
Finally, we have asked the question, how robust are mass hierarchy measurements in these future facilities in the presence of LRF? In the standard case, due to the large Earth matter effects at the Fermilab-Homestake and CERN-Pyhäsalmi baselines, both DUNE and LBNO can resolve the issue of mass hierarchy at very high confidence level. Now, if LRF exists in Nature, then for DUNE, the mass hierarchy sensitivity remains above the standard expectations provided the true value of α eµ < 5 × 10 −52 . In case of LBNO, the mass hierarchy discovery reach does not go below the 'SM' value as long as the true value of α eµ is smaller than 10 −52 .
ing the completion of this work. S.S.C. would like to thank the organizers of the XXI DAE-BRNS High Energy Physics Symposium 2014 for giving an opportunity to present the preliminary results of this work.
A Discussion at the Probability Level -Anti-Neutrino Case  Figure 16. The transition probability Pμē as a function of anti-neutrino energy. The band reflects the effect of unknown δ CP . Inside each band, the probability for δ CP = 0 • case is shown by the black dashed line. The left panels (right panels) are for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential. In the upper (lower) panels, we take α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with long-range potential.
In Fig. 16, we plot the exact numerical transition probabilityν µ →ν e as a function of anti-neutrino energy. The band shows the impact of unknown δ CP . Inside each band, the probability for δ CP = 0 • case is shown by the black dashed line. The left panels (right panels) are drawn for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential. In the upper (lower) panels, we consider α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with long-range potential.  Figure 17. The transition probability Pμμ as a function of anti-neutrino energy. The band reflects the effect of unknown δ CP . Inside each band, the probability for δ CP = 0 • case is shown by the black dashed line. The left panels (right panels) are for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential. In the upper (lower) panels, we take α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with long-range potential. Fig. 17 shows the exact numericalν µ →ν µ disappearance probability as a function of anti-neutrino energy. The thin band portrays the mild impact of unknown δ CP . Inside each band, the probability for δ CP = 0 • case is given by the black dashed line. The left panels (right panels) are drawn for 1300 km (2290 km) baseline. In each panel, we compare the probabilities for NH and IH with and without long-range potential. In the upper (lower) panels, we consider α eµ = 10 −52 (α eµ = 10 −51 ) for the cases with long-range potential.