First search for dark-trident processes using the MicroBooNE detector

We present a first search for dark-trident scattering in a neutrino beam using a data set corresponding to 7 . 2 × 10 20 protons on target taken with the MicroBooNE detector at Fermilab. Proton interactions in the neutrino target at the Main Injector produce π 0 and η mesons, which could decay into dark-matter (DM) particles mediated via a dark photon A ′ . A convolutional neural network is trained to identify interactions of the DM particles in the liquid-argon time projection chamber (LArTPC) exploiting its image-like reconstruction capability. In the absence of a DM signal, we provide limits at the 90% confidence level on the squared kinematic mixing parameter ε 2 as a function of the dark-photon mass in the range 10 ≤ M A ′ ≤ 400 MeV. The limits cover previously unconstrained parameter space for the production of fermion or scalar DM particles χ for two benchmark models with mass ratios M χ /M A ′ = 0 . 6 and 2 and for dark fine-structure constants 0 . 1 ≤ α D ≤ 1.

We present a first search for dark-trident scattering in a neutrino beam using a data set corresponding to 7.2 × 10 20 protons on target taken with the MicroBooNE detector at Fermilab.Proton interactions in the neutrino target at the Main Injector produce π 0 and η mesons, which could decay into dark-matter (DM) particles mediated via a dark photon A ′ .A convolutional neural network is trained to identify interactions of the DM particles in the liquid-argon time projection chamber (LArTPC) exploiting its image-like reconstruction capability.In the absence of a DM signal, we provide limits at the 90% confidence level on the squared kinematic mixing parameter ε 2 as a function of the dark-photon mass in the range 10 ≤ M A ′ ≤ 400 MeV.The limits cover previously unconstrained parameter space for the production of fermion or scalar DM particles χ for two benchmark models with mass ratios Mχ/M A ′ = 0.6 and 2 and for dark fine-structure constants 0.1 ≤ αD ≤ 1.
A wealth of astronomical data at different scales provide evidence for the existence of dark matter (DM): the motion of galaxies and the stars within them, gravitational lensing, the cosmic microwave background, and the large-scale structure of the universe [1].The nature of dark matter, however, remains elusive.Non-baryonic particles predicted by dark-sector models are candidates for dark matter [2].The search for their production at accelerators is a focus of the high-energy hadron collider program at the LHC [3] and of fixed-target experiments exposed to high-intensity beams [4].
The dark-trident process has been proposed as a new way to search for low-mass dark-matter particles in neutrino beams [5].In this letter, we report a first search for such dark tridents with the MicroBooNE liquid-argon time projection chamber (LArTPC) [6].In the future, similar searches can be performed with the DUNE near detector [7] and the detectors of the Fermilab shortbaseline program [8].
A pair of DM particles, χ χ, is produced through the decay of neutral π 0 or η mesons, which was created by the interactions of the protons and by secondary interactions in the neutrino target (Fig. 1a).The decays π 0 , η → γχ χ are mediated by a virtual, off-shell dark photon A ′ * .The masses of the dark photon, M A ′ , and of the dark fermion (or scalar), M χ , are parameters of the model.The energies of the DM particles χ are typically in the range 0.1-3 GeV for the mass range 10 < M χ < 400 MeV.
The DM particle χ (or χ) then travels uninterrupted to the MicroBooNE detector where it could scatter off argon nuclei through the trident process χ + Ar → χ + Ar + A ′ (see Fig. 1b).The dark photon A ′ promptly decays inside the detector into an e + e − pair.The energies and opening angles of the e + e − pairs depend on the ratio M χ /M A ′ (see Fig. 2).The χ production rate depends on the kinematic mixing parameter ε and a dark finestructure constant α D , which is defined in terms of the dark-photon gauge coupling g D as α D = g 2 D /(4π).We consider the mass ratios M χ /M A ′ = 0.6 and 2 in this search as proposed in Ref. [5].Since M χ /M A ′ > 0.5, the dark photons need to be off-shell to decay into χ χ pairs, and, when on-shell, they will exclusively decay to e + e − .The signal rate therefore scales with ε 4 α 3 D .In this letter, we study dark photon masses in the range 10 < M A ′ < 400 MeV, since the parameter space below M A ′ = 10 MeV is constrained by beam dump searches [9][10][11].The decay η → χ χγ is kinematically forbidden for M A ′ > 457 MeV (assuming M χ /M A ′ = 0.6).Other recent experimental searches cover the mass range where A ′ decays invisibly [12][13][14].
The MicroBooNE's LArTPC has an instrumented volume of 85 tonnes of liquid argon inside a cryostat.Ionization charge drifts across an electric field of 273 V/cm and is read out by one charge collection and two induction planes forming the anode.The LArTPC was simultaneously exposed to the on-axis booster neutrino beam (BNB) [15] and the off-axis beam of neutrinos from the main injector (NuMI) [16].Only NuMI data are used in this search, as the higher average energy of the NuMI beam gives access to higher values of M A ′ .
The NuMI data used for this analysis [17] correspond to 7.2 × 10 20 protons on target (POT), which were taken in two operating modes -forward horn current (FHC) with 2.2 × 10 20 POT (Run 1 from October 2015 to November 2017) and reverse horn current (RHC) with 5.0 × 10 20 POT (Run 3 from November 2017 to July 2018).This data set has previously been used to search for heavy neutral leptons [18,19] and Higgs portal scalars [19,20], and to measure neutrino cross sections [21,22].
We simulate the dark-trident process with a dedicated generator in three steps: the neutral meson flux in the beamline, the decay of the neutral mesons, and the scattering of the DM particles on argon.First, the kinematics of the π 0 and η mesons for both beam configurations, FHC and RHC, are obtained using the g4NuMI simulation [23], which is based on a full GEANT4 description of the beamline geometry.The g4NuMI simulation predicts ≈ 32 π 0 and ≈ 2.5 η mesons per POT, compared to 4.5 π 0 and 0.5 η mesons in Ref. [5], since additional mesons can be produced by secondary interactions within the ≈ 1 m long graphite target and other beamline components.
We then simulate the radiative decays π 0 , η → γχ χ with BdNMC [24].In addition to the scalar DM production supported by BdNMC, we added the option to generate fermions.We calculate the rate of the scattering process χ + Ar → χ + Ar + A ′ inside the LArTPC as a function of the energy of the DM particle and the path traveled inside the detector [25].We compare our signal simulation to the calculations of Ref. [5] and find good agreement in the kinematics, e.g., the distribution of the e + e − opening angle as a function of the energy of each lepton.The cross section of the process shown in Fig. 1b is calculated using GenExLight [25].We find an agreement better than 1% when comparing these cross sections to calculations obtained with MadGraph [26].
We use a "beam-on" data sample to search for the dark-trident signal where the event triggers coincide with the arrival time of neutrinos from the NuMI beam.The background is modeled considering three contributions.Beam-on background events that are triggered by a cosmic ray and not a neutrino interaction are modeled by a "beam-off" sample collected under identical trigger conditions but when no neutrino beam is present.The "beam-off" sample is scaled so that its normalization corresponds to the number of beam spills of the beam-on sample.Neutrino-induced background from the NuMI beam is modeled using a GENIE Monte Carlo simulation [27] embedded in the LArSoft software framework [28].The "in-cryostat ν" sample contains interactions of neutrinos with the argon inside the cryostat, and the "out-of-cryostat ν" sample describes interactions with the material surrounding the detector.
We reconstruct neutrino interactions and cosmic rays within the argon with a chain of pattern-recognition algorithms, implemented using the Pandora software development kit [29,30].The algorithms use hits that are formed from the waveforms read out by the charge col-lection plane and the two induction planes.Collections of hits are reconstructed as a track, as expected for a minimum ionizing particle, or a shower, consistent with being an electron or photon conversion.We use the results of the Pandora reconstruction to select events that are consistent with the signal hypothesis.Dark-trident events are frequently reconstructed as a single shower due to the small opening angle of the e + e − pairs and, in a few cases, as two showers arising from a common vertex.Background processes that can mimic such signal topologies are neutral current (NC) interactions ν +Ar → ν +π 0 /η+X, where the decays π 0 /η → γγ are reconstructed as an e + e − pair.
Each event is therefore required to have at least one vertex, at least one shower, and no tracks.The efficiency of this preselection for a DM signal lies in the range of (32-40)% for masses in the range (10-400) MeV.We find good agreement between the number of data events and the sum of the predictions for the background processes after this preselection (Table I).
We use a convolutional neural network (CNN) for discriminating signal and background based on the previous development of such algorithms in MicroBooNE for multiple particle identification (MPID) [31].Convolutional neural networks (CNNs) are deep learning networks that are ideally suited for images reconstructed from LArTPC data [32][33][34].The CNN architecture is based on a model for dense images with adaptations for LArTPCs.Convolution filters of size 3 × 3 allow scanning of the information contained in showers.The output layer has two neurons that correspond to the probability for signal or background.
We only consider images from the charge collection plane, as it has the best signal to noise ratio [31].Adding information from the induction planes increases the computing time significantly, but has only minimal impact on the performance of the CNN.The size of each image in pixels corresponds to 3456 wires multiplied by 6048 time ticks.To improve processing time, we first compress the time axis by a factor of 6 and then crop the images around the interaction vertex producing a region of interest (ROI) of 512 × 512 pixels.After compression, each pixel has a resolution of ≈ 3 × 3 mm 2 .
We validate the agreement of the vertex reconstruction by comparing data and the background model after the preselection (Fig. 3).The increase of beam-off events towards the top of the detector due to cosmic rays is reproduced by the background model.While we use the reconstructed vertices for the data and background samples, the true vertex location is used for the training.This prevents the CNN from training on an ROI that does not contain the interaction of interest, which can occur when a vertex is reconstructed at a large distance from the true interaction vertex.For the training of the CNN we prepare a dedicated training data set.We use a single signal sample with the parameters α D = 0.1, M A ′ = 50 MeV, and M A ′ /M χ = 0.6.As background samples we use cosmic rays simulated with CORSIKA [35] and ν interactions leading to π 0 mesons simulated with GENIE [27].In addition, we overlay the hits of cosmic rays simulated with CORSIKA to the ν interaction background and the signal samples.Further details on the training of the CNN are provided in the Appendix.
The signal MC samples and the samples listed in Tab.I are passed to the trained CNN model.The events with a CNN signal score < 0 are rejected obtaining a 99% of background rejection efficiency and signal selection efficiencies in the range (27-30)%.Figure 4 shows the signal score distribution obtained for the events contained in the signal region.A data event with a high CNN signal score is shown in Fig. 5, where the shower points in the direction of the NuMI beam.By modifying the training events, we determine that the CNN learns about the orientation of the showers arising from the scattering process.
We evaluate systematic uncertainties for each bin of the CNN score distributions separately for the different signal models and for background.For the in-cryostat ν background, we consider the impact of the neutrino flux simulation (10-20)% [23] and the neutrino-argon crosssection modeling (12-20)% [44], hadron interactions with argon (≈ 1%) [45], and detector modeling (≈ 30%) [46].The beam-off sample is taken from data and therefore has no associated systematic uncertainties other than statistical fluctuations.The impact of the normalization uncertainty on the out-of-cryostat sample and of the POT counting is negligible [19].
The sum of the detector-related systematic uncertainties on the signal is in the range (10-20)%.A form factor accounts for the spatial distribution of the argon nucleus in the χ-Ar scattering [5].Recalculating the cross sections with different form factors [47,48]  The signal rate also depends on the NuMI π 0 and η flux simulated by g4NuMI.We confirm that the ratio of π 0 production relative to π ± production in g4NuMI is consistent with expectations of isospin symmetry.We therefore use the beam flux uncertainty of 22% determined for the charged meson flux [21], which includes hadron production and beam line modeling uncertainties.
The CNN score distributions, shown for one model point in Fig. 4, are all found to be consistent with the background expectation within uncertainties.We therefore proceed to derive limits on the squared mixing parameter ε 2 as a function of M A ′ .The limit setting is performed with the pyhf algorithm [49], which is an implementation of a statistical model to estimate confidence intervals [50].Systematic uncertainties are treated through profile likelihood ratios.The results are validated with the modified frequentist CL s calculation of the COLLIE program [51].
The observed limits of Fig. 6 are shown at the 90% confidence level (CL s = 0.1) for several benchmark points.Since we use a single CNN model for all signal points, the background CNN score distributions are highly correlated between the different mass hypotheses M A ′ .All observed limits are therefore consistently within the 1 and 2 standard deviation ranges around the median expected limit.
In Fig. 7, we compare the results for a scalar dark matter particle χ with existing constraints on dark-trident processes from rare pion decays measured by the NA48/2 collaboration [36], beam dump experiments [9][10][11], and searches for promptly decaying dark photons into e + e − pairs by the BaBar [37], FASER [41], and NA64 [52] collaborations.The limits obtained by the LHCb collaboration [39] apply for higher masses M A ′ > 200 MeV.The most sensitive constraints are obtained for α D = 1 and M χ /M A ′ = 0.6.
For the fermion model, we also compare to reinterpretations of LSND results [5,40].Cosmological constraints on χ χ annihilation in the early universe are obtained using Planck measurements on the cosmic microwave background [42,43].The χ χ annihilation cross section is only relevant for a fermion χ and M χ /M A ′ = 0.6.The cosmological data constrain ε 2 from below, since the thermal relic dark-matter density becomes too small for larger ε 2 [5].
In summary, we apply convolutional neural networks to obtain first constraints on the production of dark matter in a liquid argon detector exposed to a neutrino beam.The dark matter particles are assumed to interact through the dark photon portal [5].We consider thermally-produced fermion and boson dark-matter particles χ with M χ /M A ′ = 0.6 and M χ /M A ′ = 2, and dark fine-structure constants in the range 0.1 ≤ α D ≤ 1.The constraints in the plane of the squared kinematic mixing parameter ε 2 and the dark-photon mass M A ′ exclude , BaBar [37], NA64 [38], and LHCb collaborations [39], and by beam dump experiments [9][10][11] are displayed as shaded regions.The reinterpretations of LSND results [5,40] and the unpublished FASER [41] limits are shown as dashed lines.The two isolated contours at M A ′ ≈ 200-300 MeV are also excluded by LHCb data.The upper limits on ε 2 from Planck data [42,43] apply for fermion DM with Mχ/M A ′ = 0.6.struction, and operation of the MicroBooNE detector as well as the contributions of past collaborators to the development of MicroBooNE analyses without whom this work would not have been possible.We also thank the authors of Ref. [5] for useful discussions about the darktrident model.For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) public copyright license to any Author Accepted Manuscript version arising from this submission.

