Probing the Higgs Portal at the Fermilab Short-Baseline Neutrino Experiments

The Fermilab Short-Baseline Neutrino (SBN) experiments, MicroBooNE, ICARUS, and SBND, are expected to have significant sensitivity to light weakly coupled hidden sector particles. Here we study the capability of the SBN experiments to probe dark scalars interacting through the Higgs portal. We investigate production of dark scalars using both the Fermilab Booster 8 GeV and NuMI 120 GeV proton beams, simulating kaons decaying to dark scalars and taking into account the beamline geometry. We also investigate strategies to mitigate backgrounds from beam-related neutrino scattering events. We find that SBND, with its comparatively short ${\cal O}(100\ {\rm m})$ baseline, will have the best sensitivity to scalars produced with Booster, while ICARUS, with its large detector volume, will provide the best limits on off-axis dark scalar production from NuMI. The SBN experiments can provide leading tests of dark scalars with masses in the 50 - 350 MeV range in the near term. Our results motivate dedicated experimental searches for dark scalars and other long-lived hidden sector states at these experiments.


I. INTRODUCTION
Light weakly coupled hidden sectors may play a role in addressing some of the outstanding puzzles in particle physics and cosmology, such as dark matter, neutrino masses, the matterantimatter asymmetry, the hierarchy problem, and inflation [1][2][3][4]. They also provide an interesting physics target for a variety of intensity frontier experiments. One well known example is proton beam fixed target experiments, including those used to study neutrino oscillations. In these experiments, a beam of relativistic hidden sector particles produced in the primary proton-target collisions passes through a downstream detector and decays to visible Standard Model (SM) particles, providing a distinctive signature that can be discriminated from beam-related neutrino and cosmic-ray induced backgrounds.
The Short-Baseline Neutrino program at Fermilab utilizes three liquid argon time projection chamber detectors -SBND [5], MicroBooNE [6], and ICARUS [7] -situated along the Booster Neutrino Beam (BNB) [8]. While the primary physics goals of the SBN program include searches for eV-scale sterile neutrinos (motivated by the LSND [9] and Mini-BooNE [10][11][12] anomalies) and measurements of neutrino-nucleus scattering cross sections, these experiments are also expected to have significant sensitivity to hidden sector particles in the MeV-GeV range. In this paper we study the capability of the SBN experiments to search for dark scalar particles, S, interacting via the Higgs portal.
We consider production of dark scalars from collisions of both the 8 GeV Booster and the 120 GeV Main Injector protons. Concerning the latter, MicroBooNE and ICARUS are located approximately 8 • and 6 • off-axis with respect to the NuMI beam, implying that the associated flux of dark scalars passing through these detectors can be substantial. We perform a careful simulation of the BNB and NuMI beamlines that includes the targets, magnetic focusing horns, and absorbers. While we also attempt to take account of various experimental factors such as particle identification and reconstruction efficiencies, detection thresholds, and measurement resolutions, a complete characterization of these parameters is still under active experimental study, so our assumptions in this regard must be considered preliminary. We consider two primary dark scalar signal channels in detail, 1) the S → e + e − channel, relevant for scalars below the di-muon threshold, and 2) the combined S → µ + µ − , π + π − channels which are important for somewhat heavier scalars. We study the beam-related backgrounds to these signatures and design search strategies to efficiently separate them from the signal. We also discuss the potential role of timing measurements as a discriminator for heavy O(100 MeV) scalar particles.
We find that the SBN experiments have significant sensitivity to dark scalars with masses below a few hundred MeV. In particular, SBND, being in close proximity to the BNB target, has the best sensitivity to scalar particles originating from the Booster beam. On the other hand, for scalars produced along the NuMI beamline from the Main Injector, we find that ICARUS, with its large detector volume, will have the leading sensitivity. Notably, kaons decaying at rest (KDAR) after stopping in the NuMI absorber could produce monoenergetic scalars entering ICARUS at a different angle from neutrinos produced in the target, providing a novel signature. In both cases, our projections cover regions of scalar mass -mixing angle parameter space extending beyond existing experimental limits. Our results motivate dedicated searches for dark scalars and other light weakly coupled particles at the SBN experiments. We note that previous studies have examined the sensitivity of the SBN experiments to heavy neutral leptons produced with the Booster beam [13] and the darktrident signature from dark matter produced in the NuMI beam [14]. Furthermore, a number of studies have explored the reach of accelerator-based neutrino experiments to long-lived and other exotic particles such as dark matter in recent years; see e.g. Refs.  and the recent white paper [54].
The outline of this paper is as follows: In Section II we provide a brief review of scalars interacting through the Higgs portal. Section III provides an overview of the SBN facillities.
Section IV describes our simulation pipeline for the signal and beam related backgrounds, as well as our analysis strategy for the channels under consideration. Our results are presented in Section V, and we conclude in Section VI with some discussion and outlook. An appendix contains some results of our simulation validation.

