Can INO be sensitive to flavor-dependent long-range forces?

Flavor-dependent long-range leptonic forces mediated by the ultra-light and neutral bosons associated with gauged Le − Lμ or Le − Lτ symmetry constitute a minimal extension of the Standard Model. In presence of these new anomaly free abelian symmetries, the SM remains invariant and renormalizable, and can lead to interesting phenomenological consequences. For an example, the electrons inside the Sun can generate a flavor-dependent long-range potential at the Earth surface, which can enhance νμ and ν¯μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\overline{\nu}}_{\mu } $$\end{document} survival probabilities over a wide range of energies and baselines in atmospheric neutrino experiments. In this paper, we explore in detail the possible impacts of these long-range flavor-diagonal neutral current interactions due to Le − Lμ and Le − Lτ symmetries (one at-a-time) in the context of proposed 50 kt magnetized ICAL detector at INO. Combining the information on muon momentum and hadron energy on an event-by-event basis, ICAL will be sensitive to long-range forces at 90% (3σ) C.L. with 500 kt·yr exposure if the effective gauge coupling αeμ/eτ> 1.2 × 10−53 (1.75 × 10−53). The sensitivity of ICAL towards αeμ (αeτ) is ∼ 46 (53) times better than the existing bound from the Super-Kamiokande experiment at 90% C.L., and at 3σ, the sensitivity of ICAL is comparable to the bounds obtained from the combined solar and KamLAND data.


Introduction and motivation
The confirmation of neutrino flavor oscillation via several pioneering experiments over the last two decades is a landmark achievement in the intensity frontier of the high energy particle physics [1]. All the neutrino oscillation data available so far can be accommodated in the standard three-flavor oscillation picture of neutrinos [2][3][4]. The 3ν mixing framework contains six fundamental parameters: a) three mixing angles (θ 12 , θ 13 , θ 23 ), b) one Dirac CP phase (δ CP ), and c) two independent mass-squared differences 1 (∆m 2 21 and ∆m 2 32 ). Let us briefly discuss about the present status of these oscillation parameters. According to the latest global fit of world neutrino data available till November 2017 [2,5], the best fit values of the solar parameters sin 2 θ 12 and ∆m 2 21 are 0.307 and 7.4×10 −5 eV 2 respectively.

JHEP04(2018)023
we have on the effective gauge couplings α eµ,eτ of the L e − L µ,τ symmetries from various experiments. In section 3, we study in detail how the three-flavor oscillation picture gets modified in presence of long-range potential. We present compact analytical expressions for the effective oscillations parameters in presence of LRF. Next, we show the accuracy of our analytical probability expressions (for L e −L τ ) by comparing them with the exact numerical results. In appendix A, we perform the similar comparison for the L e − L µ symmetry. In section 4, we draw the neutrino oscillograms in (E ν , cos θ ν ) plane for ν e → ν µ and ν µ → ν µ oscillation channels in presence of L e − L µ,τ symmetry. We mention the important features of ICAL detector in section 5. In section 6, we show the expected event spectra in ICAL with and without LRF. Section 7 deals with the simulation procedure that we adopt in this work. Next, we derive the expected sensitivity of ICAL in probing the long-range force parameters α eµ,eτ in section 8, and discuss few other interesting results. Finally, we summarize and draw our conclusions in section 9.
2 Flavor-dependent long-range forces One of the possible ways to extend the SM gauge group SU(3) C ×SU(2) L ×U(1) Y with minimal matter content is by introducing anomaly free U(1) symmetries with the gauge quantum number (for vectorial representations) [30,31] Q = a 0 (B − L) + a 1 (L e − L µ ) + a 2 (L e − L τ ) + a 3 (L µ − L τ ) . (2.1) Here, B and L are baryon and lepton numbers respectively. L l are lepton flavor numbers and a i with i = 0, 1, 2, 3 are arbitrary constants. Note that the SM remains invariant and renormalizable if we extend its gauge group in the above way [32]. There are three lepton flavor combinations: i) L e − L µ (a 1 = 1, a 0,2,3 = 0), ii) L e − L τ (a 2 = 1, a 0,1,3 = 0), and iii) L µ − L τ (a 3 = 1, a 0,1,2 = 0), which can be gauged in an anomaly free way with the particle content of the SM [33][34][35][36]. In this paper, we concentrate on L e − L µ,τ symmetries and the implications of L µ − L τ symmetry in neutrino oscillation will be discussed elsewhere. Over the last two decades, it has been confirmed that neutrinos do oscillate from one flavor to another, which requires that they should have non-degenerate masses and mix among each other [1]. To make it happen, the above mentioned U(1) gauge symmetries have to be broken in Nature [37,38]. It is quite obvious that the resultant gauge boson should couple to matter very weakly to escape direct detection. On top of it, if the extra gauge boson associated with this abelian symmetry is very light, then it can give rise to long-range force having terrestrial range (greater than or equal to the Sun-Earth distance) [37,39,40]. Interestingly, this LRF depends on the leptonic content and the mass of an object. Therefore it violates the universality of free fall which can be tested in the classic lunar ranging [41,42], and Eötvös type gravity experiments [43,44]. Lee and Yang gave this idea long back in ref. [45]. Later, Okun used their idea and gave a 2σ bound on α < 3.4 × 10 −49 (α stands for the strength of long-range potential) for a range of the Sun-Earth distance or more [46,47]. The coupling of the solar electron to L e −L µ/τ gauge boson leads to a flavor-dependent long-range potential for neutrinos [48][49][50], which can affect neutrino oscillations [37-40, 51, 52] in spite of such tight constraint on α as mentioned above. Here, (L e − L µ/τ )-charge of JHEP04(2018)023 ν e is opposite to that of ν µ or ν τ , which results in new non-universal FDNC interactions of neutrinos. These new interactions along with the standard W -exchange interactions between ambient electrons and propagating ν e in matter can alter the effective values of oscillation parameters in non-trivial fashion [53]. For an example, the electrons inside the Sun can generate a flavor-dependent long-range potential V eµ/eτ at the Earth surface which has the following form [37,38], where α eµ/eτ = g 2 eµ/eτ 4π is the "fine structure constant" of the new abelian symmetry and g eµ/eτ is the corresponding gauge coupling. In above equation, N e denotes the total number of electrons (≈ 10 57 ) in the Sun [54] and R SE is the Sun-Earth distance ≈ 1.5 × 10 13 cm = 7.6 × 10 26 GeV −1 . The LRF potential V eµ/eτ in eq. (2.2) comes with a negative sign for antineutrinos and can be probed separately in ICAL along with the corresponding potential for neutrinos. The LRF potential due to the electrons inside the Earth with the Earthradius range (R E ∼ 6400 km) is roughly one order of magnitude smaller as compared to the potential due to the Sun. Therefore, we safely neglect the contributions coming from the Earth [37,38].
There are already tight constraints on the effective gauge coupling α eµ/eτ of L e − L µ/τ abelian symmetry using the data from various neutrino oscillation experiments. In [37], an upper bound of α eµ < 5.5×10 −52 at 90% C.L. was obtained using the atmospheric neutrino data of the Super-Kamiokande experiment. The corresponding limit on α eτ is < 6.4×10 −52 at 90% confidence level. A global fit of the solar neutrino and KamLAND data in the presence of LRF was performed in [38]. They gave an upper bound of α eµ < 3.4 × 10 −53 at 3σ C.L. assuming θ 13 = 0 • . Their limit on α eτ is < 2.5 × 10 −53 at 3σ. In [51], the authors performed a similar analysis to derive the limits on LRF mediated by vector and non-vector (scalar or tensor) neutral bosons assuming one mass scale dominance. A preliminary study to constrain the LRF parameters in the context ICAL detector was carried out in [52]. According to [52], ICAL would be sensitive to α eµ/eτ > 1.65 × 10 −53 at 3σ using an exposure of one Mton·yr and considering only the muon momentum as observable.

