Probing Long-Range Forces in Neutrino Oscillations at the ESSnuSB Experiment

(ESSnuSB Collaboration)

The discovery of neutrino oscillations Super-Kamiokande:1998kpq ; SNO:2002tuh ; Kajita:2016cak ; McDonald:2016ixn has provided compelling evidence for physics beyond the Standard Model (SM), opening new avenues to explore new fundamental interactions and forces. The unique properties of neutrinos, including their elusive nature and tiny masses, make them an excellent probe to detect even the most subtle signatures of new physics. The neutrino experiments with their increasing precision are now sensitive to sub-leading effects due to potential non-standard interactions (NSIs), offering an indirect hint of new particles and forces not predicted by the SM.

In the standard scenario, the interaction of neutrinos with matter is described by the so-called Mikheyev–Smirnov–Wolfenstein (MSW) mechanism, which results from coherent forward scattering of neutrinos with ambient matter Wolfenstein:1977ue . In this seminal paper, Wolfenstein also proposed the possibility of NSIs¹¹1In this manuscript, we will focus on new interactions mediated by vector bosons. There exists other forms of such interactions with different Lorentz structures Kopp:2007ne ; Du:2020dwr ; Gupta:2023wct ; ESSnuSB:2023lbg ; Denton:2024upc , which physics signatures are however different from the ones discussed here., which have been extensively studied in the literature Ohlsson:2012kf ; Miranda:2015dra ; Farzan:2017xzy ; 10.21468/SciPostPhysProc.2.001 ; Huitu:2016bmb ; Chaves:2021kxe .

In this work, we focus on another kind of such a new leptonic neutrino-matter interaction known as the long-range force (LRF), which may be flavour-dependent and mediated by light vector mediators He:1990pn ; Foot:1990mn ; He:1991qd ; Foot:1994vd ; Dolgov:1999gk . This is particularly intriguing since their effects can accumulate over astronomical distances, making them distinct from other NSIs. For instance, the matter content within astrophysical objects (Sun, Earth, Milky Way, etc.) can act as a source of LRF potential. These interactions significantly modify the probabilities of neutrino oscillations by introducing new potential terms in the Hamiltonian for neutrino propagation Smirnov:2019cae . Such interactions originate by extending the SM gauge group with additional anomaly-free $U(1)$ symmetries associated with lepton numbers $L_{e},L_{\mu},L_{\tau}$ and the baryon number $B$. We consider the three possible combinations of lepton flavours of symmetries Pontecorvo:1967fh ; Gribov:1968kq ; Cirigliano:2005ck ; Altarelli:2010gt , for example, $L_{e}-L_{\mu}$, $L_{e}-L_{\tau}$, $L_{\mu}-L_{\tau}$. These symmetries are also important for generating neutrino masses Asai:2017ryy ; Asai:2018ocx ; Lou:2024fvw . The constraints on the LRF parameters have already been obtained from solar Grifols:2003gy ; Bandyopadhyay:2006uh ; Gonzalez-Garcia:2006vic , atmospheric Joshipura:2003jh and astrophysical neutrinos Bustamante:2018mzu ; Agarwalla:2023sng ²²2Most works on LRFs in neutrino oscillations take into account specific models or mediator mass ranges. The only model-independent constraints on the LRF potentials from existing experiments come from high-energetic neutrinos observed at IceCube Agarwalla:2023sng and of $\mathcal{O}(10^{-19}\,\,\mathrm{eV})$. These bounds are much tighter than the ones expected at terrestrial experiments due to energy-enhanced effects of LRFs. However, it is still worth exploring the bounds of accelerator experiments that employ a well-known and controlled neutrino beam.. In Ref. Coloma:2020gfv , a global analysis of three-flavour oscillation data has been performed in the presence of flavour-dependent long-range interactions. Furthermore, the effect of LRFs on long-baseline (LBL) neutrino experiments has been explored in Refs. Chatterjee:2015gta ; Khatun:2018lzs ; Singh:2023nek ; Mishra:2024riq .

