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 $L_e-L_\mu$ or $L_e-L_\tau$ 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 $\nu_\mu$ and $\bar\nu_\mu$ 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 $L_e-L_\mu$ and $L_e-L_\tau$ 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 can place stringent constraints on the effective gauge coupling $\alpha_{e\mu/e\tau}<1.2\times 10^{-53}$ ($1.75\times 10^{-53}$) at 90$\%$ (3$\sigma$) C.L. with 500 kt$\cdot$yr exposure. The 90$\%$ C.L. limit on $\alpha_{e\mu}$ ($\alpha_{e\tau}$) from ICAL is $\sim 46$ (53) times better than the existing bound from the Super-Kamiokande experiment.


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], strength of long-range potential of V eµ/eτ symmetry at the Earth surface generated by the electrons inside the Sun. We end this section by mentioning the current constraints that 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 constraints on α eµ,eτ from ICAL in section 8, and discuss few other interesting results. Finally, we summarize and draw our conclusions in section 9.

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 longrange force having terrestrial range (greater than or equal to the Sun-Earth distance) and without introducing extremely low mass scales [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 ν 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 "running" 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]. Using an exposure of one Mton·yr and considering only the muon momentum as observable, an expected upper bound of α eµ/eτ 1.65 × 10 −53 at 3σ was obtained for ICAL.

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 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 (3.2) 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 underline 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 5 × 10 −52 (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 "running" of oscillation parameters in matter in presence of LRF potential.

"Running" of Oscillation Parameters
The approximate analytical expressions for the effective mass-squared differences and mixing angles in presence of V CC and V eµ (due to L e − L µ symmetry) have been given in Ref. [40]. In this paper, we derive the analytical expressions for L e − L τ symmetry. Assuming δ CP = 0 • , the effective Hamiltonian can be written as  , and V eµ/eτ (fifth column) for our benchmark choices of E, L, and α eµ/eτ . We take ∆m 2 31 = 2.524 × 10 −3 eV 2 . Based on the PREM profile, the line-averaged constant Earth matter densities for 2000 km, 5000 km, and 8000 km baselines are 3.46 g/cm 3 , 3.9 g/cm 3 , and 4.26 g/cm 3 respectively. The parameter θν is the zenith angle for a given baseline. symmetry, V = Diag(V CC + V eτ , 0, −V eτ ). Considering maximal mixing for θ 23 (= 45 • ), we rewrite H f in the following way where b 11 = A + W + sin 2 θ 13 + α sin 2 θ 12 cos 2 θ 13 , (3.5) b 12 = 1 √ 2 cos θ 13 (α cos θ 12 sin θ 12 + sin θ 13 − α sin 2 θ 12 sin θ 13 ) , (3.6) b 13 = 1 √ 2 cos θ 13 (−α cos θ 12 sin θ 12 + sin θ 13 − α sin 2 θ 12 sin θ 13 ) , (3.7) b 22 = 1 2 cos 2 θ 13 + α cos 2 θ 12 − α sin 2θ 12 sin θ 13 + α sin 2 θ 12 sin 2 θ 13 , (3.8) b 23 = 1 2 cos 2 θ 13 − α cos 2 θ 12 + α sin 2 θ 12 sin 2 θ 13 , (3.9) b 33 = 1 2 cos 2 θ 13 + α cos 2 θ 12 + α sin 2θ 12 sin θ 13 + α sin 2 θ 12 sin 2 θ 13 − 2W . (3.10) In the above equations, the terms A, W , and α are defined as , and α ≡ ∆m 2 21 ∆m 2 31 . (3.11) The following unitary matrixŨ can almost diagonalize the effective Hamiltonian (H f ): 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 (3. 16) 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 and The eigenvalues m 2 i,m /2E (i = 1, 2, 3) can be written in following fashion and (3.22) To observe the "running" of oscillation parameters in presence of V CC and V eµ/eτ , we take the following benchmark values of vacuum oscillation parameters: sin 2 θ 23 = 0.5,  In each panel, we draw the curves for the following three cases 4 : i) α eµ = α eτ = 0 (the SM case), ii) α eµ = 10 −52 , α eτ = 0 iii) α eµ = 0, α eτ = 10 −52 . We repeat the same exercise for the effective mass-squared differences 5 in Fig. 2. From the extreme right panel of Fig. 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τ affect the "running" of θ m 23 significantly as can be seen from the extreme left panel of Fig. 1. As we approach to higher energies, θ m 23 deviates from the maximal mixing and its value decreases (increases) very sharply for non-zero α eµ (α eτ ). This opposite behavior in "running" 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 Fig. 1), the impact of V eµ and V eτ are same and its "running" 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 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 running 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]. E res can be obtained from the following: (3.23) 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]  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 Fig. 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 Fig. 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 "running" 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.

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 "running" values as discussed in the previous section. In Fig. 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 (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 Fig. 9 (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 Fig. 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 Fig. 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. All these different "running" of oscillation parameters are responsible to shift the location of first oscillation maximum toward lower energy and also to enhance its amplitude. 6 For both analytical and numerical calculations, we take the line-averaged constant Earth matter density based on the PREM profile [58]. 7 We obtain this formula using the general expression as given in Eq. 3.30 in Ref. [40]. 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 Fig. 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 Fig. 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 Fig. 9 in appendix A. We see a similar increase in case ofν µ →ν µ survival channel with NH (see bottom right panel of Fig. 3). We observe this behavior for other baselines as well in Figs. 5 and 6, which we discuss later.
4 Neutrino Oscillograms in (E ν , cos θ ν ) Plane The atmospheric neutrino experiments deal with a wide range of baselines and energies. Therefore, it is quite important to see how the long-range forces under discussion affect 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]. Fig. 4 shows the oscillograms for ν e to ν µ appearance channel in E ν and cos θ ν plane assuming NH. We present the oscillograms for three different cases: i) extreme left panel is for the SM case (α eµ = α eτ = 0), ii) middle panel is for the SM + LRF (α eµ = 10 −52 ), and iii) extreme right panel deals with the SM + LRF (α eτ = 10 −52 ). For the SM case, ν e to ν µ transition probability attains the maximum value around the resonance region which occurs in the range of E ∈ 4 to 8 GeV and cos θ ν ∈ -0.8 to -0.4. The resonance condition in presence of LRF (see Eq. 3.25) suggests that θ m 13 can reach 45 • at smaller energies and baselines as compared to the SM case. This feature gets reflected in the middle and right panels of Fig. 4 for non-zero α eµ and α eτ respectively. Fig. 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 "running" of θ m 23 (see left panel of Fig. 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µ .

Oscillograms for ν µ → ν µ Disappearance Channel
In Fig. 5, 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 Fig. 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 non-zero α 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 Fig. 5, we see some differences in the oscillogram patterns in ν θ cos Here, in all the panels, we assume NH.   the 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 Fig. 6, which we discuss later.

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 Ref. [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.

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 8 Preliminary calculation of the expected fluxes at the INO site have been performed in Ref. [67]. We plan to use these fluxes in future analysis once they are finalized. The horizontal components of the geo-magnetic field are different at the INO (40 µT) and Kamioka (30 µT). Due to this reason, we observe a difference in atmospheric fluxes at these sites.    [11,21] GeV ∈ µ E µ θ cos  µ − and 2369 µ + events. In Fig. 6, 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), 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 Fig. 6, which we discuss now.
In all the panels of Fig. 6, 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 Fig. 5. 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 Fig. 5) in section 4.2. 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 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 µ − 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 constraints on α 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 constrain 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 Fig. 7, we show the distribution 10 of ∆χ 2 µ − (left panels) and 10 In Fig. 7, 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. µ θ cos      (see table 2). In the upper (lower) panels of Fig. 7, 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. Fig. 8 shows the upper bound on α eµ and α eτ (one at-a-time) using 500 kt·yr exposure of ICAL if there is no signal of long-range forces in the data. We set new upper limit 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 Fig. 8 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 constraints 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 Figs. 4 (fit) α and 5). The expected upper limit on α eµ/eτ from ICAL is < 1.2 × 10 −53 (1.75 × 10 −53 ) at 90% (3σ) C.L. with 500 kt·yr exposure and NH as true hierarchy. This future limit on α eµ from ICAL at 90% C.L. is ∼ 46 times better than the existing limit from the Super-Kamiokande experiment [37]. For α eτ , the limit is 53 times better at 90% confidence level. We obtain similar constraints assuming IH as true hierarchy. We see a marginal improvement in the upper limits if we keep all the oscillation parameters fixed in the fit. In this fixed parameter case, the new bound becomes α eµ < 1.63 × 10 −53 at 3σ confidence level. We study few interesting issues in this fixed parameter scenario which we discuss now.
• 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 limit from ICAL becomes α eµ < 2.2 × 10 −52 at 3σ confidence level. This limit 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.
• 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 constrain 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 6), it is not crucial to separate these events in our analysis in constraining 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 capabilities of ICAL to constrain 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 running 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 (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 place strong limits on these parameters if ICAL do not observe a signal of LRF in oscillations. The expected upper bound on α eµ/eτ from ICAL is < 1.2 × 10 −53 (1.75 × 10 −53 ) at 90% (3σ) C.L. with 500 kt·yr exposure and NH as true hierarchy. ICAL's limit at 90% C.L. on α eµ (α eτ ) is ∼ 46 (53) times better than the existing limit from the Super-Kamiokande experiment. Here, 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.

Acknowledgment
This work is a part of the ongoing effort of INO-ICAL Collaboration to study various physics potentials of the proposed ICAL detector. Many members of the Collaboration have contributed for the completion of this work. We would like to thank A. Dighe, A.M. Srivastava, P. Agrawal, S. Goswami, P.K. Behera, D. Indumathi for their useful comments on our work. We acknowledge S.S. Chatterjee for many useful discussions during this work. We thank A. Kumar  A Oscillation Probability with L e − L µ Symmetry Fig. 9 shows 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 and NH. We repeat the same for ν µ → ν µ (ν µ →ν µ ) survival channels in bottom left (right) panel. We perform these comparisons among analytical (solid curves) and numerical (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.