Three-flavor neutrino oscillation with long-range forces
In this section, we discuss how the flavor-dependent long-range potential due to the electrons inside the sun modify the oscillation of terrestrial neutrinos. In presence of LRF, the effective Hamiltonian (in the flavor basis) for neutrino propagation inside the Earth is given by

JHEP04(2018)023
where U is the vacuum PMNS matrix [55][56][57], E denotes the energy of neutrino, and V CC represents the Earth matter potential which can be expressed as In above, G F is the Fermi coupling constant, N e is the number density of electron inside the Earth, ρ stands for matter density, and Y e ( Ne Np+Nn ) is the relative electron number density. Here, N p and N n are the proton and neutron densities respectively. For an electrically neutral and isoscalar medium, N e = N p = N n and therefore, Y e = 0.5. In eq. (3.1), ζ, ξ, and η appear due to the long-range potential. In case of L e − L µ symmetry, ζ = −ξ = V eµ with η = 0. On the other hand, if the underlying symmetry is L e − L τ , then ζ = −η = V eτ with ξ = 0. Here, V eµ (V eτ ) is the LRF potential due to the interactions mediated by neutral gauge boson corresponding to L e −L µ (L e −L τ ) symmetry. Since the strength of V eµ/eτ (see eq. (2.2)) does not depend on the Earth matter density, hence its value remains same for all the baselines. In case of antineutrino, the sign of V CC , V eµ , V eτ , and δ CP will be reversed.
It is evident from eq. (3.1) that if the strength of V eµ/eτ is comparable to ∆m 2 31 /2E and V CC , then LRF would certainly affect the neutrino propagation. Now, let us consider some benchmark choices of energies (E) and baselines (L) for which the above mentioned quantities are comparable in the context of ICAL detector. This detector is quite efficient to detect neutrinos and antineutrinos separately in multi-GeV energy range with baselines in the range of 2000 to 8000 km where we have substantial Earth matter effect. Therefore, in table 1, we show the comparison for three choices of E and L: (2 GeV, 2000 km), (5 GeV, 5000 km), and (15 GeV, 8000 km). Using eq. (3.2), we estimate the size of V CC for these three baselines for which the line-averaged constant Earth matter densities (ρ) based on the PREM [58] profile are 3.46 g/cm 3 , 3.9 g/cm 3 , and 4.26 g/cm 3 respectively. From eq. (2.2), we obtain the values of V eµ/eτ for two benchmark choices of α eµ/eτ : 10 −52 and 3 × 10 −53 (see last column of table 1). We compute the value of ∆m 2 31 /2E assuming the best fit value of ∆m 2 31 = 2.524 × 10 −3 eV 2 [2]. Table 1 shows that the quantities ∆m 2 31 /2E, V CC , and V eµ/eτ are of comparable strengths for our benchmark choices of E, L, and α eµ/eτ . It suggests that they can interfere with each other to alter the oscillation probabilities significantly. Next, we study the modification of oscillation parameters in matter in presence of LRF potential.