II. HIGGS PORTAL
We extend the SM by a new singlet scalar particle, S, which interacts through the Higgs portal. There are two possible renormalizable portal couplings in general, We will assume the dark scalar acquires a small coupling to SM particles via mass mixing with the Higgs. This always occurs if the super-renormalizable term A in Eq. (1) is nonvanishing. However, even if A = 0, the dark scalar can acquire a vacuum expectation value for appropriate choices of the scalar potential, which in turn leads to mass mixing. After diagonalization, there are just two parameters relevant for our phenomenological study, namely the physical mass m S of the dark scalar and the scalar-Higgs mixing angle θ: Given the tight experimental constraints on the mixing angle, we will be working in the regime θ 1. Additional interactions between the scalar and the Higgs (e.g., hSS), while important for higher energy collisions at the LHC, will not be relevant for the SBN experiments. As S inherits its couplings in Eq.
(2) to SM particles from the Higgs, many important phenomenological considerations can be recycled from light Higgs boson studies carried out decades ago [55].
Starting from the Lagrangian (2) we can compute the decay rates of the dark scalar. The partial width to charged leptons is given by For m S > 2m π hadronic decays of the scalar become important. The description of hadronic decays is complicated by strong interactions and resonance effects leading to sizable theoretical uncertainty in the decay width for scalar masses of order 1 GeV. The two pion decay can be described as 1 [57,58] where we follow the treatment in Ref. [59] for the form factor G π (s), which uses the results of Ref. [60]. At low energies, the form factor is given by chiral perturbation theory, G π (s) = 2 9 s + 11 9 m 2 π . Within the model described by Eq. (2), any other SM decays have negligible branching fraction. In particular, there are no other dark sector states into which the scalar S might decay. In Figure 1 we show the scalar branching ratios as a function of its mass, as well as isocontours of the scalar decay length in the m S − θ plane. 1 Results for scalars with general couplings to SM particles can be found in Ref. [56]. For the SBN experiments the most important production channel will be through kaon decays, K → πS, which proceed via the one-loop penguin and scalar bremsstrahlung processes [17,55,59,61,62]. The dominant contribution comes from the W boson -top quark penguin shown in Fig. 2, leading to the partial decay width, with Γ(K 0 L → π ± S) Γ(K ± → π ± S). The branching ratio for scalars produced in K 0 S decays is smaller by several orders of magnitude compared to those from K ± and K 0 L , owing to suppression from the small CP-violating phase and the substantially shorter K 0 S lifetime. Besides kaon decays, B mesons have an even larger branching ratio to S and furthermore allow for production of heavier scalars [17,23,[63][64][65]. However, while such decays are important for higher energy facilities such as the SPS and LHC beams at CERN, they will not be relevant for the lower energy Booster and Main Injector beams due to the substantially lower B production rate. A more promising process at these energies is scalar production through proton bremsstrahlung, which also allows for production of heavier scalars. We have estimated the sensitivity from proton bremsstrahlung using the results of Ref. [66], finding that it leads to no additional sensitivity beyond that from Kaon decays. Additional sources of low mass scalar production at these energies include rare π ± , η, η decays, although they are subdominant to kaon decays and will be neglected in our study.
There are a number of experimental and astrophysical constraints on the Higgs portal scalar parameter space, but we will defer a discussion of these to Section V below. We now move on to discuss the SBN experiments.

III. SBN EXPERIMENT SETUP
FIG. 3. Coordinates in the "BNB frame" defined by a left-handed coordinate system with the origin at the BNB target, the z-axis pointing to center of the SBN detectors, and the y-axis pointing vertically upwards. We show the coordinates of the BNB and NuMI targets, NuMI absorber, and the SBND, MicroBooNE, and ICARUS detectors in both aerial (left) and elevation (right) views.
The sensitivity of an experiment to decays of long-lived dark sector particle decays depends primarily on the number of such particles that can be made to pass through the detector and the ability of the detector to distinguish the decays from backgrounds. In both these regards, the Short-Baseline Neutrino (SBN) program at Fermilab offers strong capabilities for near future sensitivity. For Higgs portal scalars the dominant production mechanism comes from kaon decays. The Booster Neutrino Beam (BNB) and the Neutrinos at the Main Injector (NuMI) beam offer a large production of both focused and stopped kaons with several detectors in the 100 to 1000 m range of their targets. Among these detectors, the Liquid Argon Time Projection Chamber (LArTPC) detectors of the SBN program, namely the Short-Baseline Near Detector (SBND), MicroBooNE, and ICARUS-T600, will have low thresholds and excellent particle identification. We therefore focus our study on the signals of beam-produced dark scalars decaying in these detectors. Figure 3 shows aerial and elevation views of the Fermilab campus, including the locations of the BNB and NuMI targets, NuMI absorber, and the SBND, MicroBooNE, and ICARUS detectors.
In this section, we discuss the capabilities of both the beams and SBN detectors in turn. The Fermilab accelerator neutrino program currently consists of the aforementioned BNB and NuMI beams. Each has its own advantages and disadvantages for the purposes of dark scalar searches and will have its best sensitivity at different SBN detectors. We therefore discuss the two in turn. Their properties are summarized in Table I. for around 4 years at that intensity for SBND and ICARUS. We take this number, 1.3×10 21 , as our benchmark for our BNB-based studies. Charged mesons are focused by a magnetic horn to deliver a beam enhanced in either neutrinos or anti-neutrinos.

NuMI Beam
The NuMI beam is made by striking 120 GeV kinetic energy protons on a graphite target. The beamline is roughly north pointing along the axis from Fermilab to the MINOS far detector in Minnesota. The SBN detectors are not along the beam axis, but are slightly off axis. The beam does however pass at small angle with respect to both MicroBooNE and ICARUS-T600. Given the higher energy of the beam and the angular spread in the produced scalars, these detectors can have significant sensitivity to scalars produced with the NuMI beam. The MicroBooNE detector is located a distance of roughly 685 m from the target at an angle of 0.14 radians (7.8 • ), while the ICARUS-T600 is 803 m at an angle of 0.097 radians (5.5 • ). The accelerator delivered 5.7 × 10 20 POT in the 2018 fiscal year [67], for a projected total of close to 3 × 10 21 expected from running for 5 years, which we take as our benchmark for NuMI-based studies. The NuMI beam has two horns whose position and current can be adjusted depending on the desired energy and type (neutrino or antineutrino) of the beam. The actual running mode will be determined by NOνA physics requirements.
However, because an order-1 fraction of our signal will come from decays of neutral mesons which are unaffected by the magnetic horns, we do not expect the mode to significantly affect the final results. We assume the standard "medium energy" mode horn locations, which have been used for the past several years [68], and neutrino mode running with a 200 kA forward horn current configuration below.

B. SBN Detectors
The SBN program consists of the SBND, MicroBooNE and ICARUS-T600 LArTPC detectors positioned along the BNB. For the most part, the detectors function in similar manners, with similar TPC designs. Though we assume the same detection capabilities for the three SBN detectors in this study, we do account for their difference in volume and location, as outlined in Table II. MicroBooNE is currently operating, while SBND and ICARUS-T600 are expected to be collecting data by the end of 2019 and beginning of 2021 respectively [69,70]. The volume and mass are particularly relevant for signal and background rate calculation respectively as described below. Note that while we list the full dimensions for a box shaped detector volume, we make the simplifying approximation that the rates for the signal do not depend on the detailed geometry of the detector and perform our calculations using a spherical detector with the correct fiducial volume. This approximation is valid in the limit that the detector dimensions are much less than the distance that the scalar travels before decaying. Since the detector dimensions are O(1 m), while the distance the scalar travels is O(100 m), this approximation is valid at the 1% level.
We discuss detector effects and event reconstruction in detail below in Section IV C, after further discussing the relevant signals and backgrounds.  [5]. The dimensions and masses assumed are those of the fiducial volume for the ν µ analyses.

A. Signal event generation
To simulate the proton collisions from the Booster and NuMI beams, we employ the g4bnb [71] and g4numi [72] codes. Each beamline is assumed to be configured in neutrino mode. The simulations, based on Geant4 [73], incorporate the geometry of each beamline including the targets, focusing horns and decay volumes. They take into account the full evolution of the particles produced in the primary collisions of protons with the target, including not only secondary meson production but also the subsequent interactions with the beamline elements. In both codes, whenever a neutrino is produced in a particle decay, information about the neutrino momentum and position, as well as its ancestor, is stored using the dk2nu [74] format.
We use the ancestry information provided in the g4bnb/g4numi output to extract the momentum and position of each kaon that decayed to a neutrino, at the time of decay. In and enter the Booster beamline detectors at an angle that is quite different from those of their counterparts arising from kaons decaying in the target. We will discuss this possibility further In Sec. IV D.
Using the kaons, we then simulate K → πS decays for different scalar masses, yielding the phase space distribution of scalars produced in each beamline. We consider only decays of K ± and K 0 L , as the K 0 S provides a negligible contribution due to its significantly smaller branching ratio to scalars. The mixing angle θ sets the overall normalization of the distribution, i.e. the number of scalars produced per POT.
Finally, for each of the SBN detectors, we consider the scalars from each beam that would pass through the detector volume. If a scalar with velocity β and rest lifetime τ would enter and exit a given detector at distances L in and L out from its production point, the probability that it decays inside the detector is  Table II), more scalars from NuMI decay in ICARUS than in the other SBN detectors. Typically, the scalar parent is a highly boosted kaon, so the scalar and its decay products are all well-aligned with the beam direction. In Sec. IV D, we use these characteristic features of the scalar decay products to design cuts to reduce backgrounds from neutrinos and cosmic rays.
We have checked our pipeline in two different ways.

