Uncovering Secret Neutrino Interactions at Tau Neutrino Experiments

We investigate the potential of future tau neutrino experiments for identifying the $\nu_\tau$ appearance in probing secret neutrino interactions. The reference experiments include the DUNE far detector utilizing the atmospheric data, which is for the first time in probing the secret interactions, the Forward Liquid Argon Experiment (FLArE100) detector at the Forward Physics Facility (FPF), and emulsion detector experiments such as SND@LHC, AdvSND, FASER$\nu$2, and SND@SHiP. For concreteness, we consider a reference scenario in which the hidden interactions among the neutrinos are mediated by a single light gauge boson $Z'$ with a mass at most below the sub-GeV scale and an interaction strength $g_{\alpha \beta}$ between the active neutrinos. We confirm that these experiments have the capability to significantly enhance the current sensitivities on $g_{\alpha \beta }$ for $m_{Z'} \lesssim 500$ MeV due to the production of high energy neutrinos and excellent ability to detect tau neutrinos. Our analysis highlights the crucial role of downward-going DUNE atmospheric data in the search for secret neutrino interactions because of the rejection of backgrounds dominated in the upward-going events. Specifically, 10 years of DUNE atmospheric data can provide the best sensitivities on $g_{\alpha \beta}$ which is about two orders of magnitude improvement. In addition, the beam-based experiments such as FLArE100 and FASER$\nu$2 can improve the current constraint on $g_{e\tau}$ and $g_{\mu\tau}$ by more than an order of magnitude after the full running of the high luminosity LHC with the integrated luminosity of 3 ab$^{-1}$. For $g_{e\mu}$ and $g_{ee}$ the SHiP experiment can play the most important role in the high energy region of $E>few~100$ MeV.