JHEP04(2018)023
In the above equation, we neglect the off-diagonal terms which are small. Diagonalizing the (2,3) block of H f , we get the following expression for θ m 23 tan 2θ m 23 = cos 2 θ 13 − α cos 2 θ 12 + α sin 2 θ 12 sin 2 θ 13 −W + α sin 2θ 12 sin θ 13 . (3.14) We can obtain the expressions for θ m 13 and θ m 12 by diagonalizing the (1,3) and (1,2) blocks subsequently. These effective mixing angles can be written in following way In the above expressions, λ 3 , λ 2 , and λ 1 take the following forms λ 3 = 1 2 cos 2 θ 13 + α cos 2 θ 12 + α sin 2 θ 12 sin 2 θ 13 − W + (α sin 2θ 12 sin θ 13 − W ) cos 2θ m The eigenvalues m 2 i,m /2E (i = 1, 2, 3) can be written in following fashion and In each panel, we draw the curves for the following three cases: We repeat the same exercise for the effective mass-squared differences 5 in figure 2. From the extreme right panel of figure 1, we can see that θ m 12 approaches to 90 • very rapidly as we increase E. This behavior is true for the SM case and as well as for non-zero α eµ/eτ , but it is not true for θ m 23 and θ m 13 . The long-range potential V eµ/eτ modifies θ m 23 significantly as can be seen from the extreme left panel of figure 1. As we approach to higher energies, θ m 23 deviates from the maximal mixing and its value decreases (increases) very sharply for nonzero α eµ (α eτ ). This opposite behavior in the variation of θ m 23 for finite α eµ and α eτ affect the oscillation probabilities in different manner, which we discuss in next subsection. Note that θ m 23 is independent of V CC (see eq. (3.14)). Therefore, its value remains same for all the baselines and same is true for the SM case as well as for non-zero α eµ/eτ . In case of θ m 13 (see middle panel of figure 1), the impact of V eµ and V eτ are same and its variation is quite different as compared to θ m 23 . Assuming NH, as we go to higher energies, θ m 13 quickly reaches to maximal mixing (resonance point) for both the symmetries as compared to the SM case. Finally, it approaches toward 90 • as we further increase the energy. For α eµ/eτ = 10 −52 , the resonance occurs around 3.5 GeV for 5000 km baseline. An analytical expression for the resonance energy can be obtained from eq. (3.15) assuming θ m 13 = 45 • . In one mass scale dominance approximation (∆m 2 21 = 0, i.e. α = 0), the expression for the resonance energy E res can be obtained from the following: (3.23) 4 In case of non-zero αeτ , we use eq. (3.14), eq. (3.15), and eq. (3.16). For non-zero αeµ, we take the help of eq. 3.16, eq. 3.17, and eq. 3.18 as given in ref. [40]. 5 For non-zero αeτ , we obtain the effective values of ∆m 2 31,m and ∆m 2 21,m using eq. (3.20), eq. (3.21), and eq. (3.22). For finite αeµ, we derive the same using eq. 3.22, eq. 3.23, and eq. 3.24 as given in ref. [40]. Assuming α = 0 in eqs. (3.17) and (3.14), we get a simplified expression of λ 3 which appears as since at E res , the term W 2 is small compared to cos 4 θ 13 , and we can safely neglect it. Comparing eq. (3.24) and eq. (3.23), we obtain a simple and compact expression for E res : Note that in the absence of LRF, the above equation boils down to the well-known expression for E res in the SM case. Also, we notice that the expression for resonance energy is same for both L e − L τ and L e − L µ symmetries (see eq. 3.27 in [40]). It is evident from eq. (3.25) that for a fixed baseline, in the presence of V eµ/eτ , the resonance takes place at lower energy as compared to the SM case (see middle panel of figure 1). We observe from both the panels of figure 2 that in presence of LRF, the variations in ∆m 2 31,m and ∆m 2 21,m with energy are different as compared to the SM case. Interesting to note that both V eµ and V eτ modify the values of effective mass-squared differences in same fashion. In case of ∆m 2 21,m (see right panel of figure 2), it increases with energy and can be comparable to the vacuum value of ∆m 2 31 at around E = 10 GeV for both the SM and SM + LRF scenarios. For ∆m 2 31,m (see left panel of figure 2), the change with energy is very mild in the SM case, but in presence of LRF, ∆m 2 31,m gets increased substantially as we approach to higher energies. In case of antineutrino, the effective values of oscillation parameters can be obtained in the similar fashion by just replacing A → −A and W → −W in eqs. (3.14) to (3.22). Next, we compare the neutrino and antineutrino oscillation probabilities obtained from our analytical expressions with those calculated numerically.  Figure 3. ν e → ν µ (ν e →ν µ ) transition probability for 5000 km in upper left (right) panel assuming NH. In bottom left (right) panel, we show ν µ → ν µ (ν µ →ν µ ) survival probability. In all the panels, we compare our analytical expressions (solid curves) to the exact numerical results (dashed curves) for the SM and SM + LRF cases. For LRF, we consider α eτ = 10 −52 . Note that the y-axis ranges are different in the upper left and right panels.