B. Background Generation
There are several sources of potential background events for the signal we are considering, though all are in principle reducible. As the SBN detectors are operated on the surface, cosmic ray induced backgrounds are a significant worry. On the other hand, it is rare for cosmic rays to mimic the signal we are considering, as any tracks or showers should come in pairs and generally appear in the middle of the detector. Furthermore, systems are being designed specifically to tag and eliminate cosmic ray backgrounds at MicroBooNE [76], while SBND and ICARUS are also being designed with dedicated cosmic ray taggers as well [5].
Beyond such a tagger, a combination of timing and angular cuts can further reduce this background, though it could in principle still hinder reconstruction of signal events. We do not consider this background further in this study.
The remaining backgrounds are due to beam neutrinos. Both interactions in the dirt and in the detector volume can contribute, though we only consider interactions in the fiducial detector volume, which should dominate for the analysis we have developed. Even this dominant background cannot fully mimic our signal and several cuts can be used to drastically reduce it, as discussed below.
To generate the background, we use the neutrino flux predictions presented in [77]. The fluxes are used for input for generation of neutrino interactions on 40 Ar using GENIE v3.
The π 0 in the GENIE output are decayed to photons. These photons, along with any other final state particles reported to GENIE are smeared and thresholds are applied according to the description in Section IV C below.

