Longer-Lived Mediators from Charged Mesons and Photons at Neutrino Experiments

Since many of the dark-sector particles interact with Standard Model (SM) particles in multiple ways, they can appear in experimental facilities where SM particles appear in abundance. In this study, we explore a particular class of longer-lived mediators that are produced from photons, charged mesons, neutral mesons, and $e^\pm$ that arise in proton-beam fixed-target-type neutrino experiments. This class of mediators encompasses light scalars that appear in theories like extended Higgs sectors, muon(electro)philic scalars, etc. We evaluate the sensitivities of these mediators at beam-based neutrino experiments such as the finished ArgoNeuT, ongoing MicroBooNE, SBND, ICARUS, and the upcoming DUNE experiment. We realize that scalars are more enhanced while produced from three-body decay of charged mesons, especially if they are muonphilic in nature. For scenarios that contain muonphilic scalars, these experiments can probe unexplored regions of parameter space that can explain the current discrepancy in the anomalous magnetic moment of muons. The sensitivity of electrophilic scalars at the DUNE Near Detector can explore new regions. We also show that Bethe-Heitler scattering processes can be used to probe flavor-specific lepton final states even for the mediator masses below twice the lepton mass.


I. INTRODUCTION
There are compelling reasons for the existence of a particle sector (often called dark sector or hidden sector) beyond the Standard Model (SM) of particle physics, for example, dark matter, non-zero neutrino masses, massflavor hierarchy puzzle, etc.An attractive scenario consists of a new particle sector that is very weakly or feebly connected to the SM sector via portal particles that are often called mediators [1].A myriad of models with mediators has been built along this line to address these issues and similarly, many experiments are being developed towards unraveling these mysteries.These efforts are also motivated by anomalies in the experimental results such as the LSND excess [2], the MiniBooNE anomaly [3][4][5], and the discrepancy in the anomalous magnetic moments of the muon [6,7] and electron [8].A subset of these experiments are fixed-target-type experiments involving high-intensity protons on target (POT) and they are widely adopted at neutrino facilities with beam energy being O(1) to a few hundred GeV.While neutrino facilities serve as neutrino factories, copious amounts of charged mesons, (secondary) photons, electrons, positrons, and neutral mesons are also produced.Therefore, given the beam energy and intensity of neutrino facilities, they can test MeV-to-sub-GeV-scale new physics particles interacting with those SM particles that can be found at these facilities.
While the landscape of dark-sector models is vast, we will focus on a particular class of models that consti-tutes scalar mediators that couple either to all SM matter particles or a subset of flavors.We find that the flux of the above mediators produced from charged mesons inside a proton target can be quite enhanced.Scalars could abundantly appear through the three-body decay of charged mesons such as charged pions and kaons; their corresponding two-body decay is helicity-suppressed but adding another particle in the final state would evade the suppression and enhance the branching fraction of three-body decay modes [9].We also notice an enhancement when they are sourced from photons.For example, scalars can be copiously produced as a consequence of Primakoff scattering [10][11][12] of photons as they interact with atoms present in the target.This process is coherently enhanced by a factor of Z 2 from the nuclear form factor [13].Similar to scalars, flavor-specific massive vector mediators can be produced from the above sources as well.One such example we considered in this paper is a muonphilic gauge boson that appears in the context of a U (1) T 3R model [14][15][16][17].We find that these can also appear in good abundance from charged mesons, neutral mesons, photons, electrons, and bremsstrahlung.
We investigate the detection prospects of the signals induced by the aforementioned mediators at experiments along the NuMI [18] and BNB [19] beamlines at Fermi National Laboratory: ArgoNeuT [20], ICARUS [21], Mi-croBooNE [22], and SBND [23].We also probe them at the upcoming DUNE Near Detector [24] (DUNE ND), which is placed along the LBNF beam line [25].These detectors feature different baselines and angular distances with respect to their respective beam axis.The magnetic horns -which are designed to focus or deflect charged particles and, in turn, their corresponding neutrino decay products -affect the mediator production via charged mesons.As a result of the magnetic horn effect along with the position of detectors, we expect to take advantage of the multiple experiments and probe regions of parameter space in a complementary manner.
Once a mediator is produced inside the proton target of these beam facilities, it should be safely delivered to a detector of interest.Due to the feebly-interacting nature of the above-mentioned mediators, they would live longer rather than decay immediately.Once they survive, some fraction of mediators can leave detectable signatures within the detector fiducial volume.These signatures include electron-positron pairs, muon-antimuon pairs, photon pairs, electron-photon pairs, and single photons from scattering and decay processes.We again expect to benefit from different baselines and detector sizes in the search for these long-lived mediators as they provide complementarity.
We emphasize that while one can look for electronpositron pairs and muon-antimuon pairs from decay processes, the same final states can arise through the splitting process, also known as the Bethe-Heitler scattering process [26].Owing to its energy-dependent nature, these final states can also be produced from a mediator with a keV-range mass, which is not kinematically allowed for decay processes.For example, a 10 keV (muonphilic) scalar with a total energy greater than 210 MeV can split into a muon-antimuon pair through the Bethe-Heitler scattering process.From the above example, we see that the appearance of lepton-antilepton final states for all possible masses helps us to probe flavor-specific models directly.With all the above-mentioned detection channels, we investigate the parameter space of the Higgs Portal Scalar and (g − 2)-motivated parameter space, utilizing muon/electrophilic scalar mediators which can be efficiently probed by experiments operating at the sub-GeV scale.
In Sec.II, we discuss essential features of the example models that we explore.Section III is reserved for a brief overview of the benchmark short baseline experiments for which we study sensitivity reaches.We then explain how the aforementioned mediators are produced at generic proton-on-target experiments in Sec.IV and elaborate on the signals that the mediators produce at the detectors in Sec.V.In Sec.VI, we explain our analysis methodology and report our main results including sensitivity plots.In Sec.VII, we explore the above study for vector mediator scenarios using the U (1) T 3R model as an example.We finally summarize and conclude our study in Sec.VIII.