Comparison between analytical and numerical results
We obtain the analytical probability expressions in the presence of V CC and V eµ/eτ by replacing the well known vacuum values of the elements of U PMNS and the mass-squared differences ∆m 2 ij with their effective values as discussed in the previous section. In figure 3, we show our approximate ν e → ν µ (ν e →ν µ ) oscillation probabilities in the top left (right) panel as a function of E against the exact numerical results considering L = 5000 km 6 and NH. We repeat the same for ν µ → ν µ (ν µ →ν µ ) survival channels in bottom left JHEP04(2018)023 (right) panel. We perform these comparisons among analytical (solid curves) and numerical (dashed curves) cases for both the SM and SM + LRF scenarios assuming our benchmark choice of α eτ = 10 −52 . For L e − L µ symmetry, we perform the similar comparison in figure 12 (see appendix A). For the SM case (α eτ = 0), our approximate results match exactly with numerically obtained probabilities. In the presence of L e − L τ symmetry, we see that our analytical expressions work quite well and can produce almost accurate L/E oscillation patterns.
We can see from the top left panel of figure 3 that for non-zero α eτ , the location of the first oscillation maximum shifts toward lower energy (from 5.8 GeV to 3.5 GeV) and also its amplitude gets enhanced (from 0.18 to 0.64) for ν e → ν µ transition probability assuming NH. To understand this feature, we can use the following simple expression 7 for P (ν e → ν µ ) considering θ m 12 = 90 • (see right panel of figure 1): As can be seen from the previous section, θ m 23 does not "run" for the SM case, but for non-zero α eτ , it approaches toward 90 • as we increase E. As far as θ m 13 is concerned, it quickly reaches to the resonance point at a lower energy for non-zero α eτ as compared to α eτ = 0 case. Also, ∆m 2 32,m (∆m 2 31,m − ∆m 2 21,m ) decreases with energy as ∆m 2 21,m increases substantially in comparison to ∆m 2 31,m till E ∼ 4 GeV for 5000 km baseline. The modifications of mixing parameters in different fashion are responsible to shift the location of first oscillation maximum toward lower energy and also to enhance its amplitude.
In case of ν µ → ν µ survival probability (P µµ ), we can use the following simple expression assuming θ m 12 = 90 • : In the above expression, the term sin 2 2θ m 23 plays an important role. Now, we see from left panel of figure 1 that as we go to higher energies, θ m 23 deviates from the maximal mixing very sharply in presence of LRF. For this reason, the value of sin 2 2θ m 23 gets reduced substantially, which ultimately enhances the survival probability for non-zero α eτ as can be seen from the bottom left panel of figure 3. In the energy range of 6 to 20 GeV, we see a substantial enhancement in P µµ with non-zero α eτ as compared to the SM case. The same is true for non-zero α eµ as can be seen from figure 12 in appendix A. We see a similar increase in case ofν µ →ν µ survival channel with NH (see bottom right panel of figure 3). We observe this behavior for other baselines as well in figures 6 and 8, which we discuss later.     Figure 4. The oscillograms for ν e → ν µ channel in E ν , cos θ ν plane for three different scenarios: i) α eµ = α eτ = 0 (the SM case, left panel), ii) α eµ = 10 −52 , α eτ = 0 (middle panel), and iii) α eµ = 0, α eτ = 10 −52 (right panel). Here, in all the panels, we assume NH.   the neutrino oscillation probabilities for all possible choices of baseline (cos θ ν ) and energy (E ν ) which are relevant for the ICAL detector. We perform this study by drawing the neutrino oscillograms in (E ν , cos θ ν ) plane using the full three-flavor probability expressions with the varying Earth matter densities as given in the PREM profile [58]. Although in atmospheric neutrino experiments, it is not possible to measure the oscillation probabilities for ν e → ν µ and ν µ → ν µ channels separately, but to explain their features from our analytical expressions, here we present the oscillograms for appearance and disappearance probabilities separately.
baselines as compared to the SM case. This feature gets reflected in the middle and right panels of figure 4 for non-zero α eµ and α eτ respectively. Figure 4 also depicts that the value of P eµ decreases (increases) as compared to the SM case for non-zero α eµ (α eτ ). We can explain this behavior from the variation of θ m 23 (see left panel of figure 1). In presence of L e − L µ (L e − L τ ) symmetry, the term sin 2 θ m 23 in eq. (3.26) gets reduced (enhanced) as compared to the SM case, which subsequently decreases (increases) the value of P eµ . In figure 5, we present the oscillograms for ν e → ν µ appearance channel considering a smaller value of α eµ/eτ = 3 × 10 −53 .