A key objective of present and future neutrino oscillation experiments T2K:2017hed ; T2K:2019bcf ; NOvA:2019cyt ; DUNE:2020jqi is the precise determination of the leptonic CP-violating phase $\delta_{\rm CP}$. The European Spallation Source (ESS) neutrino Super-Beam ESSnuSB ESSnuSB:2021azq is a next-to-next-generation long-baseline neutrino oscillation experiment designed to achieve this goal. Located in Sweden, ESSnuSB will produce a high-intensity muon neutrino beam using a 5 MW proton beam from the upgraded ESS facility in Lund Abele:2022iml ; Alekou:2022emd . The neutrinos will be detected by a water-Cherenkov detector situated 360 km away from Lund at the mine in Zinkgruvan. By focusing on the second oscillation maximum in the appearance probability $P_{\mu e}$, ESSnuSB is uniquely positioned to provide a precise measurement of $\delta_{\rm CP}$. Currently, ESSnuSB is at the stage of preparation of a second conceptual design report to be followed by the development of a technical design report. This will help to plan the construction and data collection at a later stage ESSnuSB:2023ogw ; ESSnuSB:2024tmn . In the present work, we perform the first comprehensive study of the impact of long-range forces on the physics sensitivities of the ESSnuSB experiment. We derive bounds on the LRF potentials and the associated coupling parameters, comparing them with those achievable in the next-generation LBL experiments DUNE and T2HK. In addition, we investigate the effects of LRFs on the measurement of $\delta_{\rm CP}$ by ESSnuSB. Our analysis demonstrates that ESSnuSB’s long baseline and high precision make it an ideal facility for probing the subtle effects of LRFs, offering sensitivity that surpasses those of some existing experiments.

This paper is organized as follows. In Section 2, we provide a brief overview of the theoretical framework of LRFs in neutrino oscillations, focusing on the three $U(1)$ symmetries under consideration. Then, in Section 3, the description of the ESSnuSB experiment and other simulation details are provided. Next, in Section 4, we compute the transition probabilities and generate the event plots in the presence of LRFs for ESSnuSB. In Section 5, the sensitivity of the ESSnuSB experiment to constrain the LRF potentials and new coupling parameters are presented. Especially, Subsection 5.1 deals with some interesting correlations of LRF potentials with $\theta_{23}$ and $\delta_{\rm CP}$. Furthermore, in Section 6, the impact of LRFs on the measurement of $\delta_{\rm CP}$ is discussed, which is followed by a precision study of CP violation (CPV) in Section 7. Finally, in Section 8, we summarize our findings and conclusions.

Neutrino flavour transitions are significantly influenced by the interactions between neutrinos and the ambient matter as they propagate from the source to the detector. These interactions induce an effective potential in the Hamiltonian interaction Wolfenstein:1977ue . In standard scenario, neutrino-matter interactions occur through Charged Current (CC) and Neutral Current (NC) mechanisms. While standard NC interactions are flavour-universal and do not impact neutrino oscillations, possible Beyond Standard Model (BSM) neutrino-matter interactions could introduce new potential terms that significantly alter neutrino propagation. Long-range forces are one such case, which may affect the measurements of neutrino oscillations in long-baseline experiments.

To generate the probability spectrum, analyze event rates, and perform sensitivity studies of ESSnuSB in the presence of LRFs, we employed the GLoBES software Huber:2004ka ; Huber:2007ji . We introduced modifications to the probability engine to incorporate new potential terms due to the LRF as a new physics effect and then carried out numerical computations to obtain event rates and $\chi^{2}$ values. The experimental configuration and parameters for ESSnuSB used in our study are based on the ESSnuSB Conceptual Design Report Alekou:2022emd and were therefore implemented in GLoBES.

In particular, we considered a water Cherenkov far detector with a fiducial volume of 538 kt, positioned in the mine at Zinkgruvan, 360 km away from the neutrino source at ESS in Lund. A powerful linear accelerator (linac) will deliver $2.7\times 10^{23}$ protons on target per year, with a beam power of 5 MW and a proton kinetic energy of 2.5 GeV. We adopted updated neutrino fluxes, peaking at approximately 0.25 GeV, and applied updated migration matrices for event selection, as outlined in Refs. Alekou:2022emd ; ESSnuSB:2024yji . The energy spectrum in the [0, 2.5] GeV range was divided into 50 bins for event calculations.