Appendix
This Appendix discusses the training and performance of the CNN Classifier in more detail.The CNN model is trained during ≈ 10,000 iterations (≈ 5 epochs) with a batch size of 32 images and a learning rate of 0.001 [32].Dropout layers, regularization terms, and batch normalization are implemented during the training to prevent overfitting.The training progress is monitored with a Binary Cross Entropy (BCE) loss function and using the accuracy, which is defined as the fraction of correctly classified images over the total number of images processed by the CNN.
A test set, comprising ≈ 10% of the events included in the training set, is used to evaluate the progress of the CNN training.To determine the number of training steps where the CNN model is frozen we use the receiver operating characteristic curve.In Fig. 8a, the distributions of the signal scores obtained for the frozen CNN model are shown for the training and test sets.The distributions show no indication of overfitting.The confusion matrix, obtained after the CNN training for the test set, is given in Fig. 8b.Approximately 93% of events in the test set are labelled correctly.
The CNN model used in the analysis has been optimized with a benchmark signal point trained against the NC π 0 and cosmic-ray background samples.Before us- Fraction of events MicroBooNE Simulation, Run 3 ing the CNN in the analysis, we validate its performance over different signal points and the full background sample (see Table I).The background rejection obtained for the full background sample is not significantly different from the rejection achieved with the background training sample.We also find that the shape of the CNN score distribution does not vary much for different signal mass points (Fig. 9).