Oscillograms for ν µ → ν µ disappearance channel
In figure 6, we present the oscillograms for ν µ survival channel in the plane of cos θ ν vs. E ν considering NH. Here, we draw the oscillograms for the same three cases as considered in figure 4. First, we notice that for E ν in the range of 6 to 20 GeV and cos θ ν in the range of -1 to -0.2, survival probability P µµ gets enhanced significantly for both non-zero α eµ (middle panel) and α eτ (right panel) as compared to the SM case (see left panel). The reason is the following. As we move to higher energies, θ m 23 deviates from maximal mixing for both nonzero α eµ and α eτ . As a result, the term sin 2 2θ m 23 in eq. (3.27) gets reduced and causes an enhancement in P µµ . In figure 6, we see some differences in the oscillogram patterns in the JHEP04(2018)023 energy range of 2 to 5 GeV for L e −L µ (middle panel) and L e −L τ (right panel) symmetries. Let us try to understand the reason behind this. We have already seen that θ m 23 "runs" in the opposite directions from 45 • for L e − L µ and L e − L τ symmetries. Due to this, the only term 1 4 tan 2 θ m 23 sin 2 2θ m 13 in eq. (3.27) gives different contributions for finite α eµ and α eτ . Around the resonance region (E ∼ 2 to 5 GeV), θ m 13 attains the maximal value, and the strength of above mentioned term becomes quite significant which causes the differences in P µµ for these two U(1) symmetries under consideration. We see the effect of this feature in the top left panel of figure 8, which we discuss later. In figure 7, we show the oscillograms for ν µ → ν µ disappearance channel considering a smaller value of α eµ/eτ = 3 × 10 −53 .

Important features of the ICAL detector
The proposed Iron Calorimeter (ICAL) detector [11] under the India-based Neutrino Observatory (INO) [12] project plans to study the fundamental properties of atmospheric neutrino and antineutrino separately using the magnetic field inside the detector. The strength of the magnetic field will be around 1.5 T with a better uniformity in the central region [59]. It helps to determine the charges of µ − and µ + particles which get produced in the charged-current (CC) interactions of ν µ andν µ inside the ICAL detector. To restrict the cosmic muons, which serve as background in our case, the ICAL detector is planned to have rock coverage of more than 1 km all around. According to the latest design of the ICAL detector [11,12], it consists of alternate layers of iron plates and glass Resistive Plate Chambers (RPCs) [60], which act as the target material and active detector elements respectively. While passing through the RPCs, the minimum ionizing particle muon gives rise to a distinct track, whose path is recorded in terms of strip hits. We identify these tracks with the help of a track finding algorithm. Then, we reconstruct the momentum and charge of muon using the well known Kalman Filter [61,62] package. The typical detection efficiency of a 5 GeV muon in ICAL traveling vertically is around 80%, while the achievable charge identification efficiency is more than 95% [63]. In ICAL, the energy resolution (σ/E) of a 5 to 10 GeV muon varies in the range of 10% to 15%, while its direction may be reconstructed with an accuracy of one degree [63]. The prospects of ICAL to measure the three-flavor oscillation parameters based on the observable (E µ , cos θ µ ) have already been studied in refs. [16,17].
The hits in the RPCs due to hadrons produce shower-like features. Recently, the possibilities of detecting hadron shower and measuring its energy in ICAL have been explored [64,65]. These final state hadrons get produced along with the muons in CC deep-inelastic scattering process in multi-GeV energies, and can provide vital information about the initial neutrino. We can calibrate the energy of hadron (E had = E ν − E µ ) using number of hits produced by hadron showers [64]. Preliminary studies have shown that one can achieve an energy resolution of 85% (36%) at 1 GeV (15 GeV). Combining the muon (E µ , cos θ µ ) and hadron (E had ) information on an event-by-event basis, the physics reach of ICAL to the neutrino oscillation parameters can be improved significantly [18]. We follow the refs. [63] and [64] to incorporate the detector response for muons and hadrons respectively.    [11,21] GeV ∈ µ E µ θ cos  In each panel, we consider three different cases: i) α eµ = α eτ = 0 (the SM case, black solid line), ii) α eµ = 10 −52 , α eτ = 0 (blue dash-dotted line), and iii) α eµ = 0, α eτ = 10 −52 (red dashed line). Here, we sum over E had in its entire range of 0 to 25 GeV and show the results for 500 kt·yr exposure and assuming NH.