Our analyses included both the appearance ($\nu_{\mu}\to\nu_{e}$) and disappearance ($\nu_{\mu}\to\nu_{\mu}$) channels and their CP-conjugate transitions, and accounted for all the relevant backgrounds. We assumed systematic errors of $5\%$ for signals and $10\%$ for backgrounds unless otherwise stated. The total exposure time assumed for the far detector is 10 years, equally divided between 5 years of running the neutrino beam and 5 years for the antineutrino beam.

In this section, we first examine how the appearance and disappearance oscillations probabilities of muon neutrinos are influenced by the presence of a new interaction potential, $V_{\alpha\beta}$, sourcing the LRF at the ESSnuSB energies. Subsequently, we analyze the expected total event rates under the inclusion of LRFs in the theoretical framework. Unless stated otherwise, we adopt the best-fit values for the standard oscillation parameters from NuFIT 5.2 Esteban:2020cvm ; NuFIT5.2 , which incorporate Super-Kamiokande atmospheric data and these parameters are summarized in Table 1. For this analysis, we focus solely on the normal mass ordering (NO) for neutrinos, in line with the global fit preference for NO Capozzi:2017ipn ; Esteban:2018azc ; deSalas:2020pgw ; Capozzi:2021fjo ; Gonzalez-Garcia:2021dve , which might also be suggested by recent DESI-BAO cosmological results Jiang:2024viw .

In this section, we explore the capability of the ESSnuSB experiment to constrain the parameters of LRFs. The statistical analysis has been performed using a Poissonian $\chi^{2}$ function, defined as

where $\vec{\Lambda}$ represents the set of oscillation parameters needed to compute the rates, $\sigma_{b}$ is the normalization error, $n$ is the number of energy bins, $O_{i}$ are the observed rates and $E_{i}$ are the expected rates used for the fit. Systematic uncertainties are incorporated using the pull method Huber:2002mx ; Fogli:2002pt , implemented in GLoBES with the nuisance parameter $b$. The significance of our results in terms of standard deviations ($\sigma$) has been obtained assuming the Wilk’s theorem Wilks:1938dza ; for instance, for 1 d.o.f $\#\sigma=\sqrt{\Delta\chi^{2}}$.

In order to compute the bounds on $V_{\alpha\beta}$, we generate the true event spectrum using the hypothesis of no LRFs (i.e. $V_{\alpha\beta}$ (true) = 0) corresponding to standard three-flavour neutrino oscillations and fit the true data using the probabilities in the presence of LRFs. It should be noted that, while fitting the true data, only one LRF parameter is considered at a time in the test. This approach is justified as different potentials stem from distinct symmetries, independently affecting the oscillations. In all three cases of symmetries, we vary the potentials, $V_{\alpha\beta}$ from $10^{-15}$ eV to $10^{-13}$ eV in the test. The marginalization has been performed over $\theta_{13},\theta_{23}$ and $|\Delta m^{2}_{31}|$ by varying them within the uncertainty ranges reported in Table 1, while $\delta_{\rm CP}$ is scanned over its full $[-180\degree,180\degree]$ range. We keep the two oscillation parameters, $\theta_{12}$ and $\Delta m^{2}_{21}$, fixed at their best-fit values Esteban:2020cvm . The results are displayed in Fig. 4 where the one-dimensional $\Delta\chi^{2}$ is plotted as a function of LRF potentials $V_{\alpha\beta}$. The upper left (right) plot of Fig. 4 gives the bound on the LRF parameter $V_{e\mu}\leavevmode\nobreak\ (V_{e\tau})$ while the lower plot displays the constraint on $V_{\mu\tau}$. We also show the results for different values of the normalization systematic uncertainty, namely 2% (red curves), 5% (blue curves) and 10% (green curves). The 3$\sigma$ and 90% C.L. bounds are summarized in Table 2 for the standard 5% systematics case along with the 2% and 10% systematics cases. The main results are that ESSnuSB in the nominal conditions (i.e. 5% systematics) may be able to set the 90% limits on $V_{e\mu}<2.99\times 10^{-14}$ eV, $V_{e\tau}<2.05\times 10^{-14}$ eV and $V_{\mu\tau}<1.81\times 10^{-14}$ eV. Notably, a change between $10\%$ and $20\%$ in the bounds can also be observed by variations in systematic uncertainties, particularly when $V_{e\mu}$ and $V_{e\tau}$ are considered. The effect of systematics on $V_{\mu\tau}$ is, on the other hand, less prominent.