C. Detector Effects
In this section, we discuss our assumed detection thresholds, energy and angular resolutions, and particle identification capabilities.
The target volume consists of liquid argon, which provides a dense medium for neutrino in- The deposited energy can be determined from the charge collection, allowing a reconstruction of the energy deposited by a track per unit length, which in turn allows for strong particle identification and calorimetry. Given all of this information, the full four vector of the particles produced through an interaction in the detector volume can be reconstructed with excellent angular resolution and good energy resolution.
Many of the capabilities of the SBN LArTPC detectors are under active investigation and are not known. For the purposes of our study, we assume certain figures of merit that we believe are reasonable. As understanding of the detectors improves, we expect to have updated sensitivities. We now discuss the properties that we assume, as well as some motivation for these assumptions.
Particle identification is one of the main advantages of LArTPC technology. That said, certain particles can be rather tricky to distinguish. The particles that are stable on the scale of the position resolution of LArTPC detectors in the 100 MeV energy range of interest are e ± , γ, p, µ ± , and π ± . The first two are generally referred to as shower particles as they create electromagnetic cascade showers above 33 MeV, while the last three are referred to as track particles, as they appear as single tracks. It is worth noting that, in some instances, the distinction between track and showering particles is not 100% accurate. We assume however that this distinction is perfectly accurate in the analysis below.
The signal of interest consists of an isolated e + e − , µ + µ − or ππ pair. For the ππ case, we focus on the simpler and more prevalent π + π − topology, though a 2π 0 decay is also possible. Since no neutrino interactions from the beam produce precisely this final state, the backgrounds are, in principle, reducible. There are, however, neutrino interactions that can fake our signal in a LArTPC. Similar arguments eliminate the subdominant backgrounds due to charged-current electron production in association with a photon. Therefore, for a sufficient isolation cut, the search for e + e − pairs should be effectively background free.
For the remaining searches, the primary background comes from pions and muons. Protons, muons, and charged pions all leave tracks, but protons are distinguishable by their distinct energy deposit per unit length compared with track length. It is worth noting that some of the muons and pions from our signal have large energies and are expected to not be contained by the SBN detectors. From Geant4 simulations, we expect a muon of energy 1 GeV to have a track length of around 4 m, for example. A distribution of muon candidate energies will be shown below in Fig. 13. This does not necessarily mean that the event cannot be reconstructed, but rather makes reconstruction more challenging. For the present analysis, we do not take containment effects into account though we do not expect this to significantly affect our sensitivity estimates. Muons and charged pions are in general quite challenging to tell apart as they deposit energy at similar rates. Charged pions leave larger energy deposits at the end of their tracks due to their additional nuclear interactions, but the ability to utilize this has not yet been quantified. We therefore conservatively assume that they are indistinguishable. Note that the charge of the muon and pion should be distinguishable based on capture of µ − that is not possible for the opposite charge. The effectiveness of this procedure has not yet been determined, so we do not attempt to use the muon or pion charges. Since pions and muons are assumed to be indistinguishable, the backgrounds to searches are typically somewhat large, coming from µ ± π ∓ inelastic charge current events and π ± π ± deep inelastic events. We make use of kinematic cuts as described below to reduce these backgrounds.
Next, we discuss the kinetic energy thresholds for observing particles. Thresholds in LArTPCs are typically in the 10-100 MeV range, though they depend on the specific detector, context and particle. For on-beam analyses, the threshold is set by requiring that a sufficient number of hits are seen in the collection plane. MicroBooNE typically requires at least 20, setting thresholds of roughly 80 MeV for protons, 40 MeV for π ± /µ ± , and 20 MeV for EM showers, though muons are generally required to have longer tracks. Reconstruction of shorter proton tracks should be possible as demonstrated by ArgoNeuT [78], allowing for thresholds of 21 MeV for protons and 10 MeV for π ± /µ ± . The DUNE CDR quotes a comparable set of thresholds used for its Fast Simulation [79]. Light collection is crucial for off-beam analyses and helps with energy reconstruction in all cases. We assume that light collection thresholds are met if the tracking thresholds are met.
We assume a threshold of 30 MeV for both electrons and photons, motivated by the requirement that they have at least one hard interaction in the detector. We also assume a threshold of 30 MeV for muons and pions, as motivated by a requirement of at least 10 hits in the TPC. While this is slightly more aggressive than the current MicroBooNE assumption, improvements on reconstruction should be made in the near future. We also reconstruct protons with a threshold of 20 MeV for the purposes of vetoing. While this is lower than the nominal reconstruction threshold, it is noted that proton reconstruction is not required; we only need to determine whether a proton was present in the event.  [13]. A recent MicroBooNE thesis has found even better resolution for contained muons and good resolution on the neutral pion peak [80,81]. A study of momentum reconstruction for muons from multiple Coulomb scattering, suitable for exiting muons found resolution in the 10-25% range across different energies [82]. It is worth noting that some of the region of interest for our searches involves high energy particles that may not be contained in the detector, which would degrade both resolution and particle identification capabilities. Based on these numbers, we assume 10% resolution on electromagnetic showers, 3% momentum resolution on muons and pions below 300 MeV kinetic energy as these are likely to be contained, and 10% momentum resolution on muons and pions above this threshold as they are likely to exit. A lower resolution would require a wider mass window around the candidate scalar mass, increasing the background somewhat, but not drastically changing our results.
Finally angular resolution should be excellent for tracks, given the excellent position resolution. We adopt a universal angular resolution for particles of 0.03 rad or 1.73 • as in Ref. [5].
Our analysis assumptions are shown in Table III. We apply these thresholds and smearing at the level of final state four vectors from the event generator chains, noting that these numbers are subject to change as LArTPCs are better understood. We further veto on additional photons and electrons above 2 MeV, as well as any additional pions or muons, which will also lead to visible electrons when they decay.
We make one final comment regarding neutrons and recoiling nuclei. Both can lead to additional energy deposits in the detector, the neutrons scattering widely throughout the volume [83] and the recoiling nucleus activity centered near a neutrino interaction point [70].
Decaying mediators do not leave such signals in the detector, so vetoing such activity would be an excellent way to reduce backgrounds. Nevertheless, there is no official study of these effects, so we simply ignore neutrons and recoiling nuclei in our analyses.