Event spectrum in the ICAL detector
In this section, we present the expected event spectra and total event rates in ICAL with and without long-range forces. Using the event generator NUANCE [66] and atmospheric neutrino fluxes at Kamioka 8 [68], we obtain the unoscillated event spectra for neutrino and antineutrino. After incorporating the detector response for muons and hadrons as described in ref. [18] and for the benchmark values of the oscillation parameters as mentioned in section 3.1 (sin 2 θ 23 = 0.5, sin 2 2θ 13 = 0.0847, and NH), we obtain around 4870 (2187) µ − (µ + ) events for the SM case using a 500 kt·yr exposure. To obtain these event rates, we consider E µ in the range 1 to 21 GeV, cos θ µ in its entire range of -1 to 1, and E had in the range 0 to 25 GeV. In presence of L e − L µ symmetry with α eµ = 10 −52 , the number of µ − (µ + ) events becomes 5365 (2373). For L e − L τ symmetry with α eτ = 10 −52 , we get 5225 µ − and 2369 µ + events. In figure 8, we show the distribution of only upward going µ − (top panels) and µ + (bottom panels) events as a function of reconstructed cos θ µ in the range -1 to 0. Here, we integrate over the entire range of hadron energy (E had ∈ 0 to 25 GeV), and display the event spectra considering three different E µ bins having the ranges 1 to 5 GeV (left panels),    [11,21] GeV ∈ µ E µ θ cos   [11,21] GeV ∈ µ E Figure 9. The distributions of µ − (upper panels) and µ + (lower panels) events for three different E µ bins: 1 to 5 GeV in left panel, 5 to 11 GeV in middle panel, and 11 to 21 GeV in right panel. In each panel, we consider three different cases: i) α eµ = α eτ = 0 (the SM case, black solid line), ii) α eµ = 3 × 10 −53 , α eτ = 0 (blue dash-dotted line), and iii) α eµ = 0, α eτ = 3 × 10 −53 (red dashed line). Here, we sum over E had in its entire range of 0 to 25 GeV and show the results for 500 kt·yr exposure and assuming NH. 5 to 11 GeV (middle panels), and 11 to 21 GeV (right panels). In each panel, we compare the event distribution for three different scenarios: i) α eµ = α eτ = 0 (the SM case, black solid lines), ii) α eµ = 10 −52 , α eτ = 0 (blue dash-dotted lines), and iii) α eµ = 0, α eτ = 10 −52 (red dashed lines). We observe few interesting features in figure 8, which we discuss now.
In all the panels of figure 8, we see an enhancement in the event rates for cos θ µ ∈ -1 to -0.2 in the presence of long-range forces as compared to the SM case. This mainly happens due to substantial increase in P µµ with finite α eµ or α eτ as compared to the SM case. We have already seen this feature in figure 6. Also, we see similar event distributions for both the symmetries in all the panels, except in the top left panel (E µ ∈ 1 to 5 GeV), where we see some differences in the event spectra for L e − L µ and L e − L τ symmetries. We have already explained the reason behind this with the help of oscillogram patterns (see middle and right panels in figure 6) in section 4.2. In figure 9, we present the distributions of µ − and µ + events in ICAL considering a smaller value of α eµ/eτ = 3 × 10 −53 . Next, we discuss the binning scheme for three observables (E µ , cos θ µ , and E had ), and briefly describe the numerical technique and analysis procedure which we adopt to estimate the physics reach of ICAL.  Table 2. The binning scheme considered for the reconstructed observables E µ , cos θ µ , and E had for each muon polarity. In last column, we give the total number of bins taken for each observable.

JHEP04(2018)023
7 Simulation procedure 7.1 Binning scheme for observables (E µ , cos θ µ , E had ) Table 2 shows the binning scheme that we adopt in our simulation for three observables E µ (∈ 1 to 21 GeV), cos θ µ (∈ -1 to 1), and E had (∈ 0 to 25 GeV). In these ranges, we have total 12 bins for E µ , 15 bins for cos θ µ , and 4 bins for E had , resulting into a total of (12 × 15 × 4 =) 720 bins per polarity. We consider the same binning scheme for µ − and µ + events. As we go to higher energies, the atmospheric neutrino flux decreases resulting in lower statistics. Therefore, we take wider bins for E µ and E had at higher energies. We do not perform any optimization study for binning, however we make sure that we have sufficient statistics in most of the bins without diluting the sensitivity much. In our study, the upward going events (cos θ µ in the range 0 to -1) play an important role, where V CC , V eµ/eτ , and ∆m 2 31 /2E become comparable and can interfere with each other (see discussion in section 3). Therefore, we take 10 bins of equal width for upward going events which is compatible with the angular resolutions of muon achievable in ICAL. The downward going events do not undergo oscillations. But, they certainly enhance the overall statistics and help us to reduce the impact of normalization uncertainties in the atmospheric neutrino fluxes. Therefore, we include the downward going events in our simulation considering five cos θ µ bins of equal width in the range of 0 to 1.

Numerical analysis
In our numerical analysis, we suppress the statistical fluctuations of the "observed" event distribution. We generate 9 events using NUANCE for an exposure of 50000 kt·yr. Then, we implement the detector response and finally, normalize the event distribution to the actual exposure. This method along with the χ 2 function gives us the median sensitivity of the experiment in the frequentist approach [69]. We use the following Poissonian χ 2 − for JHEP04(2018)023 µ − events in our statistical analysis: In the above equation, N data ijk and N theory ijk denote the "observed" and expected number of µ − events in a given (E µ , cos θ µ , E had ) bin. N 0 ijk represents the number of events without systematic uncertainties. In our simulation, N E had = 4, N Eµ = 12, and N cos θµ = 15 (see table 2). We obtain N data ijk using the benchmark values of the oscillation parameters as mentioned in section 3.1 and assuming normal hierarchy as neutrino mass hierarchy. We consider five systematic errors in our analysis: 20% flux normalization error, 10% error in cross-section, 5% tilt error, 5% zenith angle error, and 5% overall systematics. We incorporate these systematic uncertainties in our simulation using the well known "pull" method [70][71][72].
In a similar fashion, we obtain χ 2 + for µ + events. We estimate the total χ 2 by adding the individual contributions coming from µ − and µ + events in the following way In the fit, we first minimize χ 2 ICAL with respect to the pull variables ζ l , and then marginalize over the oscillation parameters sin 2 θ 23 in the range 0.38 to 0.63 and ∆m 2 31 in the range 0.0024 eV 2 to 0.0026 eV 2 . While deriving the sensitivities to α eµ/eτ , we also marginalize χ 2 ICAL over both NH and IH. We do not marginalize over ∆m 2 21 , sin 2 θ 12 , and sin 2 2θ 13 since these parameters are already measured with high precision, and the existing uncertainties on these parameters do not alter our results. We consider δ CP = 0 • throughout our analysis.