Before comparing the ESSnuSB limits on LRF potentials with other experimental bounds, let us try to understand the role of appearance and disappearance channels in constraining $V_{\alpha\beta}$. In Fig. 5, we demonstrate how individual probability channels contribute for the ESSnuSB sensitivity towards LRF potentials, $V_{\alpha\beta}$. It is evident that for the $L_{e}-L_{\mu}$ and $L_{\mu}-L_{\tau}$ symmetries the major sensitivity comes from the disappearance ($P_{\mu\mu}$) probability, whereas the appearance probability ($P_{\mu e}$) plays a major role to place a bound on the $V_{e\tau}$ potential corresponding to the $L_{e}-L_{\tau}$ symmetry. This is also clear from the probability plots presented in Fig. 2, where the effect of the $e\tau$ sector is more visible in the appearance probability ($P_{\mu e}$). Although from Fig. 2, it seems that the $e\mu$ and $\mu\tau$ sectors affect both $P_{\mu e}$ and $P_{\mu\mu}$, however, due to the high statistics of $\nu_{\mu}$ the disappearance event numbers at the far detector, disappearance channel plays an important role in constraining $V_{e\mu}$ and $V_{\mu\tau}$. This explains why different oscillation channels are sensitive to different LRF potentials. In Fig. 5, we notice a dip corresponding to the disappearance-only sensitivity curves (blue solid curves) for all three cases of the LRF potentials when a marginalization over $\theta_{23}$ is performed. Similar features are also observed in other works Singh:2023nek ; Mishra:2024riq . This is because, in the disappearance probability, the octant of $\theta_{23}$ develops a degeneracy with the potential $V_{\alpha\beta}$ picking up the wrong solution in the minimum $\chi^{2}$ calculation when marginalization is performed over $\theta_{23}$. The dip vanishes for the disappearance only case if we fix $\theta_{23}$ to its best-fit value while computing the $\chi^{2}$. Also, the dip disappears when we combine both the appearance and disappearance channels (green solid curves), while marginalizing over $\theta_{23}$, highlighting the importance of the appearance channel, which is less affected by the $\theta_{23}$ octant degeneracy.

In the context of LBL experiments, the most stringent foreseen $90\%$ C.L. limits on LRF potentials have been derived by simulating the future experiment P2SO Mishra:2024riq due to its longer baseline, whereas the bounds (at $90\%$ C.L.) from the simulations of “upcoming” DUNE and T2HK experiments with the standard neutrino flux are given by Singh:2023nek :

It is worth mentioning that with a high-energy neutrino flux, the DUNE bounds on $V_{\alpha\beta}$ might become weaker as shown in Ref. Giarnetti:2024mdt .

Comparing the ESSnuSB results with other expected limits from upcoming LBL experiments in Eq. (5), we find that assuming nominal conditions ($5\%$ systematics), ESSnuSB bounds are less stringent than the DUNE ones by about a factor of 2. This is due to the higher energy and longer baseline for DUNE, so the effect of LRFs is more pronounced. However, ESSnuSB outperforms T2HK by approximately $20\%$. As mentioned earlier, systematic uncertainties play a noticeable role in placing bounds on the LRF potential by ESSnuSB, i.e., achieving a $2\%$ normalization uncertainty could improve the $V_{e\mu}$ and $V_{e\tau}$ constraints, making them comparable to future DUNE bounds. Overall, ESSnuSB is projected to set bounds on LRF parameters that are competitive with those from future LBL experiments such as DUNE and T2HK. Importantly, the complementarity of constraints from various neutrino sources, including accelerator, atmospheric, and solar neutrino data, provides a unique opportunity to significantly narrow the allowed parameter space for LRFs. By combining these results, the interplay between different datasets may uncover synergies that enhance sensitivity to LRF parameters and help elucidate the underlying physics of these new interactions.