SUPPLEMENTAL MATERIAL Additional exclusion contours and tabulated results
We show exclusion contours for additional model points and give a full list of all results in tabulated form.FIG. 10.The 90% CL limits on ε 2 as a function of M A ′ .The first row corresponds to scalar DM models and the second row to fermion DM models.The constraints provided by the NA48/2 [36], BaBar [37], NA64 [52], and LHCb collaborations [39], and by beam dump experiments [9][10][11] are displayed as shaded regions.The reinterpretations of LSND results [5,40] and the unpublished FASER [41] limits are shown as dashed lines.The two isolated contours at M A ′ ≈ 200-300 MeV are also excluded by LHCb data.The upper limits on ε 2 from Planck data [42,43] apply for fermion DM with Mχ/M A ′ = 0.6.

FIG. 1 .
FIG. 1.(a) A pair of DM particles, χ χ, is produced in a π 0 or η 0 decay; (b) in the dark-trident process, χ (or χ) scatters off an argon nucleus to produce a dark photon A ′ decaying into an e + e − pair with a branching ratio of 1.The rate depends on the kinematic mixing parameter ε and the dark fine-structure constant αD.

FIG. 2 .
FIG. 2. Simulated dark-trident interactions in the Micro-BooNE detector assuming dark photon masses of M A ′ = 50 MeV (top) and 300 MeV (bottom).The horizontal gaps are due to unresponsive wires.