D. Analysis
Below, we describe our analysis strategy for the various signal channels. For m S < 2M µ , the dominant decay channel of the dark scalar is to an e + e − pair. There is essentially no irreducible background to this signature, but it is somewhat challenging to Distributions are shown using smeared candidates above threshold in events with exactly two candidate electrons, but absent any additional kinematic cuts.
determine the reducible backgrounds. We begin by selecting events with no reconstructed muons, charged pions, or protons above the thresholds in Table III. By further selecting events with exactly two electromagnetic showers, we eliminate a large fraction of the background, as we can then require that the two showers originate at the same vertex. To estimate the effect of this requirement, we require that the signal candidates have a separation of at least 10 • so that their vertex of production is well-reconstructed. This is sufficient The second, wide angle bump is coming from scalars from KDAR, which are lower energy and lead to more widely separated decay products.
to effectively eliminate all of the background as the conversion points for diphoton events will be well separated and single photon conversions are strongly peaked at narrow opening angle, as explained above in Sec. IV C. A cut on the angle between the reconstructed scalar momentum and the beamline could be helpful in reducing the background from cosmics, but is likely not necessary for the e + e − case and we do not apply such a cut. An additional search for a bump in the invariant mass of the e + e − would lead to a measurement of the mass of the scalar. A diagram for a candidate background event is shown in Fig. 6. The event distributions in the key variables on which we cut are shown in Figs. 7, 8, and 9.
The distributions are drawn after the reconstruction and smearing procedure described in Sec. IV C before applying the isolation requirement and after requiring reconstruction of a single e + e − candidate pair. No additional kinematic cuts are applied at this level. Note that the distributions in candidate angular separation and reconstructed scalar-beamline angle show a second peak corresponding to kaons decaying at rest to scalars, which frequently occurs in the absorber.  2. µ + µ − and π + π − Above 210 MeV, muon and pion channels quickly dominate. As discussed above, distinguishing µ ± and π ± is difficult. There is some recent work along this direction using machine learning techniques, for example in Ref. [84]. Since the ability to distinguish these particles is not fully tested, we do not attempt to do so. We rather assume µ ± and π ± are Distributions are shown using smeared candidates above threshold in events with exactly two candidate muons, but absent any additional kinematic cuts. indistinguishable. If they could be distinguished, it is worth noting that there is essentially no irreducible µ + µ − background, so significant gains could be made. We veto events with any reconstructed EM shower or proton above the kinetic energy thresholds in Table III. We demand that there are exactly 2 candidate µ ± /π ± tracks. A cut of E > 2.5 GeV is imposed on reconstructed scalars as we determined that this helps to separate signal from background for target-produced scalars. We place an angular cut on the reconstructed scalar of 4 • with respect to the beamline. We lastly place a cut on the invariant mass of the reconstructed scalar of 40 MeV around the scalar mass. While the scalar mass is not known a priori, a scan over different scalar mass windows could be performed as we do in this study. A diagram for a candidate background event is shown in Fig. 10 and the event distributions in the key variables on which we cut are shown in Figs. 11 and 12. As in the e + e − case, the distributions are drawn after the reconstruction and smearing procedure described in Sec. IV C and after requiring reconstruction of a single µ + µ − candidate pair. No additional kinematic cuts are applied at this level.