II. MODELS
We apply the above production and detection mechanisms of scalars in the context of three benchmark spin-0 mediator models. 1 Although we focus on scalars in the paper, one can also look at spin-1 gauge boson mediator models.We will briefly explore these aspects in Sec.VII.

A. Higgs Portal Scalars (HPS)
This model contains a singlet dark scalar S with mass m S that interacts with the SM SU (2) scalar doublet H via a portal interaction [27]: Where A and B are free parameters.Under SU (2) L symmetry breaking, the neutral Higgs decomposes into a sum v+h.Therefore, the interaction in Eq. ( 1) induces a mass mixing between the dark scalar and the SM Higgs in two ways: (a) if A ̸ = 0, then the mass mixing is naturally induced regardless of whether the dark scalar acquires a zero or a non-zero dark scalar vacuum expectation value (vev), and (b) if A = 0, the dark scalar can acquire a non-zero vev by an appropriate choice of potential and thus induce mass mixing.After diagonalizing the masslike terms, we see that the scalar S mixes with the SM Higgs h via a small mixing angle, that is, Therefore, the dark scalar can interact with the SM particles that acquire mass via the Higgs scalar: where f runs over all quark and charged lepton flavors.
This model is of particular interest in various contexts.Examples include scalar-to-pion decays [28], Mi-croBooNE searches for the HPS to explain the KOTO excess [29,30], and a search for the HPS-induced signatures at ICARUS [31].In addition to the above-shown interactions, HPS can also couple to two photons via a fermion loop and thus widen the phenomenology [32].

B. Muonphilic scalars
These scalars (henceforth denoted by ϕ µ ) can appear in effective field theories containing singlet scalars that have minimal flavor violation [33] or in other extensions to the SM [34] with additional doublets/singlets, which contain Yukawa couplings unique to each flavor [35].The Lagrangian has the following form.
Muonphilic scalars contribute to the anomalous magnetic moment of the muon (a µ ) of the muon at the one-loop level.Therefore, their phenomenology is useful to explain the 4.2σ discrepancy between the experimental and theoretical values of a µ [6,7]: Similar phenomenology has been explored via muoncoupled axion models which bring in constraints from SN1987a data [36,37].These scalars can also couple to two photons via a muon loop and the effect of this has been studied in the context of axions [38,39].The lack of photon events at the E137 SLAC experiment imposes stringent constraints on this model as well [40].

C. Electrophilic scalars
On a similar line of thought, there are models with scalars (henceforth denoted by ϕ e ) that solely couple to electrons.One such example is an effective field theory where all heavy fermions and bosons are integrated out such that we end up with scalars that exclusively couple to electrons via a Yukawa coupling.Such a model has bounds from stellar cooling [41,42], SN1987a [43], NA64 [44], Orsay [45], E141 [46] and E774 [47] that look at electron-positron as well as electron-photon final states.The relevant Lagrangian is given by Like ϕ µ searches, we can study electrophilic scalars to explain the discrepancy in the electron anomalous magnetic moment a e , which, based on a recent measurement with 87 Rb [8], is ∆a e = a exp e − a th e = (4.8± 3.0) × 10 −13 .(7)

III. BENCHMARK EXPERIMENTS
We explore the sensitivity of these models in several neutrino experiments as mentioned earlier.We tabulate key specifications of the experiments in Table .I. For the ease of the simulation, we simplify the detector geometry to cylinders.We hence specify the dimensions for the full cube-shaped detector as well as the cylinder-shaped one that is used in our simulations.Nevertheless, we expect that our main conclusions are nearly unaffected by these changes as the modified volumes fit in the original ones.The aforementioned mediators that reach Ar-goNeuT, ICARUS, and MicroBooNE are sourced from the 120 GeV NuMI beam, those at SBND are from the 8 GeV BNB beam, and finally, those at DUNE ND are produced by the 120 GeV LBNF beam.MicroBooNE and ICARUS receive a considerable number of mediators from the BNB beam on account of being placed on its axis.In this study, however, we will consider the contributions from the NuMI beam only.
The magnetic horn system present near the target plays an integral role.They operate in either neutrino mode (focusing positive mesons) or antineutrino mode (focusing negative mesons).We remark that DUNE ND, ArgoNeuT, and SBND are located on the beamline,2 whereas ICARUS and MicroBooNE are placed off-axis.Since the magnetic horns focus charged particles along the beam axis, most of the high-energy mediators (originating from high-energy charged mesons) are boosted along the beam direction, whereas softer mediators are less focused and diverge away from the axis.Hence, onaxis detectors are more sensitive to high-energy mediators as compared to those that are off-axis.Similarly, (secondary) high-energy photons are directed more in the forward direction and softer photons are directed more away from the beam axis.Therefore, high-energy mediators are likely to travel along the beam axis.These expectations are depicted in Fig. 1 which contains the energy spectra of muonphilic scalars at DUNE ND (see Fig. 1a), which is one of the on-axis detectors, and ICARUS (see Fig. 1b), which is off-axis.We see that the energies of the scalars that reach ICARUS are much lower than those at DUNE ND.

IV. PRODUCTION OF MEDIATORS
In this section, we explain how long-lived scalars can be produced from charged mesons, photons, electrons, and positrons.