In addition to the LRF potentials, $V_{\alpha\beta}$, we also put constraints on the actual parameters of the new neutrino-matter interaction, namely, the mass of the new gauge boson $m_{Z^{\prime}}$ and the effective gauge couplings $G_{\alpha\beta}$. Following the methodology presented in refs. Bustamante:2018mzu ; Singh:2023nek ; Giarnetti:2024mdt , we use Eqs. (4) and (5) to derive the limits on $m_{Z^{\prime}}$ and $G_{\alpha\beta}$. In order to take into account all the matter content present in the Universe, we consider neutrinos from different sources ranging at distances up to $10^{3}$ Gp away from the Earth. This corresponds to the mediator mass, $m_{Z^{\prime}}$ in the range $10^{-10}-10^{-35}$ eV and the LRF potentials originating from all the matter content of the Universe can be rewritten in terms of the contributions from effective potentials relevant at different distances, i.e.,

To find the electron and neutron numbers for the LRF potentials from the Earth, an average density of a continuous distribution is modeled for the Earth such that we get $(N_{e})_{\rm Earth}=(N_{n})_{\rm Earth}\sim 4\times 10^{51}$. The Moon and the Sun are assumed to be point-like electron and neutron sources which correspond to the number of electrons and neutrons as given by $(N_{e})_{\rm Moon}=(N_{n})_{\rm Moon}\sim 5\times 10^{49}$ and $(N_{e})_{\rm Sun}\approx(N_{n})_{\rm Sun}\sim 10^{57}$ Singh:2023nek . In case of the Milky Way, the total matter content can be assumed to be distributed in the form of a thin and a thick disk, a central bulge and a diffuse gas Dehnen:1996fa ; Miller:2013nza ; Bustamante:2018mzu , yielding $(N_{e})_{\rm MW}=(N_{n})_{\rm MW}\sim 10^{67}$. For the cosmological matter content, we use $(N_{e})_{\rm Cosm}\approx(N_{n})_{\rm Cosm}\sim 10^{79}$ adopted from refs. Bustamante:2018mzu ; Singh:2023nek . Utilizing these values of electron and neutron numbers and using Eqs. (4) and (5), the contributing terms of LRF potentials from all sources can be computed provided that the values of $m_{Z^{\prime}}$ and $G_{\alpha\beta}$ are known. To constrain $m_{Z^{\prime}}$ and $G_{\alpha\beta}$, we use the $90\%$ C.L. limits on $V_{\alpha\beta}$ obtained in Table 2 and vary the free parameters. The results are presented in Fig. 6 where red, blue and green curves correspond to the $L_{e}-L_{\mu}$, $L_{e}-L_{\tau}$ and $L_{\mu}-L_{\tau}$ symmetries, respectively. We also show the interaction range $\propto 1/m_{Z^{\prime}}$ on the upper axis of the plot. It is worthwhile to mention that some astrophysical and cosmological phenomena, such as black-hole superradiance Baryakhtar:2017ngi and weak gravity conjecture Arkani-Hamed:2006emk may also exclude some parameter space of the LRFs, providing the non-oscillation exclusion limits. These regions are displayed by the grey bands in Fig. 6. From this figure, one can observe that the most stringent limit comes from the location of the causal horizon, which contains the highest number of electrons and neutrons. Therefore, the LRF potentials experienced by neutrinos from this location will be the largest.

In this section, we examine how the LRF potentials influence the CP-violation sensitivity of the ESSnuSB experiment. This analysis is crucial, as the primary aim of ESSnuSB is to achieve precise measurements of $\delta_{\rm CP}$. It is worth noting that, in the case of maximal CP violation ($\delta_{\rm CP}$ $=\pm 90\degree$), the sensitivity of ESSnuSB can reach up to 12.5$\sigma$ and it can also achieve at least 5$\sigma$ sensitivity for approximately 75% of the other possible values of $\delta_{\rm CP}$ Alekou:2022emd ; ESSnuSB:2023lbg . This surpasses the sensitivity of all upcoming next-generation LBL neutrino oscillation experiments Agarwalla:2022xdo . It is, therefore, vital to determine whether the presence of new physics, such as long-range interactions of neutrinos with matter, could jeopardize this capability or not. To do this, we generate the true event spectrum by varying $\delta_{\rm CP}$ (true) over the full range $[-180\degree,180\degree]$ and compare this with $\delta_{\rm CP}$ $=0\degree$ or $180\degree$ in the test. The same value of LRF potentials $V_{\alpha\beta}$ is considered in both true and test event spectra. The CPV-sensitivity plots are displayed in Fig. 9 in units of $\sqrt{\Delta\chi^{2}}$, where