Kaon decay at rest
One intriguing possibility using the NuMI beamline is to make use of the fact that the detectors are off of the beam axis. The portion of the beam that does not interact with the target is absorbed farther down the beamline. The line between the absorber and ICARUS is about 110 m long and forms an angle of about 46 • with respect to the beamline. This raises the possibility of searching for particles coming from kaons decaying at rest (KDAR) in the absorber. In fact, MiniBooNE has recently measured monoenergetic muon neutrino charged current events arising from KDAR in the NuMI absorber [85], and one can also envision searching for exotic particles such as dark scalars in this manner. These scalars would be monochromatic and all coming from a specific angle, leading to significant background reduction. Since the e + e − analysis above is already effectively background free, we do not search for KDAR specifically in our analysis. If the cosmic ray background proves to be more challenging than anticipated, searching for KDAR kinematics could prove helpful.
Furthermore, if a discovery of an excess is made, searching for KDAR would be an excellent tool for validation, as well as for probing the nature of the model due to the specific angle of 46 • for the reconstructed scalar and the monochromatic energy of the scalar. For µ + µ − and π + π − channels, the particles coming out of the monochromatic scalar decay would be very soft and challenging to reconstruct as they are produced near threshold in the scalar decay and the scalar itself is produced close to threshold at high mass. One could perhaps alleviate the reconstruction challenge by doing a simultaneous reconstruction of the two short candidate tracks, but development of such an analysis is beyond the scope of this work. Due to these challenges, we do not pursue the KDAR channel further for the µ + µ − or π + π − channels either. 4. Other channels: π 0 π 0 , γγ Events with two π 0 could be rather striking, with as many as 4 EM showers. The odds that all 4 showers are reconstructable is small and it is rather challenging to analyze such events. It is worth noting that the branching fraction of this channel is never dominant, either. We therefore defer further discussion of it to future analysis.
Rare decay channels of the scalar could be interesting to study as well, particularly in order to confirm the nature of the scalar in the event of a discovery. In particular, at loop level, decay to γγ is induced just as for the standard Higgs. Though this branching fraction is at the 10 −3 level, it is of great interest for confirming the nature of the scalar. Furthermore, non-minimal models could have a larger branching fraction to photons. For example, in a model with heavy charged states that do not couple to the Higgs, but do couple to the dark scalar, the diphoton coupling could be larger.