Results
We quantify the statistical significance of the analysis to obtain the sensitivity of ICAL towards the LRF parameters in the following way Here, χ 2 ICAL (SM) and χ 2 ICAL SM + α eµ/eτ are calculated by fitting the "observed" data in the absence and presence of LRF parameters respectively. In our analysis, statistical fluctuations are suppressed, and therefore, χ 2 ICAL (SM) ≈ 0. Before we present the constraints on α eµ/eτ , we identify the regions in E µ and cos θ µ plane which give significant contributions toward ∆χ 2 ICAL−LRF . In figure 10, we show the distribution 10 of ∆χ 2 µ − (left 10 In figure 10, we do not consider the constant contributions in χ 2 coming from the term which involves five pull parameters ζ 2 l in eq. (7.1). Also, we do not marginalize over the oscillation parameters in the fit to produce these figures. But, we show our final results considering full pull contributions and marginalizing over the oscillation parameters in the fit as mentioned in previous section.     Figure 10. Distributions of ∆χ 2 ICAL−LRF (per unit area) in E µ and cos θ µ plane. The left (right) panels are for µ − (µ + ) events. In upper (lower) panels, we assume non-zero α eµ (α eτ ) in the fit with a strength of 10 −52 . In all the panels, we use 500 kt·yr exposure and assume NH in both data and theory. panels) and ∆χ 2 µ + (right panels) in the reconstructed E µ and cos θ µ plane, where the events are further divided into four sub-bins depending on the reconstructed hadron energy (see table 2). In the upper (lower) panels of figure 10, we take non-zero α eµ (α eτ ) in the fit with a strength of 10 −52 . We clearly see from the left panels that for µ − events, most of the contributions (∼ 70%) stem from the range 6 to 15 GeV for E µ and for cos θ µ , the effective range is -0.8 to -0.4. We see similar trend for both the symmetries (see upper and lower panels) and for µ + events (see right panels) as well. Figure 11 shows the sensitivity of ICAL to α eµ and α eτ (one at-a-time) using 500 kt·yr exposure assuming no signal of long-range forces in the data. We find new sensitivity results on α eµ or α eτ by generating the data with no long-range forces and fitting it with some non-zero value of α eµ/eτ by means of χ 2 technique as outlined in previous section. The corresponding ∆χ 2 ICAL−LRF obtained after marginalizing over sin 2 θ 23 , ∆m 2 31 , hierarchy, and systematics parameters in the fit, is plotted in figure 11 as a function of α eµ/eτ (test). It gives a measure of the sensitivity reach of ICAL to the effective gauge coupling of LRF. For both the symmetries, we assume NH as true hierarchy. We obtain similar sensitivities for both the symmetries (one at-a-time) since α eµ and α eτ affect both P µµ and P eµ oscillation channels in almost similar fashion over a wide range of energies and baselines (see figures 4 and 6). ICAL with an exposure of 500 kt·yr would be sensitive to long-range force parameters if α eµ/eτ > 1.2 × 10 −53 (1.75 × 10 −53 ) at 90% (3σ) C.L. assuming NH as true hierarchy. This future sensitivity of ICAL to α eµ is ∼ 46 times better than the existing limit from the Super-Kamiokande experiment at 90% C.L [37]. For α eτ , the sensitivity is 53 times better at 90% confidence level. We obtain similar sensitivity assuming IH as true hierarchy. We see a marginal improvement in the sensitivity if we keep all the oscillation parameters fixed in the fit. In this fixed parameter case, ICAL would be sensitive to α eµ > 1.63 × 10 −53 at 3σ confidence level. We study few interesting issues in this fixed parameter scenario which we discuss now.

JHEP04(2018)023
• Advantage of spectral information: in ICAL, we can bin the atmospheric neutrino/antineutrino events in the observables E µ , cos θ µ , and E had . It helps us immensely to achieve hierarchy measurement at around 3σ C.L. with 500 kt·yr exposure [18]. We find that the ability of using the spectral information in ICAL also plays an important role to place tight constraint on LRF parameters. For an example, if we rely only on the total µ − and µ + event rates, the expected sensitivity from ICAL becomes α eµ > 2.2 × 10 −52 at 3σ confidence level. This sensitivity is almost 13 times weaker as compared to what we can obtain using the full spectral informations.
• Usefulness of hadron energy information: in our analysis, we use the hadron energy information (E had ) along with the muon momentum (E µ , cos θ µ ). We observe that with a value of α eµ = 1.63 × 10 −53 in the fit, ∆χ 2 ICAL−LRF increases from 5.2 to 9 when we use E µ , cos θ µ , and E had as our observables instead of only E µ and cos θ µ . It corresponds to about 73% improvement in the sensitivity.