Detectors
Beam, Energy Distance Angle off-axis Dimensions (l Table I: List of experiments at the NuMI, BNB, and LBNF baselines and key specifications of the detectors.The dimensions in the second last column are the adopted ones for simplified simulation purposes.The numbers quoted in the last column are the POT that we use in our study.

A. Charged meson decays
Charged pions and kaons dominantly decay into a charged lepton ℓ (= µ or e) and its neutrino counterpart through an off-shell intermediate W boson, for example, K + /π + → µ + ν µ and K + /π + → e + ν e .However, the above two-body decay processes are suppressed due to the required helicity of final state particles, which constrains the allowed phase space.This enables us to explore the production of long-lived mediators as the third decay product of charged mesons.Unlike the corresponding two-body decay, this three-body decay would not be limited by the helicity suppression [9,48,49].However, the branching fraction for a choice of coupling must not exceed the upper limit on three-body decay branching fractions of charged kaons [50] and charged pions [51].We use the three-body decay of kaons as our example in this section since the kinematically allowed phase space and mass range are larger than those of pions, although the same argument can be applied to pions.
HPS can emanate from the charged lepton leg as in Fig. 2a as well as from the W boson leg [52] as in Fig. 2b.Since they couple to leptons with a strength of m ℓ /v and 2m 2 W /v with the W boson, the latter contribution dominates the relevant decay matrix element.The HPS can couple to charged kaons, whose strength can be calculated from chiral perturbation theory, but since this term is subdominant in the relevant decay matrix element, we do not include their contribution to the decay width.Since we take the neutrinos to be massless leptons, we omit their contribution.HPS from K + → µ + ν µ S are kinematically restricted to the maximum mass reach m K − m µ = 388 MeV, whereas those from K + → e + ν e S can be as heavy as m K − m e = 492 MeV.
HPS can also be produced via a kaon two-body decay, i.e., K + → π + S (Fig. 2c) where the scalar couples to an intermediate top quark [31,53,54].This strong coupling to the top quark makes the branching ratio of the above two-body decay process dominate over all the three-body decays, but the HPS mass is limited to m K − m π = 354 MeV.HPS can also be produced via B + → K + S, which have been searched at LHCb [55], but the flux of B mesons is not large enough at the aforementioned neutrino facilities.Hence HPS sourced from B + do not produce a sizable signal flux here.
While looking at scalars with flavor-specific couplings such as muonphilic (electrophilic) scalars, the only possible diagrams are those where scalars emerge from muon (electron) legs (Fig. 2a).Therefore, they can appear from kaons via the process where the amplitude depends on Yukawa coupling squared y 2 22 (y 2 11 ).

B. Photons
Scalars couple to two photons at the one-loop level.This coupling can be written as [56,57], Here, the coupling strength g ϕγ is written in terms of the non-divergent one-loop factor with mass dimension −13 : where α em = 1/137 is the electromagnetic fine structure constant, p 1 and p 2 denote the momenta of the two photons, and the function I(β) carries information about the non-divergent fermion loop.It is generally expressed as where f (β) is defined as The coupling that appears in Eq. ( 9) is The dot product p 1 .p 2 equals m 2 ϕ /2 if the scalar decays into two photons, and (m 2 ϕ −q 2 )/2 if one of the photons is an off-shell propagator that appears in the Primakoff scattering Feynman diagram (Fig. 3a).From Fig. 4, we see that the loop factor is maximized when β is lying between 1 and 3, and it drops as β → 0, ∞.For a scalar that appears in the HPS model, all massive fermions contribute to the loop, but the dominant contribution is from those fermions with masses comparable to that of the scalar, as can be seen in the argument of I in Eq. (10).However, for those that appear in the muonphilic (electrophilic) scalar model, the only contribution to the fermion loop is from muons (electrons), which are proportional to y 22 (y 11 ).Through the above one-loop coupling, scalars can be produced from photons at the target via the Primakoff process, which is enhanced by a factor of Z 2 from the  nuclear form factor as mentioned earlier.These scalars are also highly forward-directed, i.e., in the same direction as the incoming photon.Despite this enhancement, HPS produced from kaon two-body decays is more than those produced via Primakoff scattering due to the presence of W boson and top quark couplings in the former scenario.Since these couplings do not exist in the case of muonphilic and electrophilic scalar models, the Primakoff production here is not as suppressed as it is in the case of HPS.We observe that this contribution exceeds that from kaon decays for electrophilic scalars with masses close to 1 MeV (twice the electron mass).
Electrophilic scalars can also appear when photons interact with electrons via Compton-like scattering (Fig. 3b).Although the enhancement factor here is only proportional to Z, unlike the Primakoff enhancement proportional to Z 2 , this process dominates over the scalar Primakoff process as it occurs at the tree level.For both the Primakoff and Compton-like scatterings, the mini- mum energy required to produce a scalar is where m T is electron mass m e for Compton-like scattering and nucleus mass m N for Primakoff scattering.

C. Electrons and positrons
When fast-moving positrons interact with electrons at the target, they can annihilate to a photon and electrophilic scalar, e + e − → γϕ e , as shown in Fig. 5a.However, if the energy of the incoming positron is resonated with a particular scalar mass, the electrophilic scalar can be produced directly, e + e − → ϕ e (Fig. 5b).The crosssection of this process is given by where E + is the energy of the incident positron.
In order to produce scalars via resonance processes, the center of mass energy of e + e − system must exactly match with the mass of the scalar (modulo its decay width) as suggested by the delta function.Since the center of mass energy √ s is ∼ 10 MeV for the NuMI and LBNF, where the peak energy of positrons is ∼ 100 MeV, scalars that are as heavy as 10 MeV could appear through this process.Similarly, at BNB where the peak energy of positrons is 10 MeV, a scalar of mass 3 MeV is preferred for resonant production.However, we find that this process is still subdominant at these resonant masses as compared to other processes.Therefore, we ignore the contributions from resonance while simulating scalars.

D. Simulation methods
We use the simulated fluxes of source particles, i.e., mesons, photons, electrons, and positrons, using the GEANT4 code package [58].The flux generated is based on a 120 GeV proton beam that impinges on a 150 cm graphite target.The distribution and normalization of charged-meson fluxes as they pass by the magnet are adjusted according to magnet specifications needed for these experiments [31,59].We also consider the secondary production of photons and electrons.To simulate mediators, we first calculate the probability of producing them in the center-of-mass (C.O.M) frame of the production process.Using this probability distribution, we simulate mediators in the C.O.M frame and then boost them to the laboratory frame.The probability functions in the C.O.M frame are explained below: 1.If a mediator is produced as a product of a twobody decay, such as K + → π + S, the branching ratio and the energy of the mediator are fixed for a given mass of the mediator.
2. If produced from a three-body decay such as K + → e + ν e ϕ e , we calculate the differential branching ratio as a function of energy in the rest-frame of the decaying particle.We generate random energy between the minimum and maximum energies in this frame which are weighted by the flux times the Monte Carlo volume.
3. If the mediator is from a 2-to-2 scattering process, for example, Compton-like scattering, we calculate (1/σ tot )dσ 2→2 /dt in the C.O.M frame of the process.In the above example, σ tot is the total scattering cross section of a photon of a given energy and σ 2→2 is the cross section of the Compton-like scattering.Also, t is one of the Lorentz-invariant Mandelstam variables.In this choice of frame, the energy and momenta of all the final state particles are fixed.
The angular distributions of the decay processes are approximately uniform (modulo the internal propagator effect).Hence, we randomly choose the angles in the restframe of the decaying particle, and then boost the fourmomentum along the direction of the decaying particle to the laboratory frame.For the scattering processes, however, the angular distribution is not necessarily uniform, but a function of the Mandelstam t.Thus we simulate t's and the appropriate Monte Carlo weights in the C.O.M frame, and then boost it back to the laboratory frame.2) polar and azimuthal angles, and (3) production point of the mediators in the laboratory frame for every possible mass.For each recorded event, we check whether the direction of the mediator is within the acceptance cone that the detector subtends at the point where it is produced.If so, we accept the event for the detector of interest and if not, we reject them.
Figure 6 depicts various production mechanisms that were discussed above in the context of the three models.We notice the dominance of each channel for different mass ranges for all three models.For HPS with masses less than 354 MeV, two-body decays of the kaons dominate over all the other sources by four orders of magnitude.Beyond this mass, however, the three-body decays of kaons and the Primakoff process become relevant in order.For muonphilic scalars, the three-body decays dominate for m ϕ < 388 MeV.Finally, in the case of the electrophilic scalars, the contribution from Compton-like scattering and the annihilation process are the most dominant.
While these plots help understand the probability of a scalar produced through a particular decay/scattering mechanism, we emphasize that the number of source particles and the energy spectra of both source particles and scalars are crucial when estimating the sensitivity.Hence, every production mechanism can be relevant and should be carefully analyzed.For example, electrophilic scalars from Compton-like scattering are more in number than those from e + e − annihilation at the aforementioned facilities since there are more photons than positrons.However, the scalars from the above two processes are softer in energy as compared to those that are produced via kaon decays.Therefore, depending on the analysis strategies (e.g., energy cut), one production channel could stand out more than the others and vice versa.

V. DETECTION OF MEDIATORS
In this section, we discuss methods to detect the abovementioned mediators at our benchmark detectors all of which adopt the liquid argon time projection chamber (LArTPC) technology.After collecting the mediators within the solid angle of the detector, they can be de-tected if their lifetime is long enough to survive up to the detector, without decaying into (in)visible particles before they reach the detector.To calculate the lifetime of the mediator for a given mass and momentum, we calculate the total decay width, which is the sum of the individual decay widths of all allowed decay channels.

A. Decay widths of mediators
Scalars can decay into a lepton and an anti-lepton if the scalar mass is greater than twice the mass of the lepton.The decay width of a scalar with mass m ϕ that couples to a lepton ℓ with strength y ℓℓ is where y ℓℓ is given in Eq. ( 12). 4 Scalars of any mass can decay into two photons with a decay width Here the coupling g ϕγ is given by Eq. ( 9).The −1 mass dimension in Eq. ( 9) is manifested in the 1/v proportionality for HPS, 1/m µ for muonphilic scalars and 1/m e for electrophilic scalars.Therefore, the inverse-mass scale of the photon mixing is the smallest for HPS, followed by muonphilic scalars and then electrophilic scalars.Additionally, HPS could also decay into two pions as well, whose decay width can be calculated from chiral perturbation theory [60]: We will now summarize all the possible decay channels in the context of the three scalar models.HPS decay into di-photons if their mass is less than 1 MeV.If they are heavier than 1 MeV, the electron-positron decay channel opens up.Above 210 MeV, the muon-antimuon decay channel dominates over the electrons, and above 276 MeV, the π + π − decay channel adds up.If we look at the two flavor-specific scalars, muonphilic (electrophilic) scalars up to 210 MeV (1 MeV), prominently decay into two photons.However, if they are heavier than 210 MeV (1 MeV), the muon-antimuon (electron-positron) decays take over.

B. Detection channels
For a given decay width of a mediator, the probability that it survives until it reaches the detector is where D is the distance between the production point of the mediator and the front end of the detector and λ L is the laboratory-frame mean decay length.λ L can be related to the lifetime in the laboratory-frame (t L ) and the total decay width (Γ 0 ) by the following equality.
where m X and p X are the mass and the momentum of the mediator respectively.After reaching the detector, it can be detected if they either decay into visible particles or scatter with an argon nucleus to give rise to visible particles.
For the decays, we require that it decays within the fiducial volume of the detector.Therefore, the probability of detecting a mediator inside a detector of length ∆ is given by Surviving mediators can give rise to a signal by scattering off electrons, nucleons, and/or nuclei at the LArTPC detector.The various possible scattering processes of scalar mediators are: 1. Scalar inverse Primakoff: Scalars can produce a single photon signal through inverse Primakoff scattering by exchanging a photon with the nucleus.This is also a coherent process that is enhanced by Z 2 from the nuclear form factor (Fig. 7a).
2. Bethe-Heitler process/splitting: This energydependent 2-to-3 scattering process gives rise to two leptons by scattering off of a nucleus.This is enhanced by the form factor and it can give rise to a lepton-antilepton signal even for mediators with masses less than twice the mass of the lepton (Fig. 7b).However, the mediator (X) requires a minimum threshold energy for this process to occur: Therefore, this scattering process in energydependent.
3. Inverse Compton-like scattering: This occurs when mediators scatter off an electron to produce a photon-electron signal at the detector.This channel is possible for HPS as well as electrophilic scalars (Fig. 7c).
The Bethe-Heitler process and the inverse Compton-like scattering occur for vector mediators too.The inverse Compton-like channel can appear as a tree-level process (if the vector mediator is electrophilic) or as a one-loop process (if not electrophilic) by mixing with the SM photon.By calculating the cross-section of the above scattering processes, we can arrive at the probability of scattering P scat within the detector length ∆.This is given by where n T is the number of target electrons/nucleons/ nuclei per unit volume and σ ϕ is the scattering cross-section of the process of interest.To clearly distinguish signal events from backgrounds at the detectors, various cuts are applied to the final-state particles.An example is the kinetic energy threshold for particle detection.According to MicroBooNE studies, we see that typical LArTPC detectors have kinetic energy thresholds around 20 MeV for leptons [61] and 30 MeV for photons [62].After discussing our results under no background assumptions, we detail the impact of the kinetic energy threshold in the context of electrophilic scalars at DUNE ND in Sec.VI B.

VI. RESULTS
In this section, we report the main results of our study.Figure 8 shows the number of events for different couplings as a function of mass.For the Higgs Portal Scalars, we depict the number of events at DUNE ND, whereas, for muonphilic and electrophilic scalars, we show the number of events at ICARUS and MicroBooNE as well.Since DUNE ND is placed on the beamline and is much closer to the target as compared to ICARUS and MicroBooNE to the NuMI beamline and its cumulative beam intensity is larger than that of ICARUS and Micro-BooNE, we see that, for a choice of mass and coupling, more events can be observed at DUNE ND than at the other two detectors.ICARUS detects more events than MicroBooNE as it is not only larger but also less away from the beamline (i.e., less off beam-axis) as compared to MicroBooNE.The flat contours depict final states that appear from the scattering of the scalars in the detector whereas the linear and curved contours represent the scalar decays to photon-photon and lepton-antilepton pairs, respectively.We also notice that the number of events for larger masses [m ϕ > 180(1) MeV for muonphilic (electrophilic) scalars] is less for larger couplings due to the shorter lifetimes.We now present the sensitivity reaches for our benchmark scalar mediator models delineated in Sec.II, and we also discuss the dependence of our findings on the background assumptions.

A. Sensitivity estimates
The sensitivities discussed in this section are obtained under zero background assumptions.Therefore, the contour for each experiment has been plotted for 3 events, which is the upper bound of the 95% confidence level (CL) interval for zero backgrounds.Any parameter that is enclosed inside the contours for each experiment yields more than 3 events.This is an all-inclusive sensitivity plot where the sensitivities include all possible events from decays and scatters of the mediators.Fig. 9 has the three sensitivity plots for the three models.The nature of the signals induced is explained in the caption below.
1. Higgs Portal Scalars. Figure 9a shows the sensitivity plot of HPS.The constraints are from recent studies such as Refs.[75][76][77].We included the limits for DUNE ND while the lines for other experiments have been taken from Ref. [31].To verify our results, we reproduced the ICARUS sensitivity given in Ref. [31].We first note that our predictions for DUNE ND differ from those in [77,78].These are rooted in differences in the number of events that define the sensitivity, and also energy/angular separation thresholds.The contribution to the scalars from two-body decays has been considered in Ref. [78].Here, however, we plot the contributions of scalars produced via two-body decays (red line) and three-body decays (blue line) separately to compare.We observe that the majority of scalars are produced from the two-body decay process, K + → π + S.This is due to the top quark coupling.However, only those scalars lighter than 354 MeV can be produced via this process.Although HPS with masses greater than 354 MeV can be produced via three-body decays, K + → e + ν e S and K + → µ + ν µ S, the flux of these scalars is suppressed.This is not only because it is a three-body decay, but also because the couplings are weaker than the previously mentioned top quark coupling.
We find that the sensitivity at DUNE ND is more enhanced in comparison to ICARUS, especially for larger couplings.For the masses and couplings of our interest, we see the number of scalars produced from Primakoff processes is subdominant.They are prominent only for θ values greater than 10 −3 , which is constrained by LHCb.For the masses and couplings of our interest, all signals produced by these scalars are from decay processes.We do not see signals from scattering processes because they are subdominant for weak couplings (θ < 10 −2 ).
2. Muonphilic scalars.Figure 9b depicts the sensitivity plot of the muonphilic scalar model.Since these scalars do not couple to quarks, W gauge bosons, or first-and third-generation leptons, their production modes are limited to kaon three-body decays K + → µ + ν µ ϕ µ by coupling to the muon, and Primakoff process γN → ϕ µ N .In the region represented by the current g − 2 discrepancy (represented by the pink band), the detected signals in this region are mostly di-photons.We also see a dip in the sensitivity plot at 210 MeV where muon-antimuon decays start to appear, thus reducing the scalar lifetimes.Stringent constraints appear from the 20 GeV electron beam experiment E137 [38].We converted the limits on the scalarphoton coupling in [38] into limits in the y 22 space and found that our results match with those in [34].
Since forward detectors such as DUNE ND, SBND, and ArgoNeuT (existing data) are more sensitive to high-energy mediators (which have longer lifetimes), the ceiling of the sensitivity curve for the above forward detectors is higher than off-axis detectors.Thus, sensitivities of forward detectors extend beyond E137 bounds, and into the g − 2 band.We also notice that big detectors are exposed to a great intensity of muonphilic scalars.Hence, they are sensitive to a large region of parameter space which is allowed by the E137.Since ArgoNeuT is smaller in size, it is challenging for them to probe couplings smaller than those excluded by E137.
Some regions of this allowed parameter space that are explored by the neutrino experiments are ruled out by SN1987a data.However, the astrophysical bounds can be avoided in light mediators models by the chameleon effect [79,80], where the masses of the mediators can depend on the background matter density.It must be noted that the environment in the core of a supernova (or a star) is extremely dense.In such a highly dense environment, the effective mass of the particle can be larger than it is on Earth causing an expansion of the allowed parameter space.In any case, it is important to probe the parameter space using laboratory experiments.
Unlike HPS, scattering channels are relevant to this model.Amongst the two scattering channels in this model, inverse Primakoff and Bethe-Heitler splitting, the latter dominates over the former despite the phase space suppression.This is because the Bethe-Heitler process occurs at the tree level.Additionally, two Feynman diagrams contribute to the matrix element of the Bethe-Heitler process.
3. Electrophilic scalars.The sensitivity of this model is depicted in Fig. 9c.The coupling of these scalars to electrons opens up many other production channels such as Compton-like scattering and associated production.While the most energetic scalars appear from K + decays, Compton-like scattering of photons with target electrons contributed to the most number of scalars.Since the sensitivity is majorly driven by Compton-like scattering, -in which the cross-section depends on the energy flux of the photons, ultimately resulting in σ Compton ∝ y 2 22 /m 2 ϕ -we see that the floor of the sensitivity plots derived from our studies varies as compared to that in the E137 bounds where the scalars are sourced by the decay of charged mesons (i.e., their branching ratio to a scalar is proportional to y 2 22 for masses lighter than the decaying meson).This results in different shapes of the floor which depicts the difference in the form of the branching ratios.Primakoff scattering majorly contributes to scalars with mass 1 MeV.The availability of multiple sources of electrophilic scalars comes at the cost of many constraints, with E137 [69] being the most stringent one.However, we find that the DUNE ND is sensitive to parameters that are still unconstrained by HB Stars [41] and E137.
These scalars are detected via electron-positron decay pairs for masses greater than 1 MeV, diphotons decays for masses between (0.01-1) MeV, and through scattering channels for masses below 0.01 MeV.There are three possible scattering channels for electrophilic scalars, inverse Compton-like, splitting contribution, and inverse Primakoff.Out of these three, the splitting contribution is the most dominant scattering channel for all detectors because the minimum threshold required for electronpositron splitting is very low, much lower than the scale of NuMI and BNB.
The splitting process, as seen in Fig. 7b, is unique because it results in lepton-antilepton final states even if the mediator mass is less than twice the lepton mass, which would not be kinematically allowed to decay.They become subdominant for higher mediator masses (masses close to twice the lepton mass) like all other scattering channels.However, they give us a unique way of identifying photon-less, purely leptonic final states.

B. Sensitivities with detection constraints
In our analysis and resulting sensitivities reported in the previous section, we envisioned the situation where associated backgrounds are sufficiently suppressed.While careful background estimates would lead to more precise sensitivity estimates, they certainly depend on signal channels and detector capabilities such as energy threshold, energy/angular resolutions, and particle identification.In particular, since the LArTPC technology is being developed and matured, higher-capability detectors would allow for rejecting more backgrounds while retaining signals.Nevertheless, in this section, we investigate how our sensitivity results are affected by the background assumption, especially in the context of DUNE ND.
The backgrounds for the HPS model at ICARUS and SBND have been investigated in Ref. [31]. Figure 10a roughly demonstrates the effect of backgrounds at DUNE ND by plotting the sensitivity contours for 100 events along with 3 events.We plot the 100 events line as it roughly corresponds to the worst-case reduction in sensitivity seen in Fig. 15 of Ref. [31].Though the parameter space coverage does not reduce drastically with increasing the number of signal events, improved background analysis is expected to minimize that reduction in the fu-(a) Higgs Portal Scalars.The dotted lines are from two-body decays at other experiments studied in [31].The 95% CL sensitivity lines for mS < 0.21 GeV correspond to the decay of scalars to e + , e − only.Beyond 0.21 GeV, the hanging lobe includes µ + , µ − signals and π + , π − as well for mS > 0.278 GeV.In this plot, we include 90% CL bounds imposed by CHARM [63], LSND [64], PS191 [65], NA62 [66], E949 [67], and 95% CL bounds from µBooNE [29], and LHCb [68].
(b) Muonphilic scalar model.The pink band shows the preferred regions of the current muon g − 2µ anomaly.Scalars lighter than 0.21 GeV dominantly scatter via spitting into µ + µ − , and decay into two photons, the latter dominating for m ϕ > 0.0001 GeV.For scalars greater than 0.21 GeV they decay to µ + µ − , which is manifested by the sharp drop in sensitivity.The 90% CL bounds from E137 [38], 95% CL bounds from B factories, g − 2 excluded regions, and constraints from SN cooling are depicted in the shaded regions.
(c) Electrophilic model.The parameters in the orange band of this model explain current electron g−2 uncertainty.These scalars give rise to two photons via decay and e + , e − via splitting for m ϕ < 0.001 GeV, and e + , e − only via decay for m ϕ > 0.001 GeV.The shaded regions represent 90% CL bounds from E137 [69], E141 [70], E774 [71], Orsay [72], NA64 [73], and HB stars [41] and SN cooling [74].ture.We see that the forward DUNE ND can continue to probe those parameters with short lifetimes even if a larger number of events are required to determine the sensitivity, i.e., the expected sensitivity reaches are not very sensitive to the underlying background assumption.
Since mediators with strong couplings tend to decay rapidly, the range of couplings in the ceiling of a sensitivity lays an upper bound to the coupling strength such that 3 of them make it to the detector.This limit depends on the mass and momentum of the mediator, and the distance between the source and detector [81].The lower edge, on the other hand, depends on statistics of mediators produced at the target which approximately scales with the square of the coupling [81].This depends on scalings, such as the number of POT, the branching ratio of the production process, etc.If we increase the required statistics by looking for the number of events greater than 3, the minimum coupling required would increase as compared to the case of 3 events, therefore pushing up the lower edge of the sensitivity.Figure 10b is an example that shows the sensitivity contours of the muonphilic scalar model at DUNE ND for 10 events and 100 events.It clearly shows that the contours do not shrink much as we increase the number of considered events.The change in number affects the lower edge of the sensitivity plot only.Since we would look for more statistics in the presence of backgrounds, they would affect the lower limits rather than the ceiling.This also implies that the presence of backgrounds would not compromise the ability of these experiments to probe parameters in the g − 2 band, as seen in the example plot.
The minimum threshold kinetic energy (E t ) required to identify signals at the detector plays a vital role in arriving at sensitivity plots of certain models at certain detectors.We notice that this constraint does not reduce the extent of the sensitivity curves of the muonphilic scalar model as the majority of the mediators are produced from processes that favor high-energy mediators.However, we see that this plays a role in the sensitivity of the electrophilic scalar model, especially for masses in the keV range as scalars produced from Compton-like scattering, being the dominant one, are lower in energy, and so are the final states.Since the sensitivity estimates for DUNE ND give us the most insight into unexplored parameter spaces, we investigate the effect of the energy thresholds at the DUNE ND in this study.Figure 10c demonstrates the effects of this by showing the contours without cuts versus those with 10 MeV and 30 MeV cuts.As expected, we see the reduction in sensitivity space coverage by a few factors toward the lower-mass regime.Similar to kinetic energy thresholds, the imposition of a minimum angular separation between the decay products can play a vital role in the sensitivity plots as they affect the ceiling of sensitivity reaches.We notice that the reduction factor of coupling values in the ceiling can vary from O (10) to O(1) based on the threshold value and the mass of the mediator.
In summary, Fig. 10 demonstrates multiple ways in which sensitivities are altered with experimental constraints.A more detailed analysis of detector responses and backgrounds for LArTPC-type detectors will certainly allow for more precise sensitivity estimates.

VII. VECTOR MEDIATORS
As discussed earlier, the above analysis is not limited only to scalar mediators.They can be extended to spin-1 mediator models (gauge bosons) as well.They can be greatly produced from charged mesons present at the target.As an example, let us consider the U (1) T 3R model that appears in the context of left-right symmetric models [82,83].This model contains a gauge boson, Z ′ , that couples to the right-handed fermions of one particular generation.Although there exist new fields in this anomaly-free model which is a spontaneously broken model at a low energy scale, e.g., low mass scalars and dark-matter candidates, we consider the effects of the gauge boson with visible decay modes only where the gauge boson couples to the right-handed muon, charm, and strange quarks.
These gauge bosons can be heavily sourced from threebody decays of charged mesons K + → µ + ν µ Z ′ , emanating from the muon leg µ + .Since these gauge bosons couple only to the right-handed component of the muon, the muon's helicity must be flipped, thereby suppressing this three-body decay by the muon mass.Despite this condition, we observe an enhancement effect, similar to the scalar three-body decay case.This is due to the existence of the longitudinal polarization mode of vector gauge bosons that allows for an enhancement proportional to 1/m 2 Z ′ [9].Since these gauge bosons can mix with the SM photon kinetically, they can couple to electrons with strength ϵe where This expands the horizon of gauge boson production including neutral meson decays, π 0 /η → γZ ′ .Similarly, gauge bosons can also be produced via Compton-like scattering, γe − → Z ′ e − [84] when photons hit the target electrons.The production mechanisms of electrophilic scalars can be applied to Z ′ gauge bosons as well, i.e., pair annihilation e + e − → γZ ′ [85] and resonance production e + e − → Z ′ [85].Additionally, Z ′ s can appear from electron/positron bremsstrahlung e ± N → e ± N Z ′ [85,86].However, it is important to note that these processes occur through γ − Z ′ mixing, which is of the order O(10 −2 ).Therefore, the flux of gauge bosons from three-body decays of charged mesons is much larger as they are produced via the direct coupling to the lepton.It is important to note that the three-body decay must satisfy the upper limit of the charged kaon/pion branching ratio.
We probe the sensitivity of the U (1) T 3R gauge boson through visible decays into electrons and positrons via the kinetic mixing loop.This decay width is given by (a) Sensitivity plot for HPS for 3 events as well as 100 events for DUNE ND and ICARUS.In this plot, we include 90% CL bounds imposed by CHARM [63], LSND [64], PS191 [65], NA62 [66], E949 [67], and 95% CL bounds from µBooNE [29], and LHCb [68].The shaded regions represent 90% CL bounds from E137 [69], E141 [70], E774 [71], Orsay [72], NA64 [73], and HB stars [41] and SN cooling [74].We see that the production of these gauge bosons from charged mesons dominates over neutral mesons despite constraints on the three-body branching ratio.The magnetic focusing horn system facilitates this production mechanism and hence, we obtain a great amount of sensitivity for weak couplings.Due to constraints on the upper limit of the three-body branching fraction, we are unable to excavate higher couplings through this production mechanism, which are mostly constrained by experiments such as NA64, E774, etc.Thus, in this U (1) T 3R model, we see that the production of these gauge bosons from charged mesons can give us insights into the lower coupling range that are unexplored parameters, as supported by our sensitivity study under the assumption of negligible backgrounds (see Fig. 11).Since the branching ratio of a charged meson into a gauge boson is proportional to g 2 T 3R /m 2 Z ′ , the floor of the sensitivity has a linear nature, with a positive slope.This leads to a unique shape of the sensitivity as compared to those in the excluded regions in which the gauge bosons are produced through two body decays of neutral mesons (where the branching ratio is roughly proportional to g 2 T 3R ).Finally, we attempt to apply these production and detection mechanisms to U (1) Li−Lj gauge bosons.Since these gauge bosons couple to neutrinos as well, their lifetimes at MeV-to-sub-GeV facilities are not long enough to suggest considerable value-added insights into unconstrained parameters.

VIII. CONCLUSIONS
In this study, we explore three dark-sector models, Higgs Portal Scalars, muonphilic scalar models, and electrophilic scalar models, at neutrino experiments by utilizing the vast flux of mesons, photons, and electrons that are produced at the neutrino target.The neutrino experiments considered in this study are the finished Ar-goNeuT, ongoing MicroBooNE, SBND, ICARUS, and the upcoming DUNE experiments.We have also demonstrated an example of a vector muonphilic mediator, which is the gauge boson in the U (1) T 3R model.The magnetic horns at the targets of these experiments have been utilized to maximize the production of the mediators from charged mesons along the beam direction.The choice of models and mediators enables us to understand the domination of one production mode over the other.Although three-body decays of mesons dominate for models with muonic couplings, and two-body meson decays dominate the HPS model, Compton-like scattering and Primakoff production from photons are most dominant for models with electron couplings.Through these production mechanisms, the potential for high-energy mediators to reach the detectors increases, especially at forward detectors.This allows us to explore the g−2 regions along with large regions of unexplored parameter space in the laboratory experiments at the ongoing/upcoming facilities, especially for the muon.
Since these mediators produce visible signals in the detector through multiple scattering and decay mechanisms, we were able to get an all-inclusive sensitivity plot demonstrating the potential of probing regions of parameter space that are still unexplored by experiments.The inclusion of the energy-dependent Bethe-Heitler scattering process and the decay process allowed us to access unique lepton-antilepton signatures for flavor-specific mediators.
A rigorous background analysis would allow us to estimate more realistic exclusion limits.A better understanding of the background would include angular separation capabilities along with exact energy thresholds.For example, Ref. [29] has done a study for electronpositron events at MicroBooNE from which we find an estimate of around 30 background events.However, it is important to note that the ceiling of sensitivities does not change appreciably, as demonstrated in Fig. 10, since they are obtained from extremely high-energy mediators.The number of events calculated in this study may have uncertainties that vary between 0.3−3, due to the charged meson fluxes.This can vary the new physics couplings in the sensitivity estimates (Figs. 9 and 10) by factors ranging from 0.7 to 2. However, these effects are visible only on the floor of the sensitivity plots.
Our studies based on the example models can straightforwardly be extended to many more scenarios with bosonic mediators.It must be noted that the positively charged pions and kaons are along the beamline due to the magnetic horn effect.As for the positrons and photons, they have an angular spread, but the most energetic ones are forward in nature.
We employ GEANT4 to obtain the charged meson fluxes for each of the targets without the effect of the magnetic field.To simulate the magnetic horn effect on charged mesons, we apply energy and angle cuts for the charged mesons to match the resulting neutrino fluxes with the reported ones.At the BNB target from the neutrino flux, we follow the prescription that only those charged mesons with a total energy greater than 750 MeV, and an angle between [0.03, 0.2] radians to the beamline, pass through the magnet [87].For LBNF, we follow the empirical models developed in Ref. [88].We first collect only those particles with kinetic energies greater than 100 MeV and with polar angles between [0.01,1] radians (z-axis is the beamline), assigning an additional weight of 3 for every particle with energy less than 10 GeV and 0.2 for those greater than 10 GeV.Using this prescription, we roughly reproduce the muon-neutrino flux [89] at the DUNE ND.Similar to the flux reported by the DUNE collaboration, our neutrino flux drops down for neutrino energies greater than 10 GeV.On comparing the fluxes for energies below 10 GeV, we notice that there are overestimates and underestimates up to a factor of 3.However, we find that these mismatches do not appreciably affect our sensitivity reaches, especially the ceiling.)) (B4) Although the kaon has the four-momentum given by (E K , ⃗ p K ), we will work in the rest frame of the decaying kaon after which we can boost the momenta to the lab frame.For a three-body decay, there are 9 (3 × 3) variables determining the momentum (and hence the energy) of the final states.The four constraints on the energy momentum of the final states reduce the degrees of freedom to 5. Out of these 5, three of them determine the plane on which the decay occurs.For the given scalar interactions, the matrix element does not depend on the choice of planes.Hence, they can be integrated, which includes a factor of 2(2π) 2 .Therefore, the remaining two degrees of dalitz variables in this case.The other five invariants, t a2 , t b2 , t a3 , t b1 and s 13 can be expressed in terms of those in Eq.C2.In order to carry this out, we define the two kinematic functions below.The phase space integral can be expressed as 3 ) 8s 23 (C6) Therefore, the differential cross-section can be expressed as , which are the angles in the rest frame of particles 2 and 3, which is l + , l − .These angles can be expressed in terms of the Dalitz invariants as below.The momenta can be obtained from the above energies.By integrating these variables, we obtain the total crosssection of the 2-3 Bethe Heitler scattering process for a given energy and mass of the incoming scalar.The total number of Bethe-Heitler scattering events is given by Where N A , ρ T , and A T are the Avogadro number, density, and atomic mass of the detector respectively.L det is the length of the detector.