In this figure, the red curve represents the sensitivity in the standard oscillation scenario ($V_{\alpha\beta}$ $=0)$. The dashed, dotted and dashdot curves are for the potentials $V_{e\mu}$, $V_{e\tau}$ and $V_{\mu\tau}$, respectively. The sensitivity curves plotted in blue colour correspond to the value of $V_{\alpha\beta}$ = $5\times 10^{-14}$ eV, which is comparable to the ESSnuSB constraints, while the curves in green are computed for the LRF potentials, $V_{\alpha\beta}$ = $5\times 10^{-13}$ eV, a much larger potential value than the ESSnuSB bounds. We notice that for small values of the new potentials, the ESSnuSB CP-violation sensitivity remains intact with some negligible impact on its sensitivity around $\delta_{\rm CP}$ $=\pm 90\degree$. However, for large values of the LRF potentials, the $\Delta\chi^{2}$ changes and the positions of the sensitivity maxima are also slightly shifted. To understand this in more detail, we compute the CPV sensitivity as a function of the LRF potentials, $V_{\alpha\beta}$. The results are displayed in Fig. 10 for two choices of $\delta_{\rm CP}$, i.e. $+90\degree$ (solid curves) and $-90\degree$ (dashed curves). We can observe that, for small values of $V_{\alpha\beta}$ ($\lesssim 10^{-14}$ eV), the $\delta_{\rm CP}$ sensitivity of ESSnuSB more or less does not change for all three cases of $V_{\alpha\beta}$. However, with increasing $V_{\alpha\beta}$, the sensitivity decreases, especially for the potential corresponding to the $L_{e}-L_{\tau}$ symmetry. The reason is that when $V_{\alpha\beta}$ are small, they appear as a correction to the standard probability and mildly affect the $\delta_{\rm CP}$ sensitivity, whereas, for large values of $V_{\alpha\beta}$, new resonances might appear, causing a significant drop in the CPV sensitivity of ESSnuSB.

In this section, we will try to understand the impact of LRF potentials on the uncertainty of $\delta_{\rm CP}$ measurement by the ESSnuSB experiment. Since the primary objective of ESSnuSB is to perform a precision measurement of $\delta_{\rm CP}$ in addition to discovering it (if next-generation LBL experiments fall short), it is imperative to see how new physics affects this capability of ESSnuSB. In Ref. Alekou:2022emd ; ESSnuSB:2023lbg , it has been shown that the optimal baseline of 360 km allows the ESSnuSB experiment to measure $\delta_{\rm CP}$ with a 1$\sigma$ uncertainty of less than $7.5\degree$ for all possible values of $\delta_{\rm CP}$. Remarkably, the experiment achieves its best precision, $\Delta\delta_{\mathrm{CP}}=5^{\circ}$, for CP-conserving values. Such a level of accuracy is unparalleled, as it surpasses the capabilities of next-generation LBL experiments, emphasizing the transformative potential of ESSnuSB in this area of research.