JHEP04(2018)023
• The role of charge identification capability: we also find that the charge identification capability of ICAL in distinguishing µ − and µ + events does not play an important role to make ICAL sensitive to the LRF parameters unlike the mass hierarchy measurements. Since the long-range forces affect the µ − and µ + event rates in almost similar fashion as compared to the SM case (see Fig 8), it is not crucial to separate these events in our analysis to obtain the sensitivity of ICAL to the LRF parameters.
Before we summarize and draw our conclusions in the next section, we make few comments on how the presence of LRF parameters may affect the mass hierarchy measurement in ICAL. To perform this study, we generate the data with a given hierarchy and assuming α eµ = α eτ = 0. Then, while fitting the "observed" event spectrum with the opposite hierarchy, we introduce α eµ or α eτ (one at-a-time) in the fit and marginalize over it in the range of 10 −55 to 10 −52 along with other oscillation parameters. During this analysis, we find that the mass hierarchy sensitivity of ICAL gets reduced very marginally by around 5%.

Summary and conclusions
The main goal of the proposed ICAL experiment at INO is to measure the neutrino mass hierarchy by observing the atmospheric neutrinos and antineutrinos separately and making use of the Earth matter effects on their oscillations. Apart from this, ICAL detector can play an important role to unravel various new physics scenarios beyond the SM (see refs. [22][23][24][25][26][27][28][29]). In this paper, we have studied in detail the sensitivity of ICAL in probing the 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. It constitutes a minimal extension of the SM preserving its renormalizibility and may alter the expected event spectrum in ICAL. As an example, the electrons inside the sun can generate a flavordependent long-range potential V eµ/eτ at the Earth surface, which may affect the effective values of oscillation parameters in presence of the Earth matter. Important point to note here is that for atmospheric 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µ/eτ even for α eµ/eτ ∼ 10 −52 , and can influence the atmospheric neutrino experiments significantly. Also, for a wide range of baselines accessible in atmospheric neutrino experiments, the Earth matter potentials (V CC ) are around 10 −13 eV (see table 1), suggesting that V CC can interfere with V eµ/eτ and ∆m 2 31 /2E, and can modify the oscillation probability substantially. In this article, we have explored these interesting possibilities in the context of the ICAL detector.
After deriving approximate analytical expressions for the effective neutrino oscillation parameters in presence of V CC and V eµ/eτ , we compare the oscillation probabilities obtained using our analytical expressions with those calculated numerically. Then, we have studied the impact of long-range forces by drawing the neutrino oscillograms in E ν and cos θ ν plane using the full three-flavor probability expressions with the varying Earth matter densities based on the PREM profile [58]. We have also presented the expected event spectra and total event rates in ICAL with and without long-range forces. As non-zero α eµ and α eτ can change the standard 3ν oscillation picture of ICAL significantly, we can expect to have JHEP04(2018)023 good sensitivity of ICAL towards the long-range force parameters. The ICAL detector will be sensitive to long-range forces at 90% (3σ) C.L. with an exposure of 500 kt·yr if the effective gauge coupling α eµ/eτ turns out to be greater than 1.2 × 10 −53 (1.75 × 10 −53 ) and NH is the true hierarchy in Nature. If there is no long-range force in Nature, the expected sensitivity from ICAL to α eµ (α eτ ) is ∼ 46 (53) times better than the existing limit from the Super-Kamiokande experiment at 90% confidence level. At 3σ, the expected sensitivity of ICAL to α eµ/eτ is comparable to the limits which have been obtained using the combined solar and KamLAND data.
In ref. [52], a preliminary study was performed to explore the sensitivity of ICAL towards the LRF parameters. In that study, only muon energy (E µ ) and muon directions (cos θ µ ) were used as observables and considered muon energy range was 0.8 GeV to 15 GeV. With an exposure of one Mton·yr, ICAl was found to be sensitive to α eµ/eτ > 1.65×10 −53 at 3σ. In our paper, we have also included the information on the reconstructed hadron energy (E had ) along with the most recent information on ICAL's response towards the muon observables. We have also considered a slightly larger muon energy range of 1 GeV to 21 GeV. Mainly, due to the inclusion of hadron energy information and larger muon energy range, we obtain almost similar sensitivity towards the LRF parameters as obtained in ref. [52], but using only 500 kt·yr exposure, which is half of the exposure that was considered in ref. [52].
At the end, we would like to mention that if the range of LRF is equal or larger than our distance from the Galactic Center, then the collective long-range potential due to all the electrons inside the Galaxy needs to be taken into account [38]. In such cases, ICAL can be sensitive to even lower values of α eµ/eτ . We hope that our present work can be an important addition to the series of interesting physics studies which can be performed using the proposed ICAL detector at the India-based Neutrino Observatory.  Figure 12. ν e → ν µ (ν e →ν µ ) transition probability for 5000 km in upper left (right) panel assuming NH. In bottom left (right) panel, we show ν µ → ν µ (ν µ →ν µ ) survival probability. In all the panels, we compare our analytical expressions (solid curves) to the exact numerical results (dashed curves) for the SM and SM + LRF cases. For LRF, we consider α eµ = 10 −52 . Note that the y-axis ranges are different in the upper left and right panels.
(dashed curves) cases for both the SM and SM + LRF scenarios considering our benchmark choice of α eµ = 10 −52 . For the SM case (α eµ = 0), the approximate results match exactly with numerically obtained probabilities. Analytical expressions also work quite well in presence of L e − L µ symmetry, and can produce almost accurate L/E oscillation patterns.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.