Figure 2 :
Figure 2: Feynman diagrams depicting the production of scalars from charged kaons.(a) Production of both HPS and flavor-specific scalars from the charged lepton leg.(b) HPS produced from the W boson leg.(c) Feynman diagram for two-body decay of K + to π + and HPS S which couples to an intermediate top quark.

Figure 3 :
Figure 3: Feynman diagrams depicting the production of scalars from photons.(a) Primakoff scattering of a photon to produce HPS as well as flavor-specific scalars.(b) Compton-like scattering of a photon to produce an electrophilic scalar.

Figure 6 :
Figure 6: Number of scalars produced from each source per unit coupling.Left: Higgs Portal Scalars, Center: Muonphilic scalars, and Right: Electrophilic scalars

Figure 7 :
Figure 7: Feynman diagrams illustrate the various scattering channels that mediators can undergo once they reach the detector.(a) Inverse Primakoff process where incoming scalars (ℓ = e, µ) scatter into a single photon.(b) Bethe-Heitler splitting process.(c) Inverse Comptonlike scattering of ϕ e into an electron and a photon.

Figure 8 :
Figure 8: Number of scalars that reach the detector as a function of mass for various couplings.Left: Higgs Portal Scalars, Center: Muonphilic scalars and Right: Electrophilic scalars

Figure 9 :
Figure 9: 95% CL lines for our three benchmark scalar mediator models under the assumption of zero backgrounds.
(b) Parameters for 10, 100 muonphilic scalar decays at DUNE ND.The 90% CL bounds from E137 [38], 95% CL bounds from B factories, g − 2 excluded regions, and constraints from SN cooling are depicted in the shaded regions.(c)Sensitivity plots for electrophilic scalars at DUNE ND where the minimum threshold kinetic energy (Et) required to identify the decay products is 0 MeV, 10 MeV, and 30 MeV.

Figure 10 :
Figure 10: Sensitivities where experimental and detection constraints are considered.

Figure 11 :
Figure 11: 95% CL sensitivity plot (with no backgrounds) of the U (1) T 3R gauge boson with only visible decays.
Appendix A: Parent particle fluxes