I. INTRODUCTION
Among the particles in the Standard Model (SM), neutrinos are unique in the sense that they play the key roles not only in determining the weak interaction structures but also in guiding a new physics beyond the SM (BSM) due to their oscillation phenomena, not explained in the context of the SM.Possible new interactions of neutrinos other than the weak interaction, therefore, can shed light on identifying the symmetrical structure of BSM.
One area of interest is the Secret Neutrino Interaction (SNI), which involves new boson(s) mediating the interactions among the active and sterile neutrinos, or involving only the sterile neutrino sector.Secret neutrino interactions might arise in the BSM theories with the neutrino masses given from the breakdown of the global symmetries of the SM such as lepton number (L) or the difference between the baryon number and the lepton number (B − L) symmetries [1][2][3][4][5][6].Other possibilities include gauging an anomaly free global symmetry [7,8], which is not technically "secret", or introducing a new gauge symmetry completely blind to the SM particles [9].
The theoretical scenarios providing SNI have been applied to explain the neutrino oscillation anomalies [10][11][12][13][14][15].Interestingly, the SNIs have been also used to resolve various issues in cosmological and astrophysical observations.The pseudo Nambu-Goldstone boson arising after the spontaneous break down of a global lepton number symmetry (L or B − L) and the electroweak symmetry, so called Majoron, can be a dark matter (DM) candidate [16].The emission of Majorons or vector bosons can also contribute to the supernova cooling [17,18].Inclusion of SNI can make the thermal sterile neutrino DM scenario [19] compatible with the astrophysical observations [20].A new gauge boson mediating the SNI can be used to resolve the small scale problems, albeit with strong cosmological constraints [21][22][23], or alleviating the Hubble parameter tension [24][25][26].
Due to its importance, the investigation of SNI has been rigorously pursued across multiple domains including cosmological, astrophysical, and laboratory experiments.Among them, various astrophysical and cosmological observational results prefer flavor non-universal secret interactions [25,27].Moreover, the laboratory experiments have provided stronger constraints on the SNIs with ν e and ν µ than those with ν τ [28][29][30][31][32].
Motivated by these observational and experimental preferences, we explore the scenarios where the SNIs are flavor non-universal and the mediators do not interact with the charged leptons in this paper.Notably, we focus on the exciting potential of a variety of future tau neutrino experiments in directly probing the ν τ SNIs with less constraints from the laboratory experiments so far, compared to the other flavor SNIs.
For concreteness, in this paper, we adopt a light vector SNI scenario where a vector with the mass below 1 GeV couples exclusively to the SM active neutrinos, described by the term α,β g αβ Z ′ µ να γ µ ν β , which provides Non-Standard Interactions (NSIs) involving light mediators [9,33,34].A viable scenario addressing these NSIs is proposed in Ref. [9], which is briefly explained in the next section.The light Z ′ can be produced via three-body rare decays of pseudoscalar mesons, if kinematically available, which can be important decay channels due to no chiral suppression.Sensitivities of meson decay experiments to secret couplings of neutrinos to Z ′ is studied in Ref. [35].We note that while meson decay experiments are sensitive to the sum of the coupling strength squares involving charged leptons produced in the decay of charged mesons, i.e., α∈{e,µ,τ } |g eα | 2 and α∈{e,µ,τ } |g µα | 2 (by identifying the produced charged lepton), neutrino detectors can detect produced neutrinos and are sensitive to each of the couplings g eα , g µα and g τ α .In order to obtain conservative sensitivities, we further assume that both of the production and the decay of Z ′ are controlled by a single SNI parameter g αβ .In particular, we discuss the possibility of using the tau neutrino flux measurement to constrain the coupling of neutrinos with the new light gauge boson.The upcoming neutrino detectors can benefit from their capability to detect high-energy neutrinos produced from heavy mesons.Additionally, their abilities of detecting ν τ directly can provide superb sensitivities on the NSI couplings [36].
As laboratory experiments probing beam-produced neutrinos, we adopt Forward Liquid Argon Experiment (FLArE100), SND@LHC, FASERν2, and SHiP for reference.The FASER [37], FASERν [38], and SND@LHC [39] detectors are currently under operations in the tunnels located in the beam forward direction nearby ATLAS and have recently announced their first phase data [40,41], which has opened up an era of intensity frontier searches for BSM at the LHC.To succeed these experiments, the Forward Physics Facility (FPF) which aims to host the next generation experiments during the running of the High Luminosity LHC (HL-LHC) is proposed [42].
The FPF neutrino experiments are expected to detect high number of neutrino interactions at the highest energies ever achieved.Thus, their measurements are crucial to uncover neutrino interactions at energies above O(100 GeV).The proposed FLArE, a liquid argon time projection chamber (LArTPC) located at FPF as well as FASERν2 are designed to detect millions of neutrino interactions, including tau neutrinos, and to search for long-lived BSM particles or dark matter.SHiP is an intensity-frontier beam dump proposed experiment which aims to explore the domain of weakly interacting hidden light particles with masses in the MeV-GeV range.Sensitivity of the currently running FASERν to secret neutrino interaction is previously studied in Ref. [43].
In addition to the laboratory experiments, we consider Deep Underground Neutrino Experiment (DUNE) far detector to analyze the atmospheric data in probing SNI for the first time.We place special emphasis on the sensitivity of DUNE atmospheric data to SNIs by probing the appearance of downward-going ν τ and discuss how this strategy is strong, as well as providing valuable insights on the flavor structure of Z ′ .Note that DUNE will have both far detector, probing high energy atmospheric neutrinos, and near detector for high intensity lower energy neutrinos; we can expect the interplay between those in probing new interactions of neutrinos.Moreover, the two detectors' excellent event reconstruction, angular resolution, and abilities in identifying ν τ directly with track reconstruction can provide key information on SNI.The potential of DUNE near detector (ND) to constrain the new interaction is studied in Ref. [44].Similar study on ν τ appearance by a neutrino-philic mediator including those for the SNI in short-baseline laboratory experiments is recently proceeded in Ref. [45].This paper is organized as follows.In Sec.II, we explain our simplified set-up from which useful analytic formulas are obtained.In Sec.III, we summarize the experimental details of our reference tau neutrino experiments and discuss the analysis strategies.In Sec.IV, we then show our analysis results for each coupling g αβ while turning off the others, Finally, we conclude our results and leave discussion in Sec.V.