Discrimination with timing
We focus in this paper on the prospects for scalar decays during the times when neutrino events would also be occurring in the SBN detectors. However, unlike neutrinos, the scalar may travel at a speed sufficiently slow compared to the speed of light that its decay is delayed relative to the neutrino arrival. While the beam spills themselves are relatively long, 1.6 µs long at BNB and 10 µs at NuMI, each spill contains multiple shorter bunches of protons.
For instance, BNB employs proton bunches that are approximately 2 ns long, separated by 20 ns. The ns-level resolution of the SBN detectors can in principle be used to search for scalars decaying in between the times when neutrino events are expected [5]. We show the timing expected for the signal as compared to the neutrino background in Fig. 14 for a longlived 300 MeV scalar. Most of the signal events occur late enough relative to the neutrino background that they could be discriminated using timing cuts. However, the delay can be much longer than the bunch spacing. There are also additional peaks at late times coming from KDAR-induced scalars that can be produced in either the target or the absorber.
It would be interesting to further investigate the use of timing to look for delayed scalar decays, as has been discussed in the context of other light hidden sector models [13,86].
Depending on the mass, the length of the delay could vary considerably, and some care would be necessary in the case when scalars produced from one proton bunch decay in a detector after neutrino events from a later proton bunch. Furthermore, in order to perform such a study, new developments would be required in the data acquisition and triggering of the SBN experiments in order to look outside of standard timing windows. Such a detailed study is beyond the scope of this work, and we do not attempt to use timing information further. Nevertheless, it should prove useful in further discriminating between BSM signals from heavy particles and neutrino events.

V. RESULTS
We display the combined results of the analyses of the previous section in Fig. 15. The regions enclosed by solid lines would be probed using our cuts, where we have conservatively assumed that sensitivity will require a minimum signal of 5 total events or 100% of the background after all cuts are applied. We have also indicated, with dashed lines, where 5 events would occur in SBND and MicroBooNE, from Booster protons, and ICARUS, from NuMI protons. These should be thought of as the ultimate possible reach for these analyses, if efficiencies could be improved and backgrounds eliminated. Similarly, we have shown where 5 scalars would decay in ICARUS from kaons decaying at rest in the NuMI absorber, as the maximum achievable sensitivity of a targeted KDAR analysis. For comparison, we also show other current limits on light Higgs portal scalars from CHARM [87,88] (this bound takes into account kaon absorption in the target), LHCb [89,90], E787/E949 [88,91] and SN 1987A [92]. The CHARM and LHCb regions are excluded at the 95% confidence level, while the E787/E949 region is excluded at 90%. The SN 1987A bound is an order of magnitude estimate. Fig. 15 shows that the SBN detectors would be sensitive to new areas of Higgs portal parameter space in the MeV-GeV region, including the gap in the E787/949 kaon decay search due to K → ππ backgrounds. In addition, the Fermilab facilities would afford a better reach than the similar search at the CHARM detector, which used the 400 GeV SPS proton beam. While SPS produced many more B mesons per POT than either of the Fermilab beams we consider, there were only 2.4 × 10 18 POT in all, orders of magnitude less than can be achieved at Fermilab. Also, the thick copper target used for the SPS beam led to many kaons being stopped before decaying, significantly reducing the potential number of scalars from kaon decays [88]. Notably, off-axis production at the SBN detectors from NuMI can lead to even stronger bounds than on-axis production from the BNB. This is again primarily due to the larger number of POT that the NuMI beam can provide.
There are several proposed or upcoming facilities that could also be sensitive to light Higgs portal scalars [4]. Those with relatively similar time scales and bounds to the SBN program include FASER [59,93], NA62 [92,94,95], and SeaQuest [40,96]. In particular, a recent study considering scalars decaying outside NA62 suggests significant improvement is possible for scalars lighter than the dimuon threshold [97]. Nevertheless, for some scalar masses we project that Fermilab facilities can test scalars with mixing angles as low as a few times 10 −5 , in a region which is only otherwise covered by detectors with longer times to construction and data-taking [4]. For m S of order of hundreds of MeV, this represents an order of magnitude improvement relative to the current bounds. Owing to the large number of POT that can be collected at the Booster and NuMI targets, the SBN detectors are well-situated to be the leading probes of the Higgs portal for sub-GeV scalars.