FIG. 3 .
FIG. 3. Distribution of the y coordinate of the reconstructed vertices for the Run 3 data after preselection compared to the background model.The positive direction of the y axis points vertically upwards.The gray band represents the systematic uncertainty in the background model.

5 (
FIG. 4. Comparison of the CNN signal score distribution for Run 3 data with the background model after the preselection.The gray band corresponds to the total systematic uncertainty on the background.The signal distribution for αD = 0.1, M A ′ = 50 MeV, and M A ′ /Mχ = 0.6 is superimposed, scaled by an arbitrary factor. yields uncertainties in the range (2-20)% in the mass range 10 ≤ M A ′ ≤ 200 MeV.

MicroBooNEMicroBooNEFIG. 6 .
FIG.6.The 90% CL observed limits on ε 2 as a function of M A ′ for αD = 0.1 and αD = 1, and (a) Mχ/M A ′ = 0.6 and (b) Mχ/M A ′ = 2, together with the 1 and 2 standard deviation bands around the median expected limits.We use a linear interpolation between the mass points.A total of 13 mass values have been simulated for Mχ/M A ′ = 0.6, equally spaced between 10-100 MeV and between 100-400 MeV and a total of 19 mass values for Mχ/M A ′ = 2, an additional 6 mass values are added at higher MA.A table of the limits at each point is provided as supplementary material.

FIG. 8 .
FIG. 8. (a) CNN score distributions for the training and test samples; (b) confusion matrix, obtained for the test set.

TABLE I .
Numbers of events that remain after preselection normalized to POT for the data and the background model.

TABLE II .
The 90% CL observed limits on ε 2 for a scalar DM particle χ as a function of M A ′ for αD = 0.1 and Mχ/M A ′ = 2, together with the 1 and 2 standard deviation bands around the median expected limits.The results are obtained with a data set of 7.2 × 10 20 POT.

TABLE III .
The 90% CL observed limits on ε 2 for a scalar DM particle χ as a function of M A ′ for αD = 1, and Mχ/M A ′ = 2, together with the 1 and 2 standard deviation bands around the median expected limits.The results are obtained with a data set of 7.2 × 10 20 POT.

TABLE IV .
The 90% CL observed limits on ε 2 for a fermion DM particle χ as a function of M A ′ for αD = 0.1 and Mχ/M A ′ = 2, together with the 1 and 2 standard deviation bands around the median expected limits.The results are obtained with a data set of 7.2 × 10 20 POT.

TABLE VI .
The 90% CL observed limits on ε 2 for a scalar DM particle χ as a function of M A ′ for αD = 0.1 and Mχ/M A ′ = 0.6, together with the 1 and 2 standard deviation bands around the median expected limits.The results are obtained with a data set of 7.2 × 10 20 POT.Scalar Dark Matter, αD = 1.0,Mχ/M A ′ = 0.6

TABLE VII .
The 90% CL observed limits on ε 2 for a scalar DM particle χ as a function of M A ′ for αD = 1 and Mχ/M A ′ = 0.6, together with the 1 and 2 standard deviation bands around the median expected limits.The results are obtained with a data set of 7.2 × 10 20 POT.