In Fig. 11, we present the projected 1$\sigma$ uncertainty in the measurement of $\delta_{\rm CP}$ for two scenarios: the CP-conserving values ($\delta_{\rm CP}=0\degree$ and $180\degree$, shown in the left panel) and the maximally CP-violating values ($\delta_{\rm CP}=\pm 90\degree$, shown in the right panel), as functions of the LRF potentials, $V_{\alpha\beta}$. In both panels, solid, dashed and dashdot curves represent LRF potentials, namely, $V_{e\mu}$, $V_{e\tau}$ and $V_{\mu\tau}$, respectively. From Fig. 11 it is evident that the effects of LRF potentials on the $\delta_{\rm CP}$ precision of ESSnuSB are negligible. Even when the values of all three potentials are large enough (almost an order of magnitude larger than the ESSnuSB bounds), the effects of $V_{\alpha\beta}$ are not significant enough to meaningfully degrade the performance of ESSnuSB. Specifically, for the maximally CP-violating values ($\delta_{\rm CP}=\pm 90\degree$), illustrated in the right panel of Fig. 11, the experiment can achieve a robust precision of $\Delta\delta_{\rm CP}<7.5\degree$, as long as $V_{\alpha\beta}$ remain below $2\times 10^{-14}$ eV. For the CP-conserving values ($\delta_{\rm CP}=0\degree$ and $180\degree$), illustrated in the left panel of Fig. 11, the precision is even better, with $\Delta\delta_{\rm CP}\lesssim 7\degree$ across the entire range of LRF potentials, $V_{\alpha\beta}$.

These results highlight the resilience of ESSnuSB in maintaining high precision in $\delta_{\rm CP}$ measurements, even in the presence of LRFs, further demonstrating its capability to probe CP violation with unprecedented accuracy.

In this paper, we explored the capabilities of the ESSnuSB experiment to set bounds on the effects of LRFs in neutrino oscillations. In the presence of additional $U(1)$ gauge symmetries in the particle physics Lagrangian, a new vector mediator $Z^{\prime}$ might be responsible for new interactions between SM particles. In the case of a very light mediator, such interactions might occur at very long distances and feebly interacting particles like neutrinos could provide valuable information about them. For instance, in neutrino oscillation experiments, LRFs modify matter effects in the neutrino oscillation probabilities introducing new terms in the Hamiltonian. We considered three different $U(1)$ symmetries, namely $L_{e}-L_{\mu}$, $L_{e}-L_{\tau}$ and $L_{\mu}-L_{\tau}$. We demonstrated how the ESSnuSB setup could provide a good environment to search for LRFs. In particular, using nominal conditions (5% systematics), we observed that ESSnuSB could be able to set 90% C.L. limits on $V_{e\mu}<2.99\times 10^{-14}$ eV, $V_{e\tau}<2.05\times 10^{-14}$ eV and $V_{\mu\tau}<1.81\times 10^{-14}$ eV. The bounds on such parameters have been obtained by means of a standard $\chi^{2}$ analysis performed using the GloBES software. Among the upcoming next-generation LBL experiments, only DUNE is expected to outperform ESSnuSB, while T2HK will set weaker limits Singh:2023nek ; Giarnetti:2024mdt . The ESSnuSB bounds might become comparable to the DUNE ones if systematic uncertainties in both the appearance and disappearance channels are reduced for the ESSnuSB experiment. We explored the correlations between the LRF parameters and the most unknown oscillation parameters, namely $\theta_{23}$ and $\delta_{\rm CP}$. We found that the octant degeneracy of $\theta_{23}$ is broken in the presence of LRFs when $\theta_{23}^{\rm true}=42.2\degree$. We also could not observe any strong correlation between $\delta_{\rm CP}$ and the LRF potentials $V_{\alpha\beta}$.

Finally, we addressed another crucial point in the context of the ESSnuSB experiment: the robustness of its most important measurement, namely the $\delta_{\rm CP}$ determination. We observed that, even in the presence of LRFs, both the CPV sensitivity and the $\delta_{\rm CP}$ precision remain unaltered except in the case of extremely large LRF potentials ($V_{\alpha\beta}\gg 10^{-13}$).

Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them.

We acknowledge further support provided by the following research funding agencies: Centre National de la Recherche Scientifique, France; Deutsche Forschungsgemeinschaft, Germany, Projektnummer 423761110; Ministry of Science and Education of Republic of Croatia grant No. KK.01.1.1.01.0001; the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 860881-HIDDeN; the European Union NextGenerationEU, through the National Recovery and Resilience Plan of the Republic of Bulgaria, project No. BG-RRP-2.004-0008-C01; as well as support provided by the universities and laboratories to which the authors of this report are affiliated, see the author list on the first page.