II. THEORETICAL SET-UP
In this section, we explain our simplified set-up of vector SNI with a new sub-GeV mass vector boson Z ′ .The relevant effective Lagrangian includes: where g αβ represents the coupling between the new light boson Z ′ and neutrinos of flavor α and β, which does not have to be flavor-diagonal .This interaction can lead to a new decay mode of meson to lepton, neutrino, and Z ′ , which is followed by a subsequent Z ′ decay. 1ote that this interaction can arise from gauging different combinations of baryon number and lepton flavor/number [46,47].However, the coupling of the electron to the new gauge boson Z ′ is subject to stringent constraints across a wide range of Z ′ mass and hence our models of interest should suppress the sizable couplings to the SM charged leptons.A possible scenario giving rise to this interaction can be obtained from adopting a new gauge symmetry U(1) ′ along with a SM singlet heavy fermion Ψ and a scalar particle both of which are charged under U(1) ′ [9,48].Then the active neutrinos can couple to Z ′ by mixing with Ψ when the new scalar particle is either a SM singlet or doublet; the active neutrinos of flavor ν α can be written as a linear combination of mass eigenstates ν i , (i = 1, 2, 3, 4): where ν 4 is the heaviest mass eigenstate that gives the main contribution to Ψ. Integrating out the heavy fourth state, the light active neutrinos receive a coupling of the form g βα Z ′ µ νβ γ µ ν α where g βα = g Ψ U α4 U * β4 and g Ψ being the U(1) ′ gauge couping of Ψ.Note that a kinetic mixing between U(1) ′ and U(1) Y can generically arise.Hence additional theoretical set-up should be assumed in such a way that the tree level mixing is turned off and the loop level mixing is induced by very heavy particles.Our process of interest is depicted in Fig. 1.In order to focus on the SNI from our reference scenario, we do not consider the baryonic couplings of Z ′ throughout the whole analysis.The flux of neutrinos in the lab frame coming from meson decay, M → Z ′ νℓ with subsequent decay of Z ′ → ν ν, is given by the equation: Here, L represents the distance from the source to the detector and P M (E M ) is the rate of the meson injection in the lab frame.
The spectrum of the neutrino in the lab frame from the decay of a meson with an energy of E M is expressed as dNν dEν lab , which is related to the spectrum of neutrinos in the rest frame of meson, The value of E ν | lab can be simply obtained by kinematics.If the direction of a neutrino reaching the detector coincides with that of the spatial momentum of the meson beam, for instance, we can write Note that dΩ r.M /dΩ lab takes care of focusing of the beam in the direction of the detector and is given by ( is the total spectrum of the (anti-)neutrino produced from both meson and Z ′ decay: where N 0 is the total number of the neutrinos produced from the M + (M − ) decay, which is ν α (ν α ) in the left panel of Fig. 1.For the electron decay mode M → eνZ ′ , we neglect the mass of electron and obtain the decay rate analytically, while the values for the heavier leptons, M → µνZ ′ and M → τ νZ ′ are obtained numerically.
As previously stated, the presence of the new light gauge boson Z ′ leads to an enhanced threebody decay rate of the pseudo scalar meson without chiral suppression compared to the two-body decay cases due to the longitudinal component of the new massive gauge boson.This enhancement arises from the new interaction between the new massive gauge boson and neutrinos, and the decay rate scales as g 2 αβ /m 2 Z ′ , which results from the summation over the new gauge boson polarizations, i.e., Σ i ϵ This phenomenon is analogous to the W boson emission in top quark decay (t → bW ), where the decay rate is proportional to 1/m 2 W .The enhancement in both cases arises from the polarization sum resulting from the spontaneously broken gauge symmetry.For the decay modes into e ± , the differential decay rate with polarization perpendicular to the Z ′ momentum (ϵ 1 , ϵ 2 ) and parallel to the Z ′ momentum (ϵ 3 ) can be expressed as: where f M is the meson decay constant.Observing the decay into the longitudinal mode, we can see that it is proportional to g 2 eα /m 2 Z ′ and will be enhanced for m Z ′ ≪ m M .For the case of pion which is dominantly produced in various accelerators, the total decay rate is given by: Figure 2 shows the branching ratio of the meson's three-body decay to electron, anti-neutrino, , as a function of m Z ′ for various mesons, including π, K, and D s assuming g ee = 0.1 (left) and g ee = 10 −5 (right).For large new gauge couplings of ≲ O(0.1), we can clearly see that the three-body branching ratios can be dominant over the conventional chiral suppressed two-body decays and would be easily constrained by various experiments.As one expects, the branching ratio decreases rapidly as the gauge boson mass approaches the mass of the charged meson which makes the process kinematically forbidden.The Z ′ gauge boson with masses of O(MeV − 100 MeV) subsequently decays into ν ν before reaching the detector producing signals at neutrino detectors over a wide range of g αβ .The total decay rate of Z ′ −→ ν α νβ for all the polarizations is given by: The number of neutrinos originating from the decay of Z ′ particles before reaching the detector is given by where N 0 is the number of produced Z ′ , L is the distance between the Z ′ production point and the detector, and γ = E Z ′ /m Z ′ is the boost factor.In Eq. ( 10), if ΓL/γ ≫ 1, almost all Z ′ particles decay before reaching the detector.Before closing this section, let us briefly comment on the cosmological effect of a light Z ′ .For Z ′ with masses m Z ′ ≲ O(MeV), it is inevitable to consider its contribution to the radiation energy density without Boltzmann suppression at the time of Big Bang Nucleosynthesis (BBN).In the presence of SNI, the new gauge boson can be generated via inverse decay ν + ν → Z ′ and neutrinoantineutrino annihilation ν + ν → Z ′ +Z ′ [49].The new gauge boson Z ′ can contribute to the extra radiation species ∆N eff when it is in thermal equilibrium with active neutrinos around T ∼ 1 MeV after neutrino decoupling era. 2 As a conservative limit, we adopt the constraints on ∆N eff with 90 % C.L., which is ∆N eff ≲ 1, from Ref. [49] assuming a flavor universal SNI.The combined constraints with the abundances of the primordial elements are similar in the reference.We expect the flavor non-universal and off-diagonal cases involving ν τ which is more proper to be applied in our analysis would provide weaker bounds but we do not pursue this direction here and leave more detailed study to a future work.Hence we simply apply the nominal bound on the vector boson for m Z ′ ≲ 5 MeV from Ref. [49].Note that possible baryonic interactions of Z ′ can provide extra constraints on the abundances of the primordial elements but we do not include those to focus on SNI here, as stated earlier.Other scenarios such as scalar SNI can have weaker bounds due to the smaller degrees of freedom.

III. ANALYSIS STRATEGIES A. Accelerator based Neutrino Experiments
In this section we explain the details of the reference experiments and the analysis strategies.The FLArE100, FASERν2, SND@LHC, and SHiP experiments allow us to probe the relevant parameter space for the relatively heavy Z ′ up to ≲ 1 GeV since their beam energies are high enough to produce heavy mesons such as charmed mesons.Moreover, these detectors have the potential to collect a large number of tau neutrino events and identify those, providing an opportunity to use tau neutrino flux measurements in probing possible new interactions of neutrinos.To compute the number of events, we have taken the neutrino cross section and the energy spectra of the charged mesons from Ref. [51].We further assume the perfect energy resolutions and 80% efficiencies for ν e and ν µ in the aforementioned experiments for simplicity.
The FPF is expected to host far-forward experiments such as FASERν2, a 20-ton emulsion detector; Advanced SND@LHC (AdvSND), a successor to SND@LHC; and FLArE, a proposed liquid argon time projection chamber with an active volume of 100 tons [42].These experiments have potential to detect millions of TeV-energy neutrinos.The AdvSND features a 5-ton fiducial mass that represents a substantial increase of 6.25 times compared to the SND@LHC experiment.Furthermore, it is expected to have the final integrated luminosity 20 times higher than SND@LHC (150 fb −1 ), resulting in a total 125 times larger data.We have taken the details of the aforementioned experiments from Ref. [52].
The Search for Hidden Particles (SHiP) experiment is a proposal of fixed target facility at the CERN Super Proton Synchrotron (SPS) which aims to search for light BSM particles with tiny interactions with the SM particles avoiding the experimental constraints thus far, so called hidden particles [53].The other main purpose of SHiP is directly observing ν τ and ντ events.Benefiting from high statistics it can perform active neutrino physics.Inside SHiP, a detector called SND@SHiP will be installed for the study of active neutrino cross-sections and angular distributions.This is expected to be located about 46 m behind the interaction point and detect about 12000 ν τ events within 5 years of operation, which is quite large compared to FASERν experiment detecting about 11 ν τ events by 2023.Interestingly, SHiP hosts a hadron absorber that light mesons such as pions or Kaons can interact with before decaying to neutrinos.Hence, the fraction of the charmed meson increases compared to the lighter mesons.Since the D s meson decays are the main sources of tau neutrino production at the SPS with the beam energy 400 GeV, it is possible to have a large tau neutrino flux.Moreover, SHiP will have the opportunity to observe the tau anti-neutrino for the first time and perform its cross section measurements.In this respect, SND@SHiP will also be an excellent experiment searching for the BSM particles interacting with tau neutrinos.

B. Atmospheric Neutrino Experiments
In addition to the beam produced neutrinos, atmospheric neutrinos can be used to set stringent bounds on the new couplings.Although the number of charmed mesons produced in the atmosphere is smaller than those in the aforementioned accelerator experiments, it is possible to effectively probe tau neutrinos by reducing most of the backgrounds remarkably.We adopt DUNE far detector to confirm our expectation.
All of the current neutrino oscillation experiments explore ν µ and ν e disappearance and ν µ → ν e appearance channels (plus anti neutrino channels).Both atmospheric and neutrino beam experiments have confirmed the ν µ → ν τ oscillation by the disappearance of ν µ .This is because the reconstruction and identification of ν τ events pose significant challenges due to the prompt and semi-visible decay of the τ leptons.In particular, the misidentified neutral current (NC) scattering of any flavor neutrinos can mimic the τ lepton signal and is hard to be rejected.Moreover, the energy threshold to detect the charged current (CC) scattering of ν τ off the matter producing τ lepton is as high as E ντ ≳ 3.35 GeV for the nuclear scattering process and E ντ ≳ 3.1 GeV for the electron scattering process, which are mostly beyond the reach of the current beam neutrino experiments.
On the other hand, DUNE is expected to have capabilities of identifying and reconstructing the τ lepton signals due to the characteristic of the Liquid Argon Time Projection Chamber (LArTPC) detector with an excellent position resolution [54,55].In particular, the Long-Baseline Neutrino Facility (LBNF) will be equipped with the 120 GeV Neutrinos at the Main Injector (NuMI) beam providing the center-of-mass energies well above 3 GeV to observe the ν τ CC processes and the near detector complex will host a variety of detectors which can reduce the backgrounds [56].
In LBNF, the first oscillation maxima for DUNE occurs around 2.5 GeV which is below but very close to the tau neutrino detection energy threshold.This will cause some ambiguity in the measurement of ∆m 2  31 and sin 2 θ 23 [54].By comparison, DUNE far detector covering a wide range of L/E and benefiting from large flux can provide a promising tool to search for ν τ .Hence atmospheric data can provide a clearer first oscillation maxima together with an excellent angular resolution (zenith angle resolution is ∼ 5  I.Estimated numbers of standard model neutrino events assuming a final integrated luminosity of 150 fb −1 for SND@LHC, while 3000 fb −1 for FASERν2 and FLArE100.For SHiP, we assume 2 × 10 20 POT in five years.We assume a data-taking period of 10 years for DUNE atmospheric neutrinos.
[57], making the atmospheric data advantageous in more accurate measurements of oscillation parameters.
Notice that upward-going atmospheric neutrinos travelling through a larger baseline can effectively oscillate into ν τ .On the other hand, we do not expect ν τ signals within the standard oscillation model for downward-going atmospheric neutrinos since their baselines are too short to oscillate.Therefore, the unexpected downward-going ν τ appearance will be a unique signal of nonstandard interaction without suffering from large background contamination in the atmospheric neutrino experiments such as DUNE far detector.In addition, the far detector of DUNE has much larger fiducial volume than those in the accelerator based experiments, increasing its capabilities in searching for non-standard interactions.Note that similar expectation of the unexpected ν τ appearance in the beam produced short-baseline neutrino experiments is studied in Ref. [45].In this study we will explore how atmospheric neutrino data of DUNE can provide wonderful sensitivities in probing the SNI.The details of the experiments we have used are given in table I.

IV. RESULTS
In this section, we show our analysis results by displaying the expected sensitivities of various reference experiments in the plane of g αβ − m Z ′ along with the current constraints.In all of our results, we assume a chosen g αβ is the only non-zero SNI coupling to make the analysis conservative and simple.From observing the ν τ events, the SNI couplings g ατ can be directly probed.
The sensitivities on g τ τ can be, in principle, dominantly obtained from the process of D s → τ ν τ Z ′ → τ 3ν τ .However, we expect those are very weak due to the small flux of D s and the phase space suppression for a given mass of Z ′ .On the other hand, g eτ and g µτ can be probed in the processes D s → eν e Z ′ and D s → µν µ Z ′ , respectively, again without chiral suppression compared to D s → τ ν τ Z ′ , providing more phase space.In addition, those SNIs can be probed from the lighter meson decays.We estimated that BR(D s → τ ν τ Z ′ ) for g τ τ = 0.1 and m Z ′ ∼ 10 MeV is about 10 −4 times smaller than BR(D s → eν e Z ′ ) in Fig. 2.
In order to show the effectiveness of our analysis strategies observing the ν τ events, we analyze the other SNI couplings, i.e., with ν e or ν µ but not ν τ .Note that our Z ′ from the reference models might be further constrained by its baryonic interactions but we do not include such a possibility as a conservative approach.FIG. 3. The upper bound on g eτ vs. m Z ′ at 90% C.L..The red curve corresponds to the sensitivity of DUNE with the ten years of atmospheric data assuming no background.The dotted red curve corresponds to the DUNE atmospheric data taking into account the background by NC scatterings of any flavor neutrinos estimated in Ref. [58], which is 70 events for 10 years.The cyan and purple curves show the sensitivity of FLAE100 and FASERν2 to constrain g eτ , respectively.The blue and curve corresponds to the SHiP, i.e., SND@SHiP experiment and the green one corresponds to the Advanced SND@LHC experiments.The dashed black curve shows the bound from DUNE near detector [44].The gray region shows the BBN constraint [49].The dark gray and light green regions shows the current constraint from Z boson decay and NA62, respectively [59,60].The light blue region indicates the constraint from core collapse supernova [18].
Figure 3 displays the 90% C.L. current constraints and future sensitivities on g eτ versus m Z ′ , while all the other SNI couplings are set to zero.The analysis takes into account all meson decays, including π, K, and D s to leptons and neutrinos (π, K, D s → l, ν).The green region indicates the current exclusion limit from NA62 [60], while the dark and light gray regions represent the constraints from Z boson decay and BBN, respectively [49,59].Note that the late decay of Z ′ → ν α ν β prior to the recombination epoch can possibly contribute to extra ∆N eff from the observation of Cosmic Microwave Background (CMB), which will be stronger than the currently displayed BBN bound.However, the detailed fitting of the CMB data in the presence of flavor nonuniversal and off-diagonal SNI is nontrivial and hence we leave more dedicated study to a future work without displaying those bounds throughout this work.We can also apply a bound from the observation of the power spectrum of CMB due to the late neutrino free streaming [27,32,61] but it is far weaker than that from NA62 for a simple universal couplings case, g ee = g µµ = g τ τ .We hence do not show the CMB power spectrum constraint here.The light blue region indicates the constraint from core-collapse supernovae [18].
The plot demonstrates that FLArE100 (cyan curve) and FASERν2 (purple curve) can set comparable and the most stringent constraints on g eτ among future beam experiments.Note that FLArE100 has the largest fiducial volume with the background events comparable to the much smaller detectors, SND@SHiP or AdvSND, as can be seen in Table I.Also, FASERν2 has a smaller fiducial volume but with much smaller number of backgrounds even compared to FLArE100.With the difference in the shape of the neutrino flux, the above advantages make FLArE100 and FASERν2 promising in probing the g eτ coupling.For m Z ′ ≳ few MeV, they can improve the current constraint by about one order of magnitude, and by more than one order of magnitude for m Z ′ ≲ 60 keV.Additionally, the plot shows that SHiP is more sensitive to the new coupling above a few MeV compared to SND@LHC.This is due to the fact that the fraction of the produced charm mesons to lighter mesons at SHiP experiment is higher than SND@LHC due to the presence of hadron absorber.It is important to note that our analysis focuses solely on the far detector, where we expect a similar neutrino flux as SND@LHC.However, it is worth mentioning that incorporating the near detector may potentially lead to significant improvements in the obtained results.As depicted in the plot, Advanced SND@LHC demonstrates a comparable sensitivity on the parameter g eτ when compared to the FASERν2 and FLArE100 experiments.
It is remarkable that the atmospheric data at DUNE with 10 years of running can provide most stringent sensitivities for m Z ′ ≳ 1 MeV and m Z ′ ≲ 60 keV.Note that the direction of ν τ is crucial, i.e., the downward-going events determine the sensitivities.The red solid line corresponds to the zero-background assumption while the red dotted line shows the inclusions of the background by the NC scatterings of any flavor downward-going neutrinos expected in Ref. [58], which is 70 events for 10 years.In performing our analysis, we use the Honda atmospheric neutrino flux model [62].We assume 100% efficiency for ν τ event reconstruction for simplicity, following the relevant study [55].Notice that the efficiency for reconstructing tau neutrino tracks in DUNE atmospheric data depends on various factors such as the energy, direction of the neutrino, the properties of the detector, and the reconstruction algorithms used.Nevertheless, DUNE is designed to have excellent spatial and angular resolutions.The detector consists of a large volume of liquid argon, which allows for precise tracking and energy measurements of particles produced in neutrino interactions.The detector is also complemented by a highly sophisticated software system for event reconstruction, which is continually being improved to increase the efficiency and accuracy of tau neutrino reconstruction.We assume angular resolution (Θ zen resolution) of 5 • for ν τ CC [58].More dedicated study in this direction in collaboration with the experimentalists is also possible in the future.Although limited, the DUNE near detector (ND) can also play roles in observing the SNI, which is shown as black dashed line following Ref.[44].Note that these experiments can provide the searches for SNIs complementary to cosmological (gray) and astrophysical (blue) probes.In some parameter regions such as m Z ′ ≳ O(MeV) or g eτ ≲ O(10 −7 ) while m Z ′ ≲ O(10 −2 MeV), the ground based experiments exhibit better sensitivities.
Figure 4 displays the upper bound at 90% C.L. on g µτ as a function of m Z ′ .As observed in the plot, FLArE100 (cyan curve) and FASERν2 (purple curve) can improve the current constraint on g µτ by more than one order of magnitude for M Z ′ > few MeV and for M Z ′ < 60 keV.Furthermore, SHiP can slightly enhance the current constraint on g µτ for masses larger than a few MeV, whereas SND@LHC can slightly improve the current bound for masses lower than few keV.Note that the NA62 bound is not applied here since g µτ is the coupling between ν µ and ν τ .Instead, we applied the experimental constraints from K → µν µ since the new gauge boson can be produced from ν µ and produce K → µν τ ν µ ν τ [63].The corresponding bound is expressed as the yellow shaded region.Again, we expect the DUNE atmospheric data with the 10 years of running can provide the best sensitivity.For comparison, we include the analysis for g ee and g eµ in Fig. 5 and 6, respectively.We can clearly see the sensitivities are weaker by about an order of magnitude than those in Figs. 3 and 4 due to the background contamination for the other flavors of neutrino events.We have taken the SM events as the background for ν e and ν µ .The number of ν e + νe is 1.6 × 10 4 .Also the number of ν µ + νµ is 2.4 × 10 4 .Note that the shape of the spectrum is important for g ee and g eµ and the direction of the neutrinos is non-critical.
For m Z ′ < few keV, FLArE100 can improve the current constraint on g ee , while SHiP shows better sensitivities for m Z ′ ≳ 400 MeV due to its higher sensitivity to neutrinos originating from heavy meson decays, such as D s .These results highlight the importance of studying heavy meson decays to further constrain the coupling of secret neutrino interactions in the higher Z ′ mass region.
Interestingly, the atmospheric neutrino data in the 10 years of running of the DUNE far detector can still provide excellent sensitivities better than the accelerator experiments in most of the parameter space of m Z ′ ≲ 500 MeV due to the size of the fiducial volume and the shape of the flux.As the mass of Z ′ approaches to GeV level, SHiP can be more sensitive since the flux of heavy mesons in the atmosphere decreases while the large backgrounds of ν µ or ν e are still not effectively rejected.
In addition, the SNI with g ee can induce the two-neutrino double beta decay (2νββ) which is also expected in the SM, even without lepton number violating interactions [31].However, the current bound (cyan shaded region in Fig. 5) is still weaker than the combined constraints from BBN and NA62.We expect the sensitivities on g eµ are even weaker that those on g ee because the production rate of ν µ in both atmospheric and accelerator data is higher.
As can be observed from Fig. 3  the most sensitive probe on g eτ even after including the NC background.Notice that, in our analysis for atmospheric neutrinos, we have fixed the flux value and did not include the uncertainty of the shape of the flux.The obtained sensitivities on g ee and g eµ can be modified significantly including this uncertainty.On the other hand, the expected sensitivities on g eτ and g µτ are quiet robust with respect to the flux uncertainty since the standard interactions do not produce downward-going ν τ .
We emphasize again that these two couplings are sensitive to the direction of the tau neutrino and the shape of the background.
Finally, it is fair to leave a comment that our Z ′ can induce lepton flavor violating rare decays such as µ → eγ in the two loop level.However, our naive estimation shows the contribution can be negligible for the couplings below g αβ ≲ O(10 −2 ) compared to the experimental limits so far.More exact calculation is beyond the scope of this paper and we leave the related study to a future work.

V. CONCLUSION
The upcoming beam and atmospheric tau neutrino experiments offer a promising avenue to explore the hidden interactions between neutrinos, whose identification is highly crucial in various fields, including neutrino physics, dark matter physics, grand unified theories, astrophysics, and cosmology.For concreteness we adopted a scenario with a light gauge boson Z ′ with coupling g αβ to ν α and ν β .Our analysis highlights the importance of DUNE atmospheric data in obtaining the best sensitivities on g αβ for the 1 MeV ≲ m Z ′ ≲ 500 MeV mass range as well as for m Z ′ ≲ O(keV), with the potential to improve the current constraint by up to two orders of magnitude.Notice that we have assumed angular resolution of 5 • to perform our analysis for the DUNE atmospheric data, which is a conservative choice.In particular, the downward-going ν τ events, together with the help of exact identification and reconstruction of tau leptons, can be highly efficient in proging g ατ couplings.We observed that including the NC background does not change our conclusions significantly.Additionally, FLArE100 and FASERν2 have the potential to significantly enhance the current bounds on g eτ and g µτ , while also slightly improving the constraints on g ee and g eµ .Notably, in the case of g ee and g eµ , above a few hundred MeV, SHiP is more sensitive in probing the couplings due to the larger number of produced D s mesons compared to the atmospheric case and the background contamination by conventional ν µ and ν e in the atmosphere.
It is worth noting that in our analysis of DUNE atmospheric neutrinos we have assumed a fixed flux shape.Inclusion of flux shape uncertainty could significantly modify the obtained sensitivities on g ee and g eµ .On the other hand, the sensitivities on g eτ and g µτ will be highly reliable on the direction of the tau neutrino and the shape of the background instead, so our results are quite robust on the flux uncertainties.Our analysis results here can guide future experimental searches for new physics beyond the Standard Model.It is also important to note that the sensitivities obtained in this study are based on the reference scenario with a sub-GeV level new gauge boson.Other theoretically well-motivated models that predict different types of secret interactions may result in different sensitivities as well as the astrophysical and cosmological constraints, which is worth to be studied in a future work.In conclusion, the currently on-going and future tau neutrino experiments such as DUNE, FLArE100, FASERν2, SND@LHC, and SND@SHiP, have great potential to search for a hidden interaction between neutrinos mediated by a new light sub-GeV gauge boson.In this regard, we emphasize on the importance of increasing the efficiency and accuracy of tau neutrino reconstruction in searching for BSM interacting with tau neutrinos and encourage the experimental colleagues in this direction.

FIG. 1 .
FIG.1.Three body meson decay M → lνZ ′ and subsequent Z ′ decays into a pair of a neutrino and an anti-neutrino.

FIG. 2 .
FIG. 2. The branching ratio of the meson three-body decay to electron, anti-neutrino and new gauge bosonΓ(M →eνZ ′ ) Γ SM total

TABLE
• for ν τ CC and ∼ 7 • for NC) and energy resolution