VI. OUTLOOK
Searches for sub-GeV hidden sectors have progressed rapidly in recent years, with implications for models of light dark matter, neutrino masses, and beyond. Intensity frontier experiments provide useful probes of these sectors, and in this work we have investigated the use of the SBN detectors at Fermilab to test the scalar portal. While the detectors are aligned with the 8 GeV Booster proton beam, the 120 GeV NuMI protons can also cause appreciable off-axis production, particularly at ICARUS. The sensitive LArTPC detectors are well-suited to observe the scalar decay products, including electrons, muons and pions.
Using Geant4-based simulations of the beamlines and GENIE to generate neutrino-induced backgrounds, we have analyzed potential signal channels using simple kinematic cuts. As the performance characteristics of the detectors are still under investigation, we have taken generally conservative assumptions in projecting sensitivity curves, which are summarized in Figure 15. Both on-axis and off-axis production of light scalars can be seen at the SBN detectors, which achieve greater sensitivity than existing limits from CHARM, LHCb, and E787/E949. Our results show an improvement of the sensitivity to the scalar-Higgs mixing angle by over an order of magnitude for masses of order of hundreds of MeV, up to the threshold m K − m π . Should our estimated sensitivities be reached, the SBN detectors would provide better probes of the scalar portal than currently planned experiments on similar timescales [4].
With greater understanding of the background, the "ideal" curves in Figure 15 could be attained. One could also consider running in beam dump mode to reduce neutrino background, as was done to search for dark matter in MiniBooNE [52,53]. For m S < 2m e , being able to observe nearly collinear electron-positron pairs while maintaining good separation between electrons and photons would be useful. For m S > 2m e , it will be necessary to improve the discrimination of muons from pions. Being able to reconstruct softer objects can also help our analyses, particularly in the KDAR case where gaps appear at masses corresponding to scalars decaying to slowly moving daughter particles. Finally, we have mentioned the use of timing information to reduce the background by looking for scalars which arrive late in the detectors. It would be useful to perform a more complete analysis of this technique, including out-of-time backgrounds such as cosmic rays, in the future.
Given the small number of renormalizable portals between the SM and a new hidden sector, models for light weakly coupled mediators are highly predictive. The upcoming SBN program at Fermilab will provide a competitive probe of the Higgs portal.

Appendix: Simulation validation
In this appendix, we validate the results of our simulation as described in Sec. IV, by reproducing predictions for neutrino fluxes from both the Booster and NuMI beams. The g4bnb and g4numi codes simulate the production of particles in fixed-target proton collisions, as well as their subsequent interactions and decays. When a neutrino is produced in a decay, information is stored about its production position and momentum, as well as that of the particles in its ancestry. Each neutrino is given an importance weight, which accounts for generator choices (such as a raw weighting towards higher momentum mesons produced in primary collisions for statistical purposes) as well as branching ratios.
For our analysis, we have used the neutrinos produced by g4bnb and g4numi which have kaons as parents. We generate scalars by starting with the kaon position and momentum, and simulating an isotropic 2-body decay, i.e. K → πS. Accounting for the relative branching fractions of K → ν + X and K → πS, we scale the importance weight of each scalar from the neutrino weight given by the Geant 4-based simulations. Then, we check whether each randomly generated scalar would pass through the detector in question, and multiply the weight by an additional factor as in Eq. 6 to include the chance that the scalar decays within the detector.
To test our framework, in Fig. 16 we have used all the neutrinos coming from mesons in the Geant 4-based simulations, i.e. pions and kaons but not muons. Then, we have replaced the K → πS decay with M → πν, in effect reproducing the last step of the Geant 4based neutrino simulations where a neutrino parent is decayed to a neutrino. We have then asked whether each neutrino would pass through various detectors, aiming to reproduce the Booster fluxes predicted at MiniBooNE [71] and MicroBooNE [98], as well as the NuMI fluxes predicted on-axis at MINERvA [72] and off-axis at MicroBooNE [99]. The 3-body muon decays which we neglect contribute negligibly to the fluxes of the neutrino flavors which we use for validation.
In all cases we see very good agreement, as most of the (anti-)muon neutrinos in (anti-)neutrino mode come from two body decays of pions or kaons, and our ignorance of muon decays and three body decays is of little consequence for the fluxes we seek to obtain. Fig. 16 confirms the validity of our procedure for extracting data about mesons from neutrinos as reported by g4bnb and g4numi, as well as for simulating their decays.