Searches for decays of new particles in the DUNE Multi-Purpose near Detector

One proposed component of the upcoming Deep Underground Neutrino Experiment (DUNE) near detector complex is a multi-purpose, magnetized, gaseous argon time projection chamber: the Multi-Purpose Detector (MPD). We explore the new-physics potential of the MPD, focusing on scenarios in which the MPD is significantly more sensitive to new physics than a liquid argon detector, specifically searches for semi-long-lived particles that are produced in/near the beam target and decay in the MPD. The specific physics possibilities studied are searches for dark vector bosons mixing kinetically with the Standard Model hypercharge group, leptophilic vector bosons, dark scalars mixing with the Standard Model Higgs boson, and heavy neutral leptons that mix with the Standard Model neutrinos. We demonstrate that the MPD can extend existing bounds in most of these scenarios. We illustrate how the ability of the MPD to measure the momentum and charge of the final state particles leads to these bounds.


Introduction
The upcoming Deep Underground Neutrino Experiment (DUNE) [1] will use liquid argon time projection chambers (LArTPCs) at its far detector site. The near detector facility is currently under design [2], but is foreseen to include a modular LArTPC [3] (ArgonCube) with a spectrometer immediately downstream, called the Multi-Purpose Detector [4,5] (MPD). The MPD will consist of a high-pressure gaseous argon time projection chamber (HPTPC) surrounded by an electromagnetic calorimeter, all situated within a magnet. The HPTPC has the same nuclear target as ArgonCube and the far detectors but has a much lower density, allowing this detector to more finely resolve the details of low-energy particles produced in neutrino interactions. The LArTPC component of the near detector, although not identical to the far detector, is critical to the DUNE mission of making precision measurements of neutrino oscillations. The MPD's main purpose is to track and momentum-analyze particles exiting the LArTPC, but it will also play an important role in reducing systematic uncertainties for the oscillation analysis. Beyond DUNE's nominal mission to perform precision measurements of neutrino properties and neutrino oscillations, both its near and far detectors will serve as powerful probes of new physics. Recently, interest in using neutrino detectors to search for new physics has grown, with proposals to search for dark matter/dark sectors [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21], millicharged particles [22][23][24], and exotic interactions among the Standard Model (SM) neutrinos [25][26][27][28][29], to name a few (see ref. [30] for a more thorough summary of beyond-the-Standard-Model searches in neutrino experiments). Several of these scenarios have been investigated by the experimental collaborations, leading to stringent limits on the new physics scenarios of interest [24,31]. The next generation of neutrino experiments will undoubtedly be able to interrogate these and other new physics scenarios with increased precision.
In most new physics searches at neutrino facilities, the dominant backgrounds for the signals of interest arise from interactions of neutrinos from the beam with the detector. 1 These scale with the mass of the detector. For new physics searches that rely on particles scattering in the detector, the signal rate is also proportional to the detector mass. However, signal rates in searches for novel particle decays scale as the active volume of the detector. For the DUNE near detectors, the LArTPC and MPD will have roughly the same active volumes, but their different densities imply that their target masses will differ by a factor JHEP02(2020)174 of 50. In searches for decaying particles, the signal-to-background ratio in the MPD will consequently be, roughly, a factor of 50 larger than in the LArTPC.
We propose using the DUNE MPD as a detector for the decays of new particles within its active volume for precisely this reason. Additionally, the excellent tracking ability and charge reconstruction, as well as the particle-identification capability, of the MPD will allow for these searches to be further optimized. For decays of new particles, the MPD, as a standalone detector, will be able to outperform the LArTPC.
As a guiding principle of which new physics searches to study, we focus on the three renormalizable portals to the SM [6,[32][33][34][35][36][37][38][39]: the vector, scalar, and neutrino portals. Each portal consists of a new particle or mediator that couples to SM particles via a renormalizable operator constructed out of Lorentz-and gauge-invariant combinations of SM fields: Here, the new physics particles are represented by F µν , V µ , S, or N while F µν is the U(1) hypercharge field strength, J SM µ is a vector current made out of SM fermion fields, H is the Higgs scalar doublet, and L is the lepton doublet. The parameter µ is a coupling with mass dimension 1. These new particles are interesting in their own right, but they might also be connected to other new particles, such as an extended dark sector. In this manuscript, we will explore the capabilities of the DUNE MPD in searching for the decays of each of these types of new particles.
The remainder of this paper is organized as follows. In section 2, we describe the DUNE Near Detector complex and the MPD. We also provide details regarding the measurement capabilities and expected background rates of interest for the MPD. Section 3 provides details of the simulation techniques we use throughout our studies. Sections 4-7 contain the specific new physics searches for which we advocate: section 4 details a search for dark vector bosons that mix kinetically with the SM hypercharge group; section 5 describes a search for vector bosons that are coupled to various combinations of lepton flavor number; section 6 explains a search for a dark scalar particle that mixes with the SM Higgs boson; and section 7 explores the capability to search for heavy neutral leptons (HNLs) that mix with the SM neutrinos. Also in section 7, we explore the capability of the MPD to determine, if a heavy neutral lepton is discovered, whether lepton number is violated in nature and therefore whether neutrinos are Dirac or Majorana fermions.
In this work, we show that the DUNE MPD will have exceptional capacity to search for new physics in a variety of scenarios. With the DUNE Near Detector design not yet finalized, we use this aspect of the MPD to advocate for its inclusion in the final design for the Near Detector suite. While the LArTPC could search for these new physics scenarios, the added capabilities of the MPD allow for a considerable increase in sensitivity and the possibility of discovering new physics.

DUNE near detector complex
In this section, we outline our assumptions for the setup of the DUNE Near Detector complex, including the proton beam target, decay volume, and the layout of detectors in the experimental hall. We direct the reader to refs. [40][41][42][43] for more detail regarding the DUNE Experiment. Figure 1 provides a schematic of the features of the DUNE target and near detector hall of interest for this work. The figure is not drawn to scale. We consider that the proton beam strikes a target and that the particles emerging from the target enter a magnetic focusing horn system. The focused particles traverse a decay volume; any particles that reach the boundaries of the decay volume are assumed to come to rest and decay. The near detector hall is positioned downstream of the decay volume and the space between them is occupied by raw earth.

Experimental parameters and assumptions
Here we specify the exact experimental parameters that we assume in our simulations, including the beam setup, magnetic focusing horns, and the detectors.
Beam: we assume that the protons in the LBNF-DUNE beam have 120 GeV of energy, and that 1.47 × 10 21 protons strike the target each year. Unless otherwise stated, we assume data collection of equal time in neutrino ("forward horn current") and antineutrino ("reverse horn current") modes. In general, we will assume ten years 2 (1.47 × 10 22 protons on target) of total data collection in our analyses. Different beam configurations, most notably the possibility of the beam protons having 80 GeV of energy or focusing higher energy SM neutrinos, have also been considered by the collaboration. In appendix A, we discuss the impact of different beam configurations on our results, which is marginal. 2 Some projections of DUNE operation include the possibility of an upgraded beam (roughly twice the number of protons on target per year) that can operate in the second half of the experiment. We assume a constant number of protons on target per year, so our ten-year projection is equivalent to a shorter operation time with such an upgraded beam.

JHEP02(2020)174
Focusing horns: consistent with simulations performed by the DUNE collaboration, we assume three aluminum focusing horns that produce toroidal magnetic fields. Depending on the polarity of the applied current, the horns focus either the positively or negatively charged particles that are produced in collisions of the protons with a graphite target embedded within the first horn. As explained in detail in section 3, we include the magnetic focusing effects on the momentum distributions of long-lived, charged mesons -π ± and K ± , specifically. Short-lived (D ± , D ± s ) and neutral mesons (π 0 , η, K 0 L , K 0 S ) produced in our simulations are not deflected by the magnetic focusing horns.
Decay volume and near detector distance: we assume a total distance of 230 m between the target and the end of the decay pipe, wherever relevant. The front of the LArTPC is 574 m from the target, and the front of the MPD is 579 m from the target.
Liquid argon near detector: the liquid argon near detector has a width (transverse to the beam direction) of 7 m, a height of 3 m, and a length (in the beam direction) of 5 m. The fiducial mass of the near detector is 50 t of argon. See ref. [2] for more detail.
Multi-Purpose Detector: the Multi-Purpose Detector (MPD) is a magnetic spectrometer with a cylindrical high-pressure gaseous argon time projection chamber (HPTPC) at its heart. The HPTPC has a diameter of 5 m and a length of 5 m. It is surrounded by an electromagnetic calorimeter (ECAL) and the HPTPC+ECAL system is situated inside a magnet with 0.5 T central field. The axis of the cylinder is perpendicular to the beam direction. However, for simplicity, when simulating our signals (see section 3), we assume the symmetry axis of the detector is aligned with the direction of the beam, since this leads to a uniform apparent detector thickness for all incoming particles. The possibility of including a muon tagger interleaved with the ECAL is under exploration by the collaboration [5]. We discuss how such a muon tagger impacts our results wherever relevant.
The HPTPC operates at a pressure of 10 atmospheres, leading to a total mass of argon in the gas TPC of roughly 1.8 tons (1 ton fiducial). See refs. [4,5,44] for more detail. Parameters associated to the detection of charged particles in the gaseous argon are given in the next subsection.

Particle thresholds and reconstruction capabilities
In this subsection, we discuss various thresholds required for particle identification in the HPTPC, as well as the capability of measuring charge and distinguishing pairs of particles. We make use of the information in ref. [44].
First, we focus on the energy thresholds for which a charged particle would leave a track long enough to be identified. We assume that 2 cm is a reasonable length to identify a track. According to the continuous-slowing-down approximation [45], this corresponds to kinetic energies of 0.14 MeV, 1.8 MeV, and 3.7 MeV for electrons, muons, and protons, respectively. These values assume a pressure of 10 atm in the gaseous argon. Since the energy deposition dE/dx is largely similar for muons and pions, we take their kinetic energy thresholds to be the same. These values are listed in that these are far lower than the thresholds for particle identification in the liquid argon near detector [40]. We assume that transverse momenta are measured at the 0.7% level if they exceed 1 GeV. Below 1 GeV, we assume them to be measured at the 1% level. Charge identification is expected to be efficient and accurate. Two particles with opposite charges emerging from a vertex will be identified as two separate tracks once they are separated by more than 1 cm. For the purposes of this work, this corresponds to nearly perfect separation of positiveparticle/negative-particle pairs, unless the vertex is very close to a detector boundary. Furthermore, there is no ambiguity in the charge assignment. Nearly all experimental signatures studied throughout this work will include at least two charged particles (with opposite signs) emerging from a common vertex. Events with only one charged particle, assuming that particle has a very high energy, could be more difficult to charge-identify.

Particle signatures in the multi-purpose detector
In the HPTPC, energy deposited in the form of ionization is sampled by a plane of chargesensitive pads at each end of the cylindrical detector. The characteristics of the charge (energy) deposited per unit length (dE/dx) for each track enable this detector to distinguish electrons from heavier particles for a large range of particle momenta, and some ability to distinguish among the tracks left by heavier particles (pions, muons, kaons, protons), especially if they come to a stop within the gas volume. The ECAL design is currently undergoing optimization, but in its current state it is only one interaction length thick. By itself, this ECAL should have some ability to separate interacting pions from (1) muons and (2) pions that punch through (all muons and roughly 30% of pions will punch through and escape the ECAL). With the addition of a muon tagger (currently under consideration [5]) outside the ECAL, more separation power would be achieved. Table 2 lists the types of signals that are expected for the new physics processes discussed in later sections, along with links to those sections. In general, the signals are oppositely charged pairs of leptons or hadrons whose reconstructed tracks should point back to the target. The reconstructed invariant mass will provide another handle for identifying some of the possible decay channels. Finally, the decaying particles are produced associated with the high-energy proton-beam-target interaction and will have O(10) GeV of energy, leading to decay products with several GeV each. In contrast, much of the background (produced in neutrino interactions) will have lower energy, giving another kinematical variable with which signal events can be distinguished. X → e ± π ∓ /µ ± π ∓ Invariant Mass, Direction Heavy Neutral Leptons
In the HPTPC, possible background pairs of e + e − arise from photon conversions. The conversion distance 3 of photons in the gas is O(10 meters); approximately 12% of photons produced by standard beam neutrino interactions in the gas will convert before reaching the ECAL. However, these converted e + e − pairs will typically be of much lower energy than the signal pairs considered in this paper and they will infrequently correlate with the beam direction. We discuss backgrounds of this nature in detail in section 4.
For momenta above a few hundred MeV/c, it is not possible to distinguish muons from pions via dE/dx alone in the HPTPC. With the addition of the ECAL (and potentially a muon tagger [5]), some fraction of the pions will be distinguishable by their hadronic interactions. We discuss the capability of distinguishing muons and pions in the relevant sections of this paper. Backgrounds could then arise from standard charged-current neutrino interactions in the gas, in which a muon and a pion of opposite charges are produced with no other visible particles, but these would still need to point back in the direction of the target to be considered background to signals discussed in this paper.
Neutral pions are only visible in the HPTPC if the photons from the pion decay convert. Most of the time, the photons will reach the ECAL, where they will convert and produce signals with a shower axis that points back to the π 0 decay point. Given the O(meters)-long photon conversion length, rarely (about 1% of the time) will both photons convert within the gas. If one photon converts in the gas and the other reaches the ECAL before converting, then the energy depositions in the two subdetectors will still point back to a common vertex, allowing reconstruction of the π 0 invariant mass. Backgrounds to the new physics processes are expected to be small, since the background processes would need to produce multiple pions and have no other activity at the neutrino interaction vertex. Multiple pion production in standard neutrino interactions is generally accompanied by additional vertex activity, which is detectable down to very low energies in the HPTPC. Table 3 lists the expected number of neutrino-related events for scattering on argon for a variety of different interaction types [44,46] Table 3. Expected number of neutrino-related events scattering on argon (liquid or gas) assuming 1 t of fiducial mass and 1 year of data collection.

Background rates from the neutrino beam
work, these events can contribute as background (assuming they have the right kinematics or could otherwise be misidentified as our desired signal). We will use these estimated rates to infer the total number of background events for each search.
As an example, consider the total number of ν µ charged-current events in one year of operation -1.64×10 6 -assuming one ton of fiducial mass. Charged-current ν µ -scattering on argon will lead to the production of a single charged pion along with the muon (and no other detectable hadronic activity at the interaction vertex) O(10%) of the time. Therefore, we naively expect O(10 5 ) of this type of background event -a µ − π + pair -per year of data collection in the MPD. If, for example, we were searching for a signal that consists of µ + µ − pairs, then this background could be further reduced using the kinematic information of the charged particles, information regarding the hadronic system, particle identification capabilities, etc. We will discuss these techniques in detail in sections 4.3 and 7.3.1.

Simulation details
All of the new physics scenarios studied in this work rely on simulating the flux of SM particles that then decay into beyond-the-Standard-Model (BSM) particles at or near the DUNE target. We will be interested in the flux of the BSM particles, both in terms of the direction and energy spectrum, in order to determine the probability for a given BSM particle to reach the DUNE MPD and decay inside it. In this section, we discuss the details of simulating the production of SM particles (predominantly mesons) that decay into BSM particles.

Production of charged and neutral mesons
The general strategy of DUNE, while operating as a neutrino beam experiment, is to produce a large flux of charged mesons (mostly π ± and K ± ) that decay leptonically, leading to a large flux of SM neutrinos (mostly ν µ or ν µ , depending on the mode of operation) at the near and far detectors. In order to produce a purer beam (i.e., one with a smaller contamination of ν µ with respect to ν µ when operating in neutrino mode, or vice versa for antineutrino mode), focusing horns are utilized immediately downstream of the target, as depicted in figure 1. This allows for focusing of positive (forward horn current) or negative (reverse horn current) mesons. Forward (reverse) horn current corresponds to an JHEP02(2020)174 Mesons/POT 2.7 2.4 0.24 0.16 3.7 × 10 −6 6.0 × 10 −6 1.2 × 10 −6 1.6 × 10 −6 Table 4. Average number of charged mesons produced per proton-on-target assuming 120 GeV protons.
enhancement of the (anti)neutrino flux, so we will refer to the modes as "neutrino" and "antineutrino" modes henceforth. We will explore scenarios in which the flux of a new particle is produced via decays of charged or neutral mesons. To simulate the production rate of a given meson, we use the software Pythia8 [47], simulating 120 GeV protons on target and using the ''SoftQCD:all'' option to generate mesons. In doing so, we obtain the overall rate of meson production from the primary proton interaction; this procedure underestimates the total meson flux as it neglects secondary production when mesons produced in the primary interaction scatter in or near the target. Our results may then be interpreted as moderately conservative.
The average numbers of different charged and neutral mesons produced per proton on target (POT), assuming a 120 GeV beam, are listed in tables 4 and 5, respectively. Pythia8 also allows us to obtain the lab-frame four-momenta of the outgoing mesons, which will be of interest in simulating the flux of different new particles. Appendix A provides more details of this simulation as well as comparisons between 80 GeV and 120 GeV proton beam results.

Focusing of charged mesons
The quantities we explicitly require in order to simulate meson decays are not their fourmomenta at production but instead their four-momenta at the time of decay. Neutral mesons (π 0 , K 0 L , K 0 S ) and charged mesons that decay promptly (D ± , D ± s ) are unaffected by the focusing horns, so, in our simulations, we use the lab-frame four-momenta obtained with Pythia8 as the four-momenta at decay, assuming they don't hit the boundaries of the decay volume. Long-lived charged mesons (namely, π ± and K ± ), on the other hand, are subject to the effects of the focusing horns. For the distributions of the long-lived charged mesons, we use output from the DUNE Beam Interface Working Group (BIWG) [46], which makes use of Geant4 [48,49] and Fluka [50,51]. Figure 2 depicts the normalized differential decay rates of the long-lived charged mesons extracted from Geant4 and Fluka (with operation in neutrino mode), as a function of the longitudinal distance z decay in the left panel and  Figure 2. Decay distributions of long-lived charged mesons as determined by the DUNE collaboration [46] using Geant4 and Fluka simulations of charged meson transport through the DUNE focusing horns. The left panel presents decay distributions in terms of z decay , the longitudinal position of the decay. The right panel presents decay distributions in terms of the ratio of transverse and longitudinal momentum, |p T /p z |, as a proxy for how well-focused the corresponding mesons are. In both panels, we show π + distributions (blue), π − distributions (red), K + distributions (green), and K − distributions (purple). This figure is for operation in neutrino mode, in which positively charged mesons are focused.
in the right panel. The left panel of figure 2 illustrates several effects of the focusing horns. First, we see that the positively charged mesons tend to decay at larger z -this is because they are the focused mesons, and travel through the horns and beam pipe to larger distances. Secondly, we also see a peak of decays for all four distributions located near z decay ≈ 230 m, corresponding to the end of the decay pipe. Mesons that reach the end of the decay pipe will stop and decay at rest. In the right panel, | p T pz | characterizes how well-focused a given meson is at the time of its decay -a perfectly focused meson would have p T = 0 and point towards the near detector. We clearly see two effects here -first, the distribution of the positively charged mesons peaks at smaller | p T pz | than that of the negatively charged ones, as the former are focused. Also, we observe that there is a greater separation, due to focusing, between the π + vs. π − distributions than for the K + vs K − distributions: pions are focused more efficiently than kaons.
In contrast to the right panel of figure 2, figure 3 depicts the normalized differential decay rates as a function of | p T pz | for a subset of the unfocused species -namely, π 0 , η, K 0 L , and K 0 S . The distributions for promptly decaying D ± and D ± s are similar to those in figure 3 and we choose not to display them. Comparing figure 3 with the right panel of figure 2, the bulk of decays of the unfocused mesons lie near | p T pz | ∼ 3 × 10 −2 , while those of the focused mesons lie near | p T pz | ∼ 10 −2 . Instead, the bulk of decays of the anti-focused, wrong-sign mesons (specifically π − in this case) lie near | p T pz | ∼ 10 −1 . As discussed, we simulate the production of neutral mesons using Pythia8 and determine the position of their decay according to the proper lifetime of the meson and its simulated energy. We insist that any meson that reaches the end of the decay pipe decays at rest at the rock surface, see figure 2 (left). Figure 3. Identical to the right panel of figure 2, except for unfocused neutral mesons: π 0 (blue), η (red), K 0 L (green), and K 0 S (purple). We include the focused/unfocused π + /π − distributions for comparison in grey solid/dashed lines.

Particle flux at the DUNE near detector
Once the parent distribution has been obtained, we must determine the flux of the BSM particle(s) of interest at the DUNE near detector. This is done on a particle-by-particle basis by Monte-Carlo simulating the rest-frame decays of the parent particle to obtain the four-momenta of the relevant daughter(s). These four-momenta are boosted from the parent's rest frame back to the lab frame. Combined with the position (x, y, z) of the meson decay location, we can determine whether the BSM particle will pass through the MPD. The fraction of produced particles that pass through the near detector defines the geometrical acceptance of the detector. A careful treatment of this effect is critical for the DUNE near detector. Since it has a small solid angle [O(10 −4 rad 2 )], the flux of new particles is highly dependent on features such as the boost of the parent meson. Because of our assumption that the detector is coaxial with the beam direction (see section 2), we have (effectively) reduced the solid angle by roughly 20%. This lower acceptance will be compensated in our signal event rates by the longer effective depth of the detector; these two effects will approximately cancel.
We are interested in decays of parent mesons into both two-and three-body final states. The two-body kinematics may be solved easily in the parent rest frame; however, the three-body kinematics is less straightforward. We discuss our recipe for handling threebody decays in the next subsection.

Three-body decays
In some cases, we are interested in the production of a particle that comes from a threebody decay. Here we discuss the assumptions made in these situations and their validity. Consider the three-body decay of a particle P into particles S 1 , S 2 , and X (with masses m P , m S 1 , m S 2 , and m X , respectively); we will be interested in the lab-frame kinematics of X (i.e., its energy and whether it is pointing in the direction of the detector). In all of the scenarios we consider in this work, P is considered to be spin-0, so the angular distribution of the outgoing particles is isotropic in its rest-frame. In contrast to two-body decays, however, X is not monoenergetic in that frame. Kinematics dictates that the rest-frame energy of X must lie in the range (3.2) Beyond this, different interaction structures governing the decay P → S 1 S 2 X determine the shape of the differential width dΓ/dE X .
Our goal in this subsection is to investigate the dependence of our analysis on the choice of dΓ/dE X , and to show that the dependence is fairly small. As an illustrative example, we focus on the decay K + → µ + π 0 N , where N is a heavy neutral lepton (see section 7) with a mass of 200 MeV.
We simulate the three-body decay (isotropic in the K + rest frame) according to three different energy distributions: (a) flat between the energy endpoints (200 MeV to 225 MeV) (b) piecewise linear, symmetric about the middle of the energy range, and (c) the true energy distribution for HNL production, as described in ref. [52]. These three energy distribution probability density functions are depicted in the inset of figure 4 (right). The decays are simulated using the focused K + distributions (for simplicity, we concentrate on operation during neutrino mode only), and we determine (1) the spatial distribution of N JHEP02(2020)174 at the location of the MPD, and (2) the energy distribution of the N particles that pass through the MPD. Figure 4 depicts the results of the simulations. The left panel depicts the radial distribution at z = 579 m (the front face of the MPD), as measured from the center of the beam axis. 4 The right panel depicts the energy distribution of accepted N (i.e., those with R (at det.) < 2.5 m). In all panels of figure 4, the blue curves correspond to the flat energy distribution, the red correspond to the piecewise-linear distribution, and the green correspond to the true energy distribution. In the left panel of figure 4, we see that the flat and true distributions both tend to agree in terms of radial distributions, so we expect that the acceptance fraction in simulating these two scenarios should be largely similar. The piecewise-linear distribution tends to prefer larger radii. Similarly, in the right panel, we see that the energy distributions of accepted events for the flat and true distributions match fairly well, where the piecewise-linear distribution prefers higher energies. Based on this evidence, when simulating three-body decays, we assume the distributions of the decay products are flat in the phase space of the final-state particle energy. For our ambitions here, this simplified analysis suffices, especially when it comes to estimating detector sensitivities.

Signals of decays in the Multi-Purpose Detector
We are, in general, interested in the event rates of different decays of new particles in the DUNE MPD. After determining the flux of a given new particle at the MPD location and determining the associated probability for a decay of interest to occur, we also simulate the kinematics of the final-state particles. All of these are model dependent and discussed in detail in sections 4-7.
Of interest are quantities such as the energy of visible particles and opening angles between different particles that leave a track in the MPD. When any threshold or efficiency effects (see table 1) are relevant, we apply them. We will discuss how the kinematic quantities will allow for an improved search for these new particles, where relevant. For instance, any new particle that decays completely visibly, such as a dark photon decaying to an electron-positron pair, will produce decay products with an invariant mass peak once the total four-momentum of the final-state particles is reconstructed. This allows for good separation of signal from background since the latter is not expected to have any such feature.
We will also be specifically focused on the decay products' kinematics when discussing disentangling different new physics scenarios. In particular, when discussing whether a newly-discovered HNL is a Dirac or Majorana fermion (see section 7.4), we will focus on the lab-frame decay product energies as a proxy for the new particle's rest-frame decay distribution. Such kinematic information can provide an interesting handle into whether lepton number is violated, and we will exploit it when possible.

JHEP02(2020)174 4 Dark photon
In this section and section 5, we explore the possibility that a new vector boson exists in nature. Here, we will focus on the case in which this vector boson acquires small couplings to SM fermions only through kinetic mixing with the SM U(1) hypercharge group -the so-called dark photon -interacting via the vector portal, where M A is the mass of the dark photon A with field strength F µν and ε characterizes the strength of kinetic mixing with the SM hypercharge gauge-group, with field strength F µν . The kinetic mixing term is a renormalizable operator and the size of ε depends upon how it is generated [54]; we will be interested in ε < 10 −4 . As discussed in section 1, many new physics search strategies utilizing neutrino detectors have been proposed, specifically in searches for dark matter and associated dark sectors. One particular intensely studied new physics model, due to its simplicity, is one in which dark matter is charged under U(1) and couples to the SM via the dark photon introduced in eq. (4.1). Such dark matter can be fermionic or scalar, and current and future neutrino experiments are capable of probing well-motivated regions of parameter space in which the dark matter is a thermal relic, symmetric between dark matter particles and antiparticles in the early universe. The search strategy adopted in these proposals typically calls for the dark matter to be produced in meson decays either via an off-shell dark photon, or a dark photon that decays promptly into dark matter pairs. The dark matter then travels to the detector and scatters, depositing energy that can be measured at the neutrino near detector (typically scattering off nuclei or electrons, although more exotic signatures have been proposed [17,55]). In this section, we focus on searching only for the dark photon without requiring the existence of dark matter, χ. Should the dark photon be part of a larger sector containing dark matter, we require that decays of the dark photon, A → 2χ, are forbidden, i.e., m χ < M A < 2m χ . Thus, any dark photon that is produced on-shell will only decay into SM particles. Depending on the strength of the kinetic mixing ε, such decays may lead to long-lived dark photons. This is the region of parameter space in which we are interested for this study.
Previous experiments -among them, E141 [56][57][58][59], Orsay [60], NuCal [20,61,62], and E137 [57] -have probed a similar region of parameter space in this fashion. 5 The main differences between the various experiments is the production mechanism and the distance between the target where the dark photons are produced and the instrumented detector region. The interplay between these means each experiment has a sweet spot where it will be best suited to search for a particular combination of dark photon mass and kinetic mixing parameter. For the DUNE MPD, we will be interested in A produced via the decays of light mesons π 0 and η into γA , as well as production from the continuum process pp → ppA . Dark photons produced in these ways have energies of several GeV. The DUNE MPD will be sensitive to regions of parameter space for which the lab-frame 5 Ref. [63] provides a thorough review of searches for dark photons and existing constraints.

JHEP02(2020)174
decay length of the A is O(100 m) or longer, where decays at the location of the MPD will be optimal. In this section, we discuss the production and decay mechanisms of dark photons in the parameter space of interest for the DUNE MPD. We discuss the associated backgrounds for searches of this type in section 4.3, and provide our estimates for the DUNE MPD sensitivity to these dark photons in section 4.4.

Dark photon production
For the DUNE MPD, the dominant production mechanism for dark photons is the decay of neutral mesons, specifically π 0 and η. We simulate the production of these as described in section 3. The branching ratio of a given neutral meson m into γA is From ref. [64], Br(π 0 → γγ) = 98.823% and Br(η → γγ) = 39.41%. Given the spectra of π 0 and η produced in the DUNE target and the kinematics of m → γA decay, we estimate (a) the fraction of produced A that reach the HPTPC detector and (b) the energy spectrum of the A flux that reaches the detector, as described in section 3.
In addition to meson decays into A , the bremsstrahlung process pp → ppA allows us to probe larger values of M A , here limited by √ s 2m p E p = 12 GeV. The calculation of the production in this process is described in refs. [14,62]. The production rate is expressed in terms of properties of the outgoing A , with momentum where p T is the transverse momentum of the outgoing A , E A is its energy, P is the initial proton incoming momentum, z is the fraction of the proton's initial momentum transferred to the longitudinal momentum of the A , and φ is an azimuthal angle. The double-differential production rate is where σ pN (s) is the cross section of a proton hitting a target nucleus N at center-of-mass energy squared s, and w ba (z, p 2 T ) is the photon splitting function, and and m p is the proton mass. In eq. (4.3), the ratio of total cross sections, evaluated at the two energies of interest, s = 2m p (E p − E A ) and s = 2m p E p , is close to one for the energies and masses of interest here. The form factor, F 1,N , in eq. (4.3) encodes the mixing of the A with SM vector mesons ρ and ω, leading to an enhancement in dark photon production when M A ≈ 800 MeV [65].
To determine the number of A that are produced in the direction of the detector, one must integrate eq.
Kinetically−Mixed Dark Photon Figure 5. Expected number of A that travel towards the DUNE Multi-Purpose Near Detector, a distance of 579 m from the production target, as a function of mass M A . We assume ten years of beam operation, ε = 1, and a beam luminosity of 1.47 × 10 21 protons per year. The red curve displays A produced via decays of π 0 , the blue line displays those from η decays, and the green displays those produced in the bremsstrahlung process pp → ppA .
geometry. The majority of the existing literature [14,66,67] focuses on detectors with a much larger solid angle (as viewed from the beam target) than we consider here. The traditional approach is to integrate over a specific range of z and |p T |, e.g., z ∈ [0.1, 0.9] and |p T | < 1 GeV for SHiP [14]. In contrast, here, we transform the variables of eq.
In these variables, the requirement that the A enters the detector leads to the constraint where θ det. is half the angular size of the detector, as viewed by the target. For the DUNE MPD, θ det. ≈ 4 × 10 −3 . Note that this relationship holds regardless of the detector's solid angle, and we advocate for its use whenever considering dark photon production in this process. As explained in refs. [62,68], the form of eq. (4.4) is only valid in certain limits, requiring Some of these restrictions are automatic once eq. (4.5) is satisfied. We include these restrictions in our calculation and only include contributions to the flux from regions in which eqs. (4.5) and (4.6) are satisfied and thus eqs. (4.3) and (4.4) are valid. Figure 5 displays the expected number of A that traverse the HPTPC in 10 years of operation at DUNE, normalized to ε = 1; this flux scales with ε 2 . As explained above, bremsstrahlung production allows us to reach larger M A . The peak in the bremsstrahlung curve for M A 800 MeV comes from enhanced mixing near the ρ meson mass. In figure 5, we see that the bremsstrahlung process (green) is comparable, even at M A 400 MeV, to production via η meson decay. This is in contrast with projections for this process at, for instance, SHiP [66] or SeaQuest [67], where the prediction is that the A flux due to η decays is an order of magnitude or so larger than that from bremsstrahlung. This difference JHEP02(2020)174 is due to the fact that the DUNE near detector is much further from its target than SHiP or SeaQuest, and has a significantly smaller solid angle than either. This leads to a relatively larger geometrical acceptance for the bremsstrahlung process relative to meson decays for DUNE, resulting in a comparable A flux, as seen in figure 5.
Direct production of dark photons via qq → A may contribute for high center-ofmass energy proton collisions. However, such a calculation requires robust knowledge of the proton parton density functions for s = M 2 A . This mechanism could increase sensitivity, but only for M A 2 GeV -we do not include this in our calculations, and expect that it would not affect the sensitivity we present here.

Dark photon decay channels and lifetime
We consider three decay modes of A : A → e + e − , A → µ + µ − , and A → hadrons. For the A masses of interest, the decay channel A → hadrons consists largely of the final state π + π − . However, depending on M A , more complicated final states exist (see, e.g., ref. [69], for a discussion of the different decay channels for M A 1 GeV). We reproduce the partial widths of A here for clarity. For decays into lepton pairs + − , the partial width is In order to express the partial width of A into hadron pairs, we rely on the R ratio, Combining these expressions, we display the branching fraction of A into each final state as a function of M A in figure 6.
Since we are interested in cτ ≈ 100 m, we can use these widths to determine the expected parameter range for which DUNE's sensitivity is maximized. We find, roughly, that the event rate peaks for ε 2 ≈ 10 −15 (1 GeV/M A ). Other factors, especially the boost of A , modify this relationship slightly.

Backgrounds
In this section we discuss the expected background rates for the different search channels A → e + e − , A → µ + µ − , and A → hadrons. All of the backgrounds we discuss here come from neutrinos in the DUNE beam scattering off an argon nucleus in the HPTPC. 6 The signal comes from a decaying particle (that originated in or near the target), compared to JHEP02(2020)174 a light neutrino scattering off a heavy nucleus, so we can take advantage of the fact that the outgoing charged particles point back to the target to greatly reduce background rates.
A → e + e − : the dominant background for this channel is ν + Ar → ν + Ar + π 0 , i.e., neutral-current single-pion (NCπ 0 ) events. This also includes π 0 produced in coherent production, where the recoiling target is too low-energy to be observable. Events of this class may be misinterpreted as signals with only e + e − pairs if the following occur: • No other hadronic activity occurs at the vertex of the neutrino interaction.
• Of the two photons coming from π 0 decay, one is either missed by the ECAL, or is, for some reason, not attributed to a common vertex with the second photon.
• The other photon converts in the gas (which occurs for roughly 12% of photons).
• The converted photon's electromagnetic shower is identified as an e + e − pair, with a direction consistent with the direction from the target/decay pipe.
In all, O(10 5 ) NCπ 0 events are expected per ton-year of exposure. We expect that the neutrino-related background of O(10 6 ) events (for the exposure we consider) may be drastically reduced by the requirements listed above. Specifically, we expect that 80% of NCπ 0 events will have detectable hadronic activity (at least one charged particle with kinetic energy above 5 MeV) in the MPD. This means that only 20%, or 2 × 10 5 events in ten years, will be NCπ 0 with no significant hadronic activity. As for the second requirement, we estimate (conservatively) that the DUNE MPD will miss 10% of photons coming from π 0 → γγ decay. This reduces our background sample to 2 × 10 4 events. The third requirement, that the other photon converts in the gas TPC (before reaching the ECAL) will occur in 12% of events -this brings our background down to 2400 events. Finally, we estimate the probability that a photon from π 0 decay will be pointing in the direction of the target. Because most of the π 0 are coming from nuclear emission, they will be emitted JHEP02(2020)174 nearly isotropically. If we require that the angle of the photon is within 4 mrad from the beam direction (the angular resolution of the gas TPC), then roughly 0.01% of photons will pass this cut. A more conservative cut of 1 • (5 • ) increases this to 0.1% (2.6%). Given all of these considerations, we expect an optimistic background rate of < 1 event in ten years. Even with a very conservative 5 • angular cut, this is a background rate of O(50) events. If such a conservative cut is required (perhaps to retain a signal from decays downstream of the beam target, such as that discussed in section 7), then we still have additional handles to separate signal from background. For the decay A → e + e − , the electron/positron pair will have an invariant mass consistent with M A , where the electron/positron pair from π 0 → γγ, γ conversion to e + e − will have zero invariant mass. Additionally, these background electrons/positrons will tend to have lower (around 500 MeV) energy, where the signal electrons/positrons will have several GeV of energy. These kinematical separations may be exploited if a 1 • angular cut turns out to be too optimistic.
In addition to π 0 faking the A → e + e − signal, there will be an irreducible background of events with e + e − being the only visible particles in the final state: neutrino trident events. Refs. [70,71] estimate O(20) neutrino trident e + e − events in the MPD (we obtain this by scaling their LArTPC projections to the mass of the HPTPC) in ten years. Their kinematics should be different from our signature (in terms of invariant mass, direction, etc.) and should not contribute sizably to the search proposed here.
A → µ + µ − and A → hadrons: we consider these two channels in the same category as they share the same background, predominantly ν µ charged-current-single-pion (CC1π) events, with an expected rate of O(10 5 ) events per ton-year. The authors of ref. [70] explored a similar set of signals and backgrounds, focusing on the search for µ + µ − neutrino trident events in the DUNE LArTPC near detector (where the background rates are larger by a ratio of the fiducial mass in each detector). They discussed certain background rejection techniques, including the rejection of events with hadronic activity 7 (around 10% of background events survive this cut), and a variety of kinematical cuts that would apply to our search as well (roughly 1% of background events survive). We also note that the ECAL surrounding the HPTPC will have moderate ability to identify charged pions vs. muonsapproximately 70% of charged pions will interact hadronically in the ECAL, allowing for a positive identification of those that do interact. Additionally, the DUNE collaboration is considering adding a muon tagger [5] (likely consisting of alternating layers of scintillator and a high-density material like steel) that will further reduce pion/muon misidentification.
Finally, we note that the same angular cuts discussed regarding A → e + e − can be applied here -when a muon and a pion fake this signature, the direction of the pair will very rarely be in the direction of the beam. This allows for even further background reduction than when trying to reduce the backgrounds to neutrino trident events [70]. Overall, we expect that ten signal events will constitute a significant excess over background, especially after considering that our signal events will share a common A invariant mass, where the background events will be smoothly distributed in that variable.

Dark photon sensitivity
Given the energy spectrum of the A that reach the DUNE MPD, we calculate the expected number of A decays as a function of M A and ε 2 for a given amount of data collection. In figure 7, we display the regions where we expect at least 10 (blue) and 100 (red) such decays, assuming ten years of data collection. In addition, we show the excluded region of this parameter space constrained by E141, Orsay, NuCal, and E137 [56][57][58][59][60][61][62] and from emission during Supernova 1987A [72] in shaded gray. These regions correspond to a total number of decays being 10 or 100 -in order to determine the number of a certain type of decay, one must include the branching fraction into e + e − , µ + µ − , or hadrons shown in figure 6. For instance, if M A ≈ 300 MeV and ε 2 ≈ 10 −15 , a total of 10 decays are expected in ten years. According to the branching fraction of A , these events should be roughly 40% µ + µ − decays and 60% e + e − . Searching for the correct combination of different final states can also improve the capabilities of a search for dark photons, and disentangle signals of dark photons from those of the other new physics scenarios we discuss in later sections. In producing figure 7, we assume 100% efficiency in identifying particle/antiparticle pairs -this approximation likely breaks down for lighter A , but in those regimes, existing limits are more powerful than the DUNE MPD search for dark photons.

JHEP02(2020)174
Overall, we find that the DUNE MPD will be able to extend sensitivity to a dark photon in the sub-GeV mass regime for kinetic mixing in the ε 2 ≈ 10 −16 − 10 −14 range. This complements other future experiments, such as FASER [73,74], SHiP [66], and SeaQuest [67], which will be sensitive to larger values of ε 2 , as the detectors associated with these proposals are significantly closer to their production targets. If dark photons exist and are produced in this type of environment, then those produced at CERN (for detection at FASER/SHiP) would be much more boosted than those produced at DUNE. These detectors are then sensitive to dark photons with relatively shorter proper lifetimes, meaning that DUNE is sensitive to smaller ε 2 for similar mass dark photons.

Leptophilic gauge bosons
While the dark photon couples to the Standard Model purely via kinetic mixing, as discussed in section 4, new vector bosons may also interact with the SM if they are gauge bosons of a new group under which some SM particles are charged. Such a gauge symmetry should be anomaly free, which often requires extending the field content of the SM. However, if SM fields are charged in a vectorlike fashion, for instance differences of baryon/lepton number between two generations, then no additional fermions are necessary. We expect that DUNE will have the strongest sensitivity to such new physics when the new symmetry is leptophilic, and therefore focus on the case of a L α − L β gauge boson with α, β = e, µ, τ [75][76][77][78][79][80][81][82][83][84][85][86][87][88].
The leptophilic gauge boson 8 V of mass M V , associated with gauging the difference between lepton numbers L α − L β , will couple directly to neutrinos or charged leptons in generations α and β. This coupling means that the gauge boson will be produced in charged meson decays through final state bremsstrahlung and can decay to neutrinos and charged leptons, if kinematically accessible. In addition to this direct coupling, kinetic mixing between V and the SM photon is induced at one loop via SM particles that are charged under both SM hypercharge and the new U(1). This loop-induced kinetic mixing, as defined in eq. (4.1), is related to the L α − L β gauge coupling g αβ as [63,83,89] where m α is the mass of charged lepton with flavor α and q µ is the four-momentum of the boson V that is mixed with the SM photon. This contribution is in addition to any kinetic mixing that might be present in the UV Lagrangian. From now on, we will assume that the bare kinetic mixing is zero and that the effective mixing is determined by eq. (5.1). This kinetic mixing allows for V production through neutral meson decays and for V to decay to a pair of charged leptons that are not in either generation α nor β. Hence, even L µ − L τ gauge bosons lighter than 2m µ will have a visible decay channel. We are interested in situations in which a V is emitted on-shell due to this mixing, so q 2 ≡ M 2 V . In the limit where the momentum transfer (i.e., the V mass) is well below the JHEP02(2020)174 lighter of the two charged lepton masses, the kinetic mixing becomes a constant, In the opposite limit, where the momentum transfer (or M V ) is much larger than both of the charged lepton masses, We show the kinetic mixing for each leptophilic U(1) symmetry in figure 8. Assuming V is produced on-shell, we show the kinetic mixing squared relative to the new gauge coupling squared for L e − L µ (red), L e − L τ (blue), and L µ − L τ (green). The change of behavior of the loop integral at particle thresholds is apparent. With any kinetic mixing, the production processes discussed in section 4 will apply here, suppressed according to eq. (5.1) and figure 8. In addition, V may be emitted in decays involving charged leptons and neutrinos by being radiated in the final state; we provide a derivation of the branching ratio for such decays in appendix B. Which of these processes is most relevant depends on how strong the kinetic mixing is, as well as the lepton flavor(s) of interest. We discuss the production mechanisms in the following, and then the decay signatures in section 5.2. We provide our sensitivity estimates in section 5.3.

Production of leptophilic gauge bosons
In addition to the production via kinetic mixing, V may be radiated in charged meson decays with final-state leptons. The greatest abundance of these will be from charged pion JHEP02(2020)174  Figure 9. The number of L e − L µ gauge bosons V traveling towards the DUNE MPD as a function of mass M V assuming ten years of data collection, and equal run time in neutrino and antineutrino modes. Solid lines display expected number due to charged meson decays -purple (π ± ) and orange (K ± ). Dot-dashed lines display expected number due to processes relying on kinetic mixing -red (π 0 decay), blue (η decay), and green (pp → ppV bremsstrahlung). and kaon decays, π ± /K ± → ± α ν α V . In appendix B, we derive the branching ratios of these three-body decays. The rate for these three-body decays is proportional to the square of the final-state lepton mass [90,91].
For a L e − L µ gauge boson, the processes π ± → e ± ν e V and π ± → µ ± ν µ V (and the equivalent kaon decays) will contribute to V production; the decays involving electrons, however, are helicity suppressed relative to those involving muons. Production via π ± → e ± ν e V allows for searches of heavier V , since the kinematic threshold will be M V = m π −m e rather than M V = m π − m µ , albeit at a suppressed rate. Figure 9 displays the expected number of V passing through the DUNE MPD, normalized to g 2 eµ = 1, assuming ten years of data collection with equal run time in neutrino and antineutrino modes. 9 We show contributions from charged mesons (π ± in purple and K ± in orange), as well as neutralmeson contributions (π 0 in red and η in blue) via loop-induced kinetic mixing, as well as the bremsstrahlung process pp → ppV in green. These rates assume that there is only loop-induced kinetic mixing, at the rate depicted in figure 8. The features apparent in the pp → ppV and η → γV production mechanisms (mostly near M V ≈ 200 MeV) are reflections of the variation of the kinetic mixing with V mass, as depicted figure 8.
For a L e − L τ vector boson, charged meson decays to electrons are chirally suppressed while production via charged meson decays with final-state muons rely on kinetic mixing. The loop-induced kinetic mixing is largest for L e − L τ over most of the V mass range, see figure 8, and is larger than (m e /m µ ) 2 for M V < ∼ 5 GeV. Thus, despite the small kinetic mixing, the processes depending on kinetic mixing dominate the production of V across all masses. Specifically, we note that the loop-suppressed process K ± → µ ± ν µ V is larger JHEP02(2020)174  Figure 10. The number of L e − L τ gauge bosons V traveling towards the DUNE MPD as a function of mass M V assuming ten years of data collection, and equal run time in neutrino and antineutrino modes. Solid lines display expected number due to charged meson decays -purple (π ± ) and orange (K ± ). Dot-dashed lines display expected number due to processes relying on kinetic mixing -red (π 0 decay), blue (η decay), green (pp → ppV bremsstrahlung), purple (π ± decays to muons), and orange (K ± decays to muons).
than the helicity-suppressed K ± → e ± ν e V . Figure 10 depicts the expected number of V passing through the DUNE MPD for L e − L τ gauge bosons as a function of mass M V from various channels. Finally, figure 11 depicts the number of V produced for a L µ − L τ gauge boson. Here, the small kinetic mixing (figure 8), causes the meson production mechanisms to final states containing muons to dominate as long as they are kinematically accessible. When the gauge boson mass is close to 770 MeV and it can mix with the SM ρ, the proton-proton bremsstrahlung process causes significant production, even with the small loop-induced kinetic mixing.

Decays of leptophilic gauge bosons
At tree level, L α − L β gauge bosons will decay to pairs of neutrinos (flavor α or β), or pairs of charged leptons V → + α − α and V → + β − β , assuming such a decay is kinematically accessible. Since we will not be able to observe neutrinos emerging from the decays of V , we sum over the flavors. The decay width, assuming massless neutrinos, is The decay into charged leptons, similar to those of dark photons in eq. (4.7), is When relevant, we also consider V decays via loop-induced kinetic mixing. Specifically, we will be interested in L µ − L τ gauge bosons that may decay into e + e − via the small kinetic JHEP02(2020)174  Figure 11. The number of L µ − L τ gauge bosons V traveling towards the DUNE MPD as a function of mass M V assuming ten years of data collection, and equal run time in neutrino and antineutrino modes. Solid lines display expected number due to charged meson decays -purple (π ± ) and orange (K ± ). Dot-dashed lines display expected number due to processes relying on kinetic mixing -red (π 0 decay), blue (η decay), and green (pp → ppV bremsstrahlung).

Sensitivity to leptophilic gauge bosons
Here we present the expected sensitivity of DUNE to leptophilic gauge bosons decaying to e + e − or µ + µ − final states within the DUNE MPD. In all cases, we assume ten years of data collection and equal run time in neutrino and antineutrino modes (in agreement with the production rates depicted in figures 9-11). For the e + e − final state, we choose to show regions of parameter space (the gauge boson mass and the new gauge coupling) for which greater than three or greater than ten signal events are expected in ten years. We discussed the backgrounds to a e + e − final state in section 4 and found that it should be background-free. As with the dark photon search, we assume here that the detector efficiency 10 is 100%. We reproduce existing limits from experimental searches from ref. [63]. Over the region of parameter space we are interested in, these are largely from electron beam dump experiments (E141, Orsay, and E137 providing the strongest constraints), as well as neutrinoelectron scattering via the exchange of the new gauge boson (with the strongest constraints coming from the TEXONO experiment). Ref. [89] recently explored the impact of a light L µ − L τ gauge boson on the cosmological history of the universe, and found that for certain combinations of masses and gauge couplings, large contributions to the number of effective relativistic species ∆N eff.   Figure 12 displays our expected DUNE MPD sensitivity to L e − L µ gauge bosons with masses between 1 MeV and 1 GeV. For points inside the blue (red) region, more than three (ten) signal events are expected. Given the discussion regarding backgrounds for A → e + e − in section 4.3, we expect the three event contour to represent the MPD sensitivity fairly well. For comparison, existing limits are depicted in grey. We see here that DUNE presents sensitivity that is nearly as powerful as electron beam dump experiments, but it is unlikely to improve on them. This is not terribly surprising since electron bremsstrahlung at existing experiments leads to a large flux of V . As stated above, the existing literature does not display limits below roughly 2 MeV -we expect that DUNE is no more powerful than the existing limits below this mass either, and include a dashed black line to extend the existing literature limits as we expect they would extend.
The existing limits on L e − L τ are nearly identical to those from L e − L µ , with the only exception being near M V ≈ 200 MeV, where decays into muons modify the lifetime of L e − L µ gauge bosons but not L e − L τ ones. Our expected DUNE MPD sensitivity, shown in figure 13, is slightly weaker, due to the reduced production discussed cf. figure 10. This results in the expected sensitivity being slightly weaker than, but comparable to, the existing limits.  Lastly, we depict the expected sensitivity to L µ −L τ gauge bosons in figure 14. Here, we are hindered in such a search for V → e + e − due to the kinetic-mixing-suppressed branching fraction to electrons. Regardless, such a search is complementary to the limit derived in ref. [89], and we emphasize that the DUNE MPD approach would be the first terrestrial search for L µ − L τ gauge bosons in this region of parameter space. As discussed around figure 11, for M V ≈ m ρ , the bremsstrahlung process is enhanced due to ρ − V mixing, leading to significant V production. At this mass, V may decay into µ + µ − unsuppressed. We find that for a very narrow range of parameter space, the DUNE MPD may expect to observe O(1) signal event with ten years of data collection. We indicate the region in which this occurs in green in figure 14 -we encourage a more detailed study of this region of parameter space for future work, specifically in determining the neutrino-related backgrounds to such a search. However, combining searches for V → µ + µ − both in the MPD and the LArTPC, as well as allowing the V to decay upstream of both detectors, could enhance the overall DUNE sensitivity in this region, allowing for searches for L µ − L τ gauge bosons in completely new areas of parameter space.

Dark Higgs boson
In the previous sections, we discussed sensitivity to semi-long-lived vector bosons that couple to the Standard Model via kinetic mixing with the SM hypercharge group or as a result of a new gauge symmetry under which some SM fermions are charged. Those scenarios are a subset of the commonly-studied renormalizable dark sector portals, where new particles may be added to the SM and the interaction between the dark sector and the SM is via a renormalizable operator.  Figure 14. Expected sensitivity to L µ − L τ gauge bosons at the DUNE MPD. The region in red (blue) indicates more than ten (three) signal events (V → e + e − ) expected. In grey, the region for which ∆N eff > 0.5 at the time of big bang nucleosynthesis [89]. In green is the region of parameter space for which V → µ + µ − may yield O(1) signal events, see text for further discussion.
A second portal option is the Higgs portal, where a new scalar ϕ interacts with the SM via one of two operators 12 L ⊃ δ |H| 2 ϕ, or L ⊃ λ |H| 2 |ϕ| 2 . (6.1) In the case of the first operator, mixing between H and ϕ occurs after electroweak symmetry breaking, while for the second operator mixing requires both H and ϕ to acquire a vacuum expectation value. This mixing may be expressed in terms of Lagrangian parameters. Here, however, we choose to use the phenomenological mixing parameter, which we denote as sin ϑ, and the ϕ mass, M ϕ , as the free parameters. When considering new vector bosons in sections 4 and 5, the couplings to the SM were given by the electric charge of the SM particle, reweighted by the parameter that characterizes the kinetic mixing between the SM and new U(1) gauge bosons. This led to roughly equivalent decays of A to electrons and muons, up to kinematic effects. For a dark Higgs, the coupling between ϕ and SM particles proceeds through the SM Higgs bosons, so couplings are related to SM particle masses, and are therefore hierarchical. This will drive production mechanisms to favor situations involving heavy quarks, which we will discuss in the following subsections. In particular, the dominant production modes are K → π ϕ and the decay modes are ϕ → e + e − , µ + µ − , π π with the ϕ decaying preferentially into the heaviest final state available. 12 The second option is more common in the literature -the first is super-renormalizable, but forbidden if ϕ carries any charge, SM or new-physics related. For this study, they are phenomenologically equivalent, as long as ϕ and the SM Higgs boson may mix. JHEP02(2020)174

Dark Higgs boson production
Refs. [95][96][97][98][99][100], among others, discuss the production of ϕ. Because ϕ couples via mixing with the SM Higgs scalar, its production will be largest in processes that allow for couplings to heavy quarks. At DUNE, the dominant contribution comes from the decays of K 0 L and K ± : the one-loop penguin diagram involving a W boson and top quark introduces ϕ − t couplings, which are large. While D-meson production is not negligible at DUNE, contributions of D meson decay into ϕ will be suppressed by elements of the CKM matrix and the bottom-quark Yukawa coupling. The meson with the largest branching fraction into ϕ is the B meson. However, B-meson production at DUNE is too small to be competitive with existing experiments, since √ s ≈ 12 GeV.
The matrix element for K decay into a light meson and a dark Higgs ϕ is [95][96][97] where γ 1 and γ 2 are both negligibly small. Here, and throughout, v ≈ 246 GeV is the vacuum expectation value of the neutral component of the Higgs field. The dominant contribution comes from the third term in eq. (6.2), given by contributions with a W boson and a loop involving up-type quarks. The largest contribution will come when i = t, the contribution from the top quark. Following refs. [95][96][97][98][99][100], we may write the branching fractions of interest as where ρ ϕ (x, y) is a dimensionless function, ρ ϕ (x, y) = 1 2 1 + x 2 + y 2 − 2(x + y + xy). (6.5) The branching fraction for K L is larger than for K ± due predominantly to its smaller total width. We disregard the contribution from K S → π 0 ϕ, which has a small branching fraction due to the large total width of K S . One may also consider ϕ produced via decays of other mesons. In comparison with eqs. (6.3)-(6.4), the branching fraction of a B meson into a strange meson and a ϕ, for small ϑ, is [100] Br(B → X s ϕ) 5.
Dark Higgs Boson π ± → e ± νϕ, are as small as 10 −9 sin 2 ϑ. The dominant sensitivity at the DUNE MPD then, will come from kaon-related production of dark Higgs bosons.
In figure 15, we present the fraction of ϕ that are produced in a given decay and are traveling in the direction of the MPD when they are produced. Several features are worthy of note. Reversing the horn current and switching from neutrino to antineutrino mode nearly perfectly interchanges the acceptance fraction of ϕ coming from K + and K − decays -dashed lines are coincident with solid lines of the opposite color. In the limit M ϕ → 0, the acceptance fraction of the focused distribution (K + in neutrino mode, K − in antineutrino mode) is almost, but not quite, double that of the defocused distribution. If analyzing decays of charged pions instead of charged kaons, this ratio would be significantly larger; DUNE is designed to focus π ± and defocus π ∓ , while kaon focusing is not the main priority of the experiment. Naïvely, we expect that the K 0 L -produced ϕ would have an acceptance fraction somewhere between the focused and defocused distributions, as they are neither focused nor unfocused. This is not the case, as evident by the green line in figure 15, and is because a significant fraction (roughly 25%) of K 0 L live long enough to reach the rock at the end of the DUNE beam decay pipe. Those K 0 L decay at rest, and a very small fraction of the isotropic decay products are pointing in the direction of the MPD.
With this information about the direction of travel of the dark Higgs, we simulate the decays with the branching fractions from eqs. (6.3) and (6.4) and determine the number of ϕ particles from these decays that pass through the DUNE MPD, assuming five years of JHEP02(2020)174 data collection in both neutrino and antineutrino mode. This is depicted in the right panel of figure 15. Here, solid (dashed) lines indicate (anti)neutrino mode, and blue, red, and green lines indicate K + , K − , and K 0 L decays, respectively. Because K 0 L are unaffected by the focusing horns, there is no distinction between this production mechanism for neutrino and antineutrino modes. The K 0 L distribution also extends to slightly larger masses, as m K ± − m π ± ≈ 354 MeV whereas m K 0 L − m π 0 ≈ 394 MeV.

Dark Higgs boson decays
For very light ϕ, M ϕ < 2m π , the decay widths of ϕ are given by For masses above 2m π , but below a few GeV, the width to mesons is hard to determine due to strong QCD effects and the effects of meson resonances [101]. For the mass range of interest here, M ϕ 0.5 GeV, the only meson decay open is to a pair of pions. This decay width can be expressed in terms of a hadronic matrix elements as, Here |G(M 2 ϕ )| is the (dimensionless, in contrast with some results in the literature) transition amplitude for the process ϕ → ππ. To estimate G(M 2 ϕ ), we follow the discussion of ref. [102], which we now briefly outline. This transition amplitude may be expressed, in the limit of isospin conservation, in terms of three independent form factors, which are functions of s, the invariant mass-squared of the outgoing pion pair: θ π (s) + 7 9 [Γ π (s) + ∆ π (s)] . (6.9) At next-to-leading order in chiral perturbation theory, 13 θ π (s) = 1 + 2m 2 π s (1 + ψ(s)) + b θ s, (6.10) In these formulae, the dimensionful quantities b θ , b Γ , and b ∆ are extracted from derivatives of the form factors at s = 0. The dimensionless parameter d F is related to the strangeness content of the pion. Finally, the function ψ(s) is defined as (6.13) 13 Note that these quantities are, like G, dimensionless, and differ in definition from those in ref. [102]. where F π = 130 MeV is the pion decay constant and κ ≡ (1 − 4m 2 π /s) 1/2 is a kinematic factor. Ref. [102] extracts the relevant parameters from measurements of pion-pion scattering. In our calculations, we use d F = 0.09, b θ = 2.7 GeV −2 , b Γ = 2.6 GeV −2 , and b ∆ = 3.3 GeV −2 . As discussed in ref. [102], this leads to an enhancement of the partial width into pions by a factor of approximately 4.4 at M ϕ = 500 MeV, relative to the leadingorder result from ref. [103]. This leading order result is obtained by setting b θ,Γ,∆ = 0, d F = 0, and ϕ(s) = 0 in eqs. (6.10)-(6.12).
Instead of depicting branching fractions of different final states as a function of mass, we choose to present instead the total lifetime of ϕ as a function of its mass in figure 16, keeping in mind that, given its hierarchical couplings to SM fermions, ϕ decay is dominated by the heaviest allowed final state. For light masses, the lifetime (assuming sin 2 ϑ = 10 −8 ) is significantly longer than the distance to the DUNE near detector. For masses between roughly 200 − 400 MeV, the lifetime is of order of the distance between the DUNE target and the MPD, meaning a ϕ decay in the detector becomes significantly more likely.

Backgrounds and sensitivity
In section 4, we discussed the possible backgrounds for searches of final states involving e + e − , µ + µ − , or pion pairs and argued that, especially when utilizing the invariant mass peak of µ + µ − , we may safely consider the expected backgrounds in the DUNE MPD to be small. The signal for a dark Higgs scalar decaying into charged leptons will appear practically identical to that for a dark photon decaying into charged leptons. We discuss here two possible ways of disentangling the two hypotheses.
In contrast to the dark photon scenario studied in section 4, where, for a particular mass M A , one expects a certain ratio of different final states, such as 60% of events being A → e + e − and 40% being A → µ + µ − , here we expect nearly every event to be in a single search channel because the individual partial widths for ϕ decay are so hierarchical. For lighter dark Higgs/vector bosons where decays into µ + µ − are kinematically forbidden, this strategy will not allow us to disentangle the two hypotheses. Secondly, because the ϕ flux JHEP02(2020)174 is coming from decays of both charged and neutral kaons, where A comes from decays of neutral mesons (or continuum processes), the A spectrum should be less focused than the ϕ one. One proposal for the DUNE Near Detector Complex is that the liquid and gas TPCs move off axis (the DUNE-PRISM proposal [53]) to measure different portions of the neutrino spectrum. In principle, one can calculate the distribution of ϕ and A as a function of how off-axis the near detector is, and search for such a shape when operating on-and off-axis, making it possible to distinguish the different new physics hypotheses.
Combining the fluxes and the decay width of ϕ, figure 17 depicts in blue (red) regions of M ϕ vs. sin 2 ϑ parameter space for which we expect at least ten (one hundred) total decays (into e + e − , µ + µ − , or ππ) in ten years of data collection at DUNE. For comparison, existing limits are depicted in grey [97,99,104]. We also include a limit from searches for charged kaons decaying to a charged pion and an invisible particle, reported in ref. [105], translated into M ϕ vs. sin 2 ϑ parameter space.
We see here that the DUNE MPD will allow for an improvement on current limits between roughly 20 − 350 MeV, improving on the BNL search and CHARM for this wide range of masses. A DUNE MPD search may also allow us to comprehensively search in the window still open between searches of this type (displaced decays) and limits from considerations of supernova luminosity from Supernova 1987A [99]. A recent study [21] of the sensitivity of Fermilab's short baseline neutrino detectors (e.g., ICARUS and SBND) to dark Higgs models projects bounds on sin 2 ϑ a factor of a few weaker than those anticipated at DUNE.

JHEP02(2020)174 7 Heavy Neutral Leptons
In this section, we consider the existence of new fermions that are not charged under the SM gauge interactions. Assuming these are not charged under any other gauge symmetries, they couple, at the renormalizable level, to lepton doublets and the Higgs doublet. After electroweak symmetry breaking, for one new fermion N , the neutrino weak-interaction eigenstates ν e , ν µ , ν τ can be written as linear combinations of the neutrino mass eigenstates ν 1 , ν 2 , ν 3 , N : Here, we refer to N as a heavy neutrino or a heavy neutral lepton (HNL). 14 We are interested in 1 MeV M N 1 GeV, when the HNLs can be produced by the same mechanisms -meson and lepton decays -as the light neutrinos in the DUNE target and decay pipe and are still heavy enough to decay into charged SM particles. The production and decay rates of such HNLs are governed by the weak interactions and the mixing parameters U αN . Furthermore, the same parameters allow for a variety of decay modes for the HNLs, depending on the HNL mass. These production and decay mechanisms allow us to search for HNLs in the DUNE MPD. For simplicity, we will consider that N couples to only one charged lepton flavor at a time, i.e., we will assume that at most one of U eN , U µN , U τ N is nonzero for a given analysis. Section 7.1 discusses the production channels for N considered here and section 7.2 discusses N decays that could be detected in the DUNE MPD. In section 7.3, we discuss backgrounds associated with each signal and present the expected sensitivity for HNLs coupling to the SM via nonzero |U eN | 2 , |U µN | 2 , or |U τ N | 2 . Finally, section 7.4 discusses how, if such an HNL were detected in DUNE, measurements of the decay products could determine whether N is a Dirac or Majorana particle, hence whether lepton-number symmetry is violated in nature.
Ref. [106] recently explored the DUNE near detector's sensitivity to HNLs in great detail. Their analysis utilized the combined decay volume of the LArTPC and the MPD, and did not separate them for the sake of background reductions, as we do. Nevertheless, our results and theirs, when asking similar questions, are qualitatively equivalent (specifically, for example, figures 24, 25, and 26 in our work and figure 5 of ref. [106]). 15 Whereas ref. [106] focused on connecting such HNLs to the origin of the light neutrino masses, we choose a different lens and focus on diagnosing whether such HNLs are Dirac or Majorana particles.
14 Such gauge-singlet fermions which mix with light neutrinos are often referred to as "sterile neutrinos." We avoid this terminology, as it often refers to mixing that coherently modifies the active neutrino oscillations. 15 We also emphasize that the simulation of HNL production in ref. [106] and our work is different for the bulk of results -ref. [106] simulated HNL fluxes by reweighting the well-studied light neutrino fluxes at the DUNE near detector hall, where we simulate HNL fluxes from the parent mesons. The fact that the results agree validates the two different approaches.

Production of Heavy Neutral Leptons
Refs. [52,66,[106][107][108][109][110][111] have explored the possibility of HNL production in fixed-target experiments such as DUNE in some detail. In such a scenario, the relevant production mechanisms are leptonic decays of mesons, either via two-body decays (such as π + → µ + N ) or three-body decays (such as K + → π 0 e + N ). Section 3 discusses how we simulate the generation of relevant mesons (in this case, π ± , K ± , D ± , and D ± s ) and the two-and threebody decays. We also include secondary production from the decays of τ ± produced in two-body D ± and D ± s decays, discussed in some detail below. For large enough beam energy, the direct production of N via qq → W * → N (where q and q are partons) may be considered. Additionally, B mesons can also be produced in the target, contributing further to HNL production. Given that DUNE uses 120 GeV protons on target, however, HNL production from both of these channels is negligibly small.
The possible decay channels resulting in a flux of N will depend on which coupling (U eN , U µN , or U τ N ) is considered -we discuss each case in turn below. In general, we express the two-body meson decay branching fractions into N via where Br(m + → + α ν) is the two-body branching fraction into SM neutrinos, and ρ N (x, y) is a function that accounts for the available phase space in the decay, including helicity suppression. 16 This function is These decays are isotropic in the rest frame of the parent meson. We also calculate the polarization asymmetry of N stemming from this decay; we will discuss this and how it affects the determination of whether N is a Dirac or Majorana fermion in section 7.4. The expressions for three-body decays are more complicated than eq. (7.2) -we direct the reader to refs. [111,112] for details. See section 3 for the approximations made in our three-body decay simulations and their validity.
Production with electron coupling: for HNLs that couple only to electrons via |U eN | 2 , the N production processes of interest are as follows. For simplicity, we focus our discussion on positively charged mesons decaying, which would be the dominant production channels in neutrino mode at DUNE. The charge-conjugated processes dominate in antineutrino mode. The meson decays considered are π + → e + N , K + → e + N , K + → π 0 e + N , D + → e + N , D + → π 0 e + N , D + → K 0 e + N , D + s → e + N , the most dominant channels for every possible M N considered. Figure 18     Production with muon coupling: as in the electron-coupled channel, we consider decays when only |U µN | 2 is nonzero. The list is as follows: π + → µ + N , K + → µ + N , The N flux from each of these contributions is depicted in figure 19. Again, we separate HNL production in neutrino mode (left panel) from antineutrino mode (right panel).  Dark colors indicate positively-charged particle decays, where faint ones indicate negatively-charged particle decays.

Neutrino or Antineutrino Mode
Production with tau coupling: lastly, we list the meson decays in which HNLs are produced when the only nonzero coupling is |U τ N | 2 . Because these decays involve the τ lepton, we need only consider decays of mesons with mass greater than m τ : D + → τ + N , D + s → τ + N . Because D s is heavier than D and its branching fraction into τ + ν τ is larger, D s decays contribute more than D decays; see the orange and green curves in figure 20. We also consider decays of secondary τ ± produced in two-body D/D s decays. Here, we assume that both the D/D s and subsequent τ decays are sufficiently prompt that the magnetic horns focus neither. We consider the following decays of τ leptons: τ → N π, τ → N K, τ → N ρ, τ → N ν e e, τ → N ν µ µ. The number of τ -produced HNLs corresponds to the purple curves in figure 20. These decays provide sensitivity to significantly heavier τ -coupled HNL than D/D s decays on their own -nearly up to the τ mass. Because the focusing has no effect on the D mesons or τ leptons, figure 20 represents the N flux for five years of data collection in either neutrino or antineutrino mode.

Decays of Heavy Neutral Leptons
Depending on the flavor to which an HNL couples and M N , it will be able to decay into a variety of final states. At the DUNE MPD, several of these final states will either be fully invisible or swamped by backgrounds. We classify these decay channels as irrelevant and will combine their branching fractions/decay widths for simplicity. Our relevant final states are those with at least one charged lepton.
Refs. [52,112,113] give the partial widths of HNL decays of interest. For a given |U αN | 2 and M N , we determine the available decay channels, the total rest-frame lifetime JHEP02(2020)174 Figure 21. The branching ratio of an electron-coupled HNL N into different final states as a function of its mass M N . The grey curve is a collection of decays considered to be irrelevant, as described in the text.
of N and its branching fractions into channels of interest. We list the decay channels of interest for each coupling below.
If HNLs are Dirac particles so that lepton number is conserved, then the N produced from positively charged meson decay (along with + α ) have lepton number +1 and their decays must conserve lepton number. For instance, the decay N → µ − π + would exist, but N → µ + π − would be forbidden. If N is a Majorana fermion so that lepton number is violated, then both decays would occur, each with the same rate as that for N → µ − π + when N is a Dirac particle. When presenting sensitivities, 17 we assume N is a Dirac fermion, which corresponds to a more conservative reach to lower |U αN | 2 than if it were a Majorana fermion. We will discuss the nuances of the differences between Dirac and Majorana HNLs in section 7.4.

Relevant and irrelevant decay channels
In creating the following list of possible decays, we have continued to assume for simplicity that only one of the mixing matrix elements |U eN |, |U µN |, and |U τ N | is nonzero at a time. In the list, the daughter ν is one of the three light neutrino mass eigenstates, ν 1 , ν 2 , or ν 3 . Of course, in practice, this daughter neutrino will not be observed. This list also assumes that N is a Dirac particle, and only gives the decays for the L = +1 particle N . We include the decays of N in our simulations.
If |U eN | is nonzero: Figure 22. The branching ratio of a muon-coupled HNL N into different final states as a function of its mass M N . The grey curve is a collection of decays considered to be irrelevant, as described in the text.
• Irrelevant: If |U µN | is nonzero: • Relevant: If |U τ N | is nonzero: • Relevant: We expect all of the decays we deemed irrelevant to be background-dominated; for instance, the N → νπ 0 channel will have a large NCπ 0 background from incoming beam JHEP02(2020)174 Figure 23. The branching ratio of a tau-coupled HNL N into different final states as a function of its mass M N . The grey curve is a collection of decays considered to be irrelevant, as described in the text. Because we do not consider any M N > m τ , we do not present any branching fractions in this region.
neutrinos. Additionally, since all irrelevant decays have at least one outgoing neutrino, in general the invariant mass of N cannot be accurately reconstructed, meaning that this analysis handle is not available to distinguish signal from background.

Sensitivity to Heavy Neutral Leptons
Here we discuss the sensitivity as a function of the HNL mass M N and mixing |U αN | 2 for the different possible N decay channels. We also discuss existing limits in this parameter space and how DUNE will be able to improve on them.

Backgrounds for relevant HNL decay channels
We will discuss each of our relevant decay channels in turn, and their associated backgrounds.
Searches for N → νe + e − : this final state is relevant regardless of which |U αN | 2 is nonzero and will contribute nontrivially to the resulting sensitivity of each coupling. Because of the final-state neutrino, the reconstructed e + e − pair will neither point perfectly in the beam direction nor reconstruct the invariant mass of N . However, as discussed in section 4, the expected number of neutrino-related e + e − background events is very small for ten years of data collection. 18 In the DUNE MPD, roughly 20 e + e − -trident events are JHEP02(2020)174 also expected [70,71]. This rate, and the kinematics of such events, may be predicted in a data-driven way by using the sample of ∼ 1000 events expected in the liquid argon TPC. Therefore, we expect that 10 signal events of N → νe + e − are a statistically significant sample for ten years of data collection.
Searches for N → νe − µ + or νe + µ − : as with N → νe + e − , these decays are expected to have very low backgrounds, mostly associated with the misidentification of a charged pion as an electron or muon. Particle identification would help reduce this background, as would the background reduction strategies discussed in section 4. We expect that 10 signal events will correspond to a statistically significant sample in these channels.
Searches for N → νµ + µ − : like the dark photon search with µ + µ − , this signal channel will be competing against ν µ charged-current backgrounds with single pion emission where the pion is misidentified to be a muon. In the case of the dark photon, the µ + µ − pair reconstructed the dark photon mass; here, they will not reconstruct the N mass. Therefore, we would not expect to find a peak on top of a smooth background in the µ + µ − invariant mass, as we would for a dark photon. We expect O(100) background events in ten years of data collection with the DUNE MPD. For this reason, we expect that at least 20 signal events, instead of 10 events, will correspond to a statistically significant sample in this channel.
Searches for two-body final states: the remainder of our relevant channels consist of decays of N with two particles in the final state, one being the charged lepton and the other being a charged meson. While certain channels will have relatively large neutrino-related scattering backgrounds (for instance, ν e CC1π scattering as a background for N → e − π + ), the invariant mass reconstruction of N and the direction of the total final-state particle momentum will allow for significant background suppression. Channels with a final-state meson that decays within the detector volume may allow for more distinguishing signatures. For instance, the channel N → µ − ρ + will result in the ρ + decaying to π + π 0 . This final state would be relatively easy to reconstruct. In the DUNE MPD, π 0 s decay to photons which deposit their energy in the ECAL, allowing the π + π 0 system to reconstruct the ρ + mass; the µ − π 0 π + system would have an invariant mass peak corresponding to the N mass of interest. As with most other relevant channels, we expect 10 signal events to be statistically significant here.

Existing limits
For each nonzero |U αN | 2 , we compare the sensitivity reach of the DUNE MPD against existing limits. We assume that U αN is the only coupling between the HNL and the SM; HNL production and decay are then completely determined by M N and |U αN | 2 . Additional interactions between the HNL and SM particles or interactions between N and other new particles modify the N lifetime, and therefore its probability of decaying in the detector, in a nontrivial way.
should still be within a degree or so of the direction of the beam, allowing us to reduce neutrino-related backgrounds significantly.
For muon-coupled HNL, many of the same probes apply. In addition to the probes discussed for e-coupled HNL, the NuTeV experiment [135] contributes limits for 400 MeV M N 2 GeV. The Brookhaven E949 experiment provides constraints in the ∼175 MeV -∼300 MeV region with precision measurements of the decay K → µν [105]. Measurements of the µ Michel decay spectrum [124] and searches for π → µN at PSI [136] provide constraints at low M N , K → µν decays provide constraints for masses between roughly 70 MeV and 175 MeV [137]. Recently, ref. [138] proposed a search using atmospheric-produced HNL that can travel to and decay inside Super-Kamiokande and ref. [139] proposed a similar search using IceCube. Moreover, the MicroBooNE collaboration has recently reported a constraint in ref. [140]. These last three limits are slightly weaker than those discussed above.
HNLs that couple purely to the tau are relatively less constrained. Existing searches from CHARM [141] and DELPHI [133] provide the most stringent constraints for the mass ranges of interest. IceCube could have sensitivity to τ -coupled HNL via double-bang signatures [142]; however, no such analysis has yet been performed by the collaboration. Current B factories could also reanalyze data to search for such HNLs using kinematical arguments [143]. The existence of such HNLs can also be investigated in the NA62 experiment [109], the forthcoming FASER experiment [74,144], or the proposed SHiP experiment [145,146].

Expected DUNE MPD sensitivity
We present the expected sensitivity of the DUNE MPD assuming only one nonzero |U αN | 2 at a time. The search channels are as listed in section 7.2.1 and the number of signal events deemed statistically significant is as discussed in section 7.3.1. For simplicity, we show only the search channels which drive the overall sensitivity at DUNE -the collection of channels that dominate the others for some range of M N . For a given M N of interest, a combined search for all kinematically-accessible channels would lead to slightly greater sensitivity than what we present.
As discussed above, in determining the number of signal events for a given point in parameter space, and therefore the DUNE MPD sensitivity, we assume that the HNL is a Dirac fermion. If it were a Majorana fermion, then its total decay width would be a factor JHEP02(2020)174  Figure 24. Expected sensitivity of the DUNE MPD to an electron-coupled HNL, assuming ten years of data collection (equal time in neutrino and antineutrino modes), assuming that the HNL is a Dirac fermion. Specific decay channels are considered: N/N → νe + e − (red), N/N → νe ± µ ∓ (blue), N/N → e ± π ∓ (purple), and N/N → e ± ρ ∓ (orange). Each channel corresponds to ten signal events, as discussed in section 7.3.1. See text for more detail.
of two larger, leading to twice as many signal events for a given point in parameter space. 19 The difference in sensitivity under the Dirac and Majorana assumptions was explored in ref. [106]. In section 7.4, we will explore how a discovered HNL could be diagnosed to be Dirac or Majorana. Figure 24 depicts the expected DUNE MPD sensitivity to an electron-coupled HNL assuming ten years of data collection with equal time in neutrino and antineutrino modes. The channels included here are N, N → νe + e − (red); N → νe − µ + or N → νe + µ − (blue); N → e − π + or N → e + π − (purple); and N → e − ρ + or N → e + ρ − (orange). Three of these search channels are sensitive to regions of the parameter space that lie beyond what is constrained by existing limits (shown in grey) for some range of M N . Most notably, in the absence of a discovery, the N → e ± π ∓ channel (which should produce a relatively clean signal in the DUNE MPD) will improve limits on |U eN | 2 , compared to those from CHARM and BEBC, by roughly an order of magnitude. Figure 25, depicts the expected DUNE MPD sensitivity if the HNL couples to the muon. The channels of interest are similar to those in figure 24 -N → νe + e − (red), N → νµ ± e ∓ (blue), N → µ ± π ∓ (purple), and N → µ ± ρ ∓ (orange). As with the electroncoupled HNL, in the absence of a discovery, searches in the DUNE MPD for a muon-coupled HNL will improve on existing limits for a wide range of M N . We highlight the sensitivity JHEP02(2020)174  Figure 25. Expected sensitivity of the DUNE Multi-Purpose Near Detector to a muon-coupled heavy neutral lepton, assuming ten years of data collection (equal time in neutrino and antineutrino modes). We assume that the HNL is a Dirac fermion. Specific decay channels are considered: N/N → νe + e − (red), N/N → νµ ± e ∓ (blue), N/N → µ ± π ∓ (purple), and N/N → µ ± ρ ∓ (orange). Each channel corresponds to ten signal events, as discussed in section 7.3.1. See text for more detail.
to low-mass M N driven by the search for N → νe + e − . We discuss this search channel and how an HNL decaying in this way could be diagnosed to be a Dirac or Majorana fermion in section 7.4.
Finally, figure 26 displays the expected sensitivity to τ -coupled HNL at the DUNE MPD. Because no production mechanism exists for HNLs with M N > m τ , we have no access to any final states with τ ± and a charged meson. The DUNE MPD sensitivity to a tau-coupled HNL is independent of how the data collection time is divided between neutrino and antineutrino modes since the production mechanisms are not impacted by the magnetic focusing horns. The two channels we display are N → νe + e − (red) and N → νµ + µ − (green). Due to background considerations discussed above, the N → νµ + µ − curve corresponds to 20 signal events. In the absence of a discovery, the N → νe + e − channel will provide a more powerful constraint than existing CHARM constraints [141], and a combination of both channels will allow us to search for tau-coupled HNL in a window where the CHARM and DELPHI [133] searches do not have much sensitivity. This same window can be probed, as discussed above, by IceCube [142], B factories [143], NA62 [109], and FASER [144]. Such searches would be complementary in probing this gap with the DUNE MPD.

Deducing the nature of the HNL -Dirac vs. Majorana
Above, we assumed that HNLs are Dirac fermions and that lepton number is conserved. In that case, when N is assumed to be the HNL produced in association with a positively JHEP02(2020)174  Figure 26. Expected sensitivity of the DUNE MPD to a τ -coupled heavy neutral lepton, assuming ten years of data collection (how this time is apportioned between neutrino and antineutrino modes is irrelevant). We assume that the HNL is a Dirac fermion. We show sensitivity to the channels N/N → νe + e − (red) and N/N → νµ + µ − (green), compared against existing limits from CHARM [141] and DELPHI [133]. For the final state νe + e − , the red region encompasses parameter space for which ten signal events are expected, whereas for the final state νµ + µ − , the green region indicates twenty signal events, due to background considerations. See text for more detail.
charged lepton and having its own lepton number L = 1, its decays always contain negatively charged leptons or neutrinos. On the other hand, if HNLs are Majorana fermions, then they will also violate lepton number in their decays, i.e., the N produced via positively charged meson decay may then decay, for instance, into either N → µ − π + or N → µ + π − , with equal probability. Since the magnetized DUNE MPD has charge-and particle-identification capabilities, it will be able to search for differences between the rates of decay into these final states. We will show in section 7.4.2 that perfect charge/particle identification are not required to perform this analysis.
As explored in refs. [146][147][148], Majorana HNLs also differ from Dirac ones in their decay kinematics. When the HNLs are produced in weak-interaction meson decays, they will in general have a net polarization, potentially leading to asymmetry in their rest-frame decay angular distributions. Different rest-frame angular distributions will be reflected in different decay-product energy spectra in the lab frame. Even with a charge-blind detector, these different spectra could be distinguished.
We explore the potential of determining whether HNLs that are discovered by DUNE are Dirac or Majorana fermions. Two classes of HNL final states are considered: those with unobservable lepton number due to the presence of a neutrino in the final state (e.g., N → νe + e − or N → νπ 0 ) and those with only charged particles, where the total lepton JHEP02(2020)174 number may be measured (e.g., N → µ ± π ∓ ).
Because of leptonic mixing, if we can establish experimentally that an HNL N is a Majorana particle, then it will follow that all neutrinos, including the light ones, are Majorana particles. For, if N is a Majorana particle, its exchange will lead, via mixing, to diagrams that give Majorana masses to all the neutrino mass eigenstates, and once these mass eigenstates have Majorana masses, it follows that they will be Majorana particles. Thus, establishing that an HNL N is a Majorana particle would have far-reaching implications.

Final states with unobservable lepton number
If there are SM neutrinos in the final state of an HNL decay, then it is not possible to determine the state's total lepton number, and counting these events will not readily indicate whether the HNL is a Dirac or Majorana fermion. Instead of identifying the final state's lepton number, one must rely exclusively on kinematic information.
Searches for N → νπ 0 : consider, for instance, the decays N → νπ 0 and N → νπ 0 for a Dirac HNL, and the decay N → νπ 0 for a Majorana HNL. In principle, the distributions of π 0 lab-frame energies and angles relative to the beam direction are different between these two hypotheses; a large enough sample of decays could allow these distributions to be distinguished. For illustrative purposes, we explore this decay channel in detail here, although, as discussed in section 7.2.1, we expect this channel to have a large background from neutral-current events.
If the neutrinos, including N , are Majorana particles, the decays N → νπ 0 will produce an isotropic distribution of π 0 s in the rest frame of N [148]. In contrast, for Dirac N , the decay N → νπ 0 is maximally anisotropic [148]. However, the DUNE beam is not a perfect environment for such searches -since there will be some contamination of N particles in the beam. Even if N is a Dirac fermion, the N decay to νπ 0 , with the opposite anisotropy, will hinder the separation between the Majorana and Dirac fermion hypotheses. For concreteness, we explore the kinematics of such decays assuming M N = 350 MeV. We determine the flux (and polarization) of HNLs from K + and K − decays when the DUNE facility is operated in neutrino mode. 20 When simulating the HNL decays, we assume they are either Majorana (i.e., the decay products are emitted isotropically in the rest frame) or Dirac (i.e., the decay products are emitted anisotropically in the rest frame) fermions, and so determine the lab-frame distribution of π 0 energies. These resulting lab-frame π 0 energy distributions are depicted in figure 27 for M N = 350 MeV; the Majorana (Dirac) hypothesis results are depicted in red (blue). The inset depicts the lab-frame distribution of N energies as they travel to the detector. We note that the HNLs are highly boosted; most of them have energies between 5 and 10 GeV, corresponding to boost factors of ∼10−20. It is apparent in figure 27 that an extremely large sample of events (even neglecting background events) would be required to use this sample to distinguish between Dirac and Majorana HNLs.
The authors of ref. [148] explored this decay channel as a possible means of distinguishing Dirac from Majorana HNLs in a toy scenario. Our analysis differs from theirs in two JHEP02(2020)174 main ways. First, we considered the decays of both N and N , instead of either one or the other. We did, however, determine that the impact of this wrong-lepton-number contamination is not very significant. Second, and more importantly, the toy analysis performed in ref. [148] considered a flux of N with boost factors between, roughly, 2−4, significantly lower than what would occur for a 350 MeV HNL at DUNE. More strongly boosted HNLs yield predicted energy distributions that are substantially more difficult to distinguish.
We conclude from this that, if DUNE discovers an HNL that is produced predominantly from charged kaon decay, then such beam-neutrino environments are not ideal for studying the intrinsic nature of the HNL mass in the final state νπ 0 . Scenarios in which N are generated from only one charge of kaon, and in which the N are significantly less boosted, would be significantly more propitious. We leave the study of such scenarios to future work.
Searches for N → νe + e − : in section 7.3, we showed that the DUNE MPD is sensitive to a large region of currently unexplored parameter space by searching for the decay channel N → νe + e − . This is particularly true for the scenarios in which N couples only to the muon (figure 25) or only to the tau (figure 26). Compared to the final state νπ 0 , this has significantly lower background. While it suffers from having unidentifiable lepton number, its kinematics are different for Dirac and Majorana N . Ref. [106] explored this decay in more detail. We anticipate that, because the relative angular distributions for Dirac/Majorana N are not as discrepant (in the N rest frame) as in the decay N → νπ 0 , such Majorana vs. Dirac fermion differences would be difficult to measure in the DUNE environment. We leave a more detailed study of these three-body decays to future work.

Final states with only charged particles
In the previous subsection, we focused on decay modes in which the final-state lepton number is unobservable due to the presence of SM neutrinos. We now explore final states with only charged particles where, if particles and their charges are identified, the lepton number will be determined. In a perfect world, assuming HNLs are Dirac fermions, an experiment operating in neutrino mode would produce HNLs with L = +1 directed toward a detector and any HNLs with L = −1 would be diverted away. Then, using measurements of charge/kinematics, one could deduce that the decaying HNLs have L = +1. Obviously, an experiment like DUNE is not perfect; even with the light neutrinos, a certain level of antineutrino contamination exists in neutrino mode, and vice versa.
Assuming the HNLs are Dirac fermions, and using the production rates depicted in figure 18-20, we calculate the ratio of two fluxes as a function of M N (for each different coupling): the number of HNLs with L = +1 reaching the MPD to the number of HNLs with L = −1 also reaching the MPD. In a perfect world, this ratio is ∞ in neutrino mode and 0 in antineutrino mode. Figure 28 presents the ratios in neutrino mode (blue) and antineutrino mode (red) for e-(left), µ-(center), and τ -(right) coupled HNLs. When the distinction between neutrino and antineutrino modes is irrelevant (because the parent particles are unfocused), we use purple lines. If HNLs were Majorana fermions, then they would decay as if they have L = +1 and L = −1 in equal abundance -the measured ratio would be 1.
Several features in figure 28 are notable. First, the expected ratio typically deviates from 1 more in neutrino mode than in antineutrino mode, making such a measurement more favorable. This is not because of the ability of the beam to focus one polarity of mesons over the other; it is simply because positively charged light mesons are produced more frequently in the beam than negatively charged ones (see table 4  abundance. This explains why, for all couplings, for heavy enough M N , both neutrino and antineutrino mode predict more L = −1 HNLs than L = +1 ones. Additionally, the ratio in neutrino mode is typically O(few). Even with perfect event-by-event identification of whether the HNL has L = +1 or L = −1, to determine whether the detected HNL is a Dirac (where the ratio is predicted by figure 28) or Majorana fermion (where the ratio should be 1) will require many events.
To demonstrate the power of measuring this ratio, consider a muon-coupled HNL with a mass M N = 250 MeV. We calculate the expected number of N → µ ± π ∓ events (at this mass, the branching ratio N → µ ± π ∓ is approximately 20% and this is the largest of the relevant channels) for five years of neutrino mode operation as a function of the mixing |U µN | 2 . For a given |U µN | 2 , we perform a number of pseudoexperiments where we assume that the HNL is a Majorana fermion which therefore decays to µ − π + and µ + π − with equal probability. These event rates are defined as N (L = +1) and N (L = −1), respectively, and are drawn from a Poisson distribution. We then determine the 68.3%, 95%, and 99% credible regions for the ratio N (L = +1)/N (L = −1) after performing these pseudoexperiments. Figure 29 displays the result of this procedure assuming five years of data collection in neutrino mode. As |U µN | 2 increases, so does the expected number of events, and it becomes easier to measure that the ratio is 1, as predicted by the Majorana HNL hypothesis. If HNLs are Dirac fermions, and if we are operating in neutrino mode, then the expected ratio is approximately 3 (see the center panel of figure 28). The Dirac hypothesis can be robustly excluded (if every event can be properly identified as L = +1 or L = −1) for JHEP02(2020)174 |U µN | 2 2 × 10 −9 . On the top axis of figure 29, we give the number of expected signal events for a Majorana HNL given the corresponding value of |U µN | 2 . Figure 29 demonstrates the ideal ability of the DUNE MPD to determine whether HNLs are Dirac or Majorana fermions if, on an event-by-event basis, the signal events can be perfectly identified as having lepton number +1 or −1 (by identifying the sign of the charged lepton in the final state). In the real world, a detector will misidentify these final states a finite fraction of the time. This is tantamount to redistributing the events between the samples, limiting the ability of the detector to separate the Dirac and Majorana hypotheses.
To study the extent to which the MPD could discriminate between the Dirac and Majorana hypotheses using these types of final states (without using kinematic information), we separately study the decays N → µ ∓ π ± , N → e ∓ π ± , and N → µ ∓ ρ ± . For each final state, we calculate the expected number of decays to each charge state assuming HNLs are Majorana particles as a function of the HNL mass and mixing. We perform a frequentist analysis to address how distinct the Majorana-assumed data are from data generated under the Dirac hypothesis. 21 Since such an analysis would take place only after an initial discovery of N , we optimistically assume that, with enough events, M N will be nearly perfectly measured. If N decays completely visibly, then its mass may be measured from the invariant mass of the final state particles. While the detector is not perfect, we assume that it can attain O(10 MeV) precision, corresponding to nearly perfect knowledge of M N as far as this analysis is concerned.
For the Majorana hypothesis, the expected numbers of L = +1 and L = −1 events are equal in both neutrino and antineutrino modes; we denote these expectations as n ν,M ± and nν ,M ± . For the Dirac hypothesis, there is a (mass-dependent) ratio between the expected numbers of L = +1 and L = −1 events in both running modes (see figure 28). We denote the numbers of events predicted by the Dirac hypothesis, as a function of the mixing U D = U αN , as n ν,D ± (U D ) and nν ,D ± (U D ). We use a likelihood ratio to compare the two hypotheses; in practice, we consider the difference in their log-likelihoods. For observed numbers of events N ν ± and predictions n ν ± , the log-likelihood is log L ν (n ν + , n ν − ) = −n ν + + N ν + log n ν and similarly for the antineutrino mode. For the Majorana hypothesis, the likelihood is clearly maximized for n ν,M ± = N ν ± , since this is our null hypothesis. For the Dirac hypothesis, we maximize the likelihood with respect to U D and compare the likelihood at this best-fit point,Û D , to the Majorana hypothesis value. Our test statistic for distinguishing between the Dirac and Majorana fermion hypotheses is −2∆L, with We consider the Majorana fermion hypothesis preferred over Dirac at 3σ confidence level if −2∆L > 9.

JHEP02(2020)174
The test statistic grows as r ν and r ν deviate from 1 -it is easier to rule out the Dirac neutrino hypothesis when the expected ratio of positively charged and negatively charged leptons is more discrepant than what is expected for the Majorana neutrino hypothesis. As expected, the test statistic vanishes if r ν,ν → 1, where there is no measurable difference between the two hypotheses. So far, we have assumed that the charges of the final-state particles are always correctly identified. To take into account that this may not be the case, we introduce an efficiency factor a that corresponds to the fraction of decays that are correctly characterized in terms of particle species and charges. In practice, a = 1 corresponds to a perfect detector and a = 0.5 corresponds to pure guessing. For the Majorana hypothesis, and thus the expected numbers of events, this has no effect. In the Dirac interpretation, this will wash out any expected asymmetry. To account for this imperfection, we make the replacements in the log-likelihood (likewise for antineutrino mode). This substitution modifies eq. (7.6), which becomes We present our results in terms of contours of the test statistic for a given assumption regarding the M N prior and the efficiency parameter a, and compare against the relevant parameter space from figures 24 and 25. We focus on two cases: perfect identification (a = 1.0), corresponding roughly to the toy analysis presented in figure 29; and imperfect identification (a = 0.75). Note that the worst-case scenario is not a = 0 but a = 0.5, in which one is effectively guessing the lepton number for each signal event. Our test statistic in eq. (7.8) vanishes for a → 0.5. Figure 30 depicts the results of this analysis for the decay N → e ± π ∓ for an electroncoupled HNL. The purple curve shows our expected DUNE MPD sensitivity using the N → e ± π ∓ channel. In contrast to what was presented in figure 24, we now assume (for consistency with this analysis) that N is a Majorana fermion. 22 In various colors, we show other existing constraints in the same parameter space (discussed above). The dashed, black curve shows the 3σ contour from our analysis assuming a perfect efficiency of a = 1.0; the solid, black curve shows the same for a = 0.75. That the 3σ contours lie relatively close to the ten-event sensitivity is unsurprising. We note here that a ≈ 1 should be realistic with the DUNE MPD -there should be little confusion between identifying electron and pion tracks in the HPTPC.
Clearly, while more than ten events are required to make this distinction at moderate significance, having ∼ 40 events is sufficient; ∼ 30 decays of Dirac HNLs with L = +1 and 22 In practice, this means that we are sensitive to values of |UαN | 4 a factor of two (or |UαN | 2 a factor of √ 2) lower than in the Dirac fermion scenario. ∼ 10 decays with L = −1 can look quite distinct from the Majorana HNL expectation of 20 for each. We note that the separation between the black and purple curves shrinks as M N increases towards the kaon-production threshold near 500 MeV, owing to the increased ratio N (L = +1)/N (L = −1). Lastly, it is particularly exciting that the parameter space occupied by these 3σ contours is as of yet largely unconstrained. Not only would the DUNE MPD have sensitivity to HNLs where other experiments have not, but it could have strong resolving power between Dirac and Majorana HNLs in this hitherto unexplored parameter space. Figure 31 displays the results of this procedure for the final state N → µ ± π ∓ , assuming ten years of data collection, equal amounts of time in neutrino and antineutrino modes. The purple curve is the expected sensitivity curve for the DUNE MPD if N is a Majorana fermion (again, assumed for consistency, in contrast to the curve presented in figure 25). Other existing constraints in the same parameter space are also included. As in figure 30, the dashed (solid), black curve shows the 3σ contour from our analysis assuming an efficiency of a = 1.0 (a = 0.75). Our analysis still indicates that ∼40 HNL decays are required in order to make a meaningful distinction, but more of the corresponding parameter space is ruled out by, in particular, T2K. For heavier, D-produced HNLs, we would require efficiency of at least a = 0.75 in order to discover the nature of N . In contrast to the N → e ± π ∓ search, a ≈ 1 is significantly harder to attain because muon and pion tracks look nearly identical in the HPTPC. This result provides another reason for the need of a muon tagger for the DUNE MPD [5]. We have also explored a similar study in the channel N → µ ± ρ ∓ . In figure 25 we saw that this channel can improve on existing limits from CHARM and BEBC. However, we find that the regions of parameter space for which DUNE could discover an HNL and also determine its nature are already excluded. This is mostly because the predicted µ + ρ − to µ − ρ + ratio in the Dirac-fermion case is close to 1 (see the middle panel of figure 28 where M N ≈ 1 GeV), meaning many events are required to differentiate between the two hypotheses.
Unfortunately, the case of a τ -coupled HNL does not have a fully visible decay mode available for analysis at DUNE. Study of the τ -coupled HNL would require production of N with mass M N > m τ , which is not possible with the DUNE beam energy. If B-mesons could be produced in abundance, then it would be possible to probe the Majorana versus Dirac hypotheses in the τ -coupled case.
In principle, one could use additional kinematic information -lab-frame energies, angular distributions, etc. -to provide additional discriminating power between the Dirac and Majorana hypotheses. Specifically, the final states N → e ± π ∓ and N → µ ± π ∓ share many of the same features regarding anisotropy as the final state N → νπ 0 (see section 7.4.1) when one combines the distributions of the charge-conjugated final states [148]. This anisotropy leads to different predicted energy spectra for the lab-frame electron/muon/pion distributions. We have simulated these differences and found them to be small -such information would only assist in the analyses of figures 30 and 31 if the number of signal events is already quite large.
It is proposed that the DUNE near detector facility will include a magnetized multi-purpose detector consisting of a high-pressure gaseous argon TPC with a surrounding electromagnetic calorimeter. The components of the MPD allow for good particle identification, tracking, and momentum resolution, and the HPTPC enables measurements of particles down to low energies. The primary purpose of the MPD is to control the systematic uncertainties associated with neutrino oscillation measurements but these capabilities also lend it great power in searching for new light physics that is weakly coupled to the Standard Model. The DUNE MPD is ideally suited to search for new particles that are light enough (m < ∼ 1 GeV) to be copiously produced in meson decays and sufficiently weakly coupled that they have a long lifetime (γβτ ∼ 500 m) and can reach the MPD. In addition to the enhanced measurement properties of the MPD, the reduced target density means that the signals suffer from fewer beam-induced backgrounds than in the LArTPC.
New physics connected with neutrinos or dark matter is typically weakly coupled to the SM, and is long-lived. If dark matter is GeV-scale in mass and thermally produced, then it must sit in a more complicated dark sector that contains additional light mediator states that can also be long-lived. In this work, we have investigated the sensitivity of the MPD to sub-GeV dark mediators or heavy neutral leptons. We have focused on models where these new states are coupled to SM degrees of freedom through one of the so-called renormalizable portals: the vector, scalar, or neutrino portals of eq. (1.1). Light particles with these couplings can be produced in the decays of mesons, predominantly neutral and charged pions and kaons, and then subsequently decay in the gaseous detector, often into charged final states. In each case analyzed, we have also estimated the expected background rates and determined regions of parameter space where the new physics signal would be large enough to be seen above backgrounds.
In the case of kinetic mixing of a dark photon with the SM (section 4), we find that after ten years of data taking, the DUNE MPD will have sensitivity to presently untested regions of parameter space. Due to different target-detector separations, DUNE is sensitive to a complementary region of parameter space to other future experiments like SHiP or SeaQuest. For dark photon masses ranging from 200 MeV to 1 GeV and kinetic mixing parameter 10 −8 < ∼ < ∼ 10 −7 , the DUNE MPD will observe at least 10 decays of dark photons to e + e − , µ + µ − , or π + π − final states, and, for some parts of parameter space, as many as 100 such decays. The excellent measurement capabilities of the gaseous argon detector, coupled with the surrounding ECAL and magnet, allow the backgrounds from neutrino scattering events to be adequately suppressed.
A related model that also has a sizable decay rate into charged leptons is a leptophilic gauge boson (section 5). Unlike the case of a kinetically-mixed U(1), a leptophilic gauge boson has additional production through the leptonic decays of charged mesons. There is a complicated interplay between which pair of lepton-flavor numbers is gauged and the mass of the vector which determines the expected number of observable decays in the DUNE near detector. Over the accessible mass range, only the gauging of L µ − L τ is expected to probe a region of parameter space not presently probed by terrestrial experiments. For JHEP02(2020)174 L µ −L τ with gauge coupling between 10 −6 and 10 −4 and vector mass below 10 MeV, DUNE could see 3-10 V → e + e − events. Should a leptonic resonance be discovered at DUNE, the ability to distinguish electrons from muons/pions will give excellent model discrimination.
The scalar portal coupling between the SM Higgs boson and a SM-singlet scalar ϕ leads to scalar-Higgs mixing and allows for production of light scalars (section 6) through a loop decay of kaons. Although the kinematics of production and decay are similar to the vector models, the scalar decays to the heaviest state available and thus, again, particle identification capabilities will be a great benefit for model discrimination. For 50 MeV < ∼ M ϕ < ∼ 400 MeV, DUNE can make considerable improvements over existing bounds by searching for a resonance in the invariant mass distributions of charged leptons or pions.
The final portal we considered was heavy neutral leptons N mixing with SM neutrinos through the neutrino portal (section 7). The possible search channels for a heavy neutral lepton (HNL) include two-and three-body decays, and depend on whether N is a Dirac fermion (and thus lepton number is conserved) or a Majorana fermion. Considering HNLs with couplings to only one lepton flavor at a time, we find that the DUNE MPD has the potential to improve upon existing bounds by up to an order of magnitude, in the case of e or τ coupling, and by up to two orders of magnitude for µ coupling. The magnetic field acting on the MPD allows good charge identification and we find that, from the relative numbers of N → + π − and N → − π + , ( = e, µ) decays, it is possible to achieve 3σ discrimination between the Dirac and Majorana fermion hypotheses, after five years of data taken in neutrino mode. We predict that learning the Dirac/Majorana nature of neutrinos from kinematics alone, as is necessary for final states involving SM neutrinos, will be more challenging.
All the portal models that we have studied have signals that can be separated from backgrounds by taking advantage of a combination of final-state particle energies, invariant mass cuts, pointing and particle identification. The ability of DUNE to distinguish between model hypotheses, should an excess be found, or to rule out more parameter space, would be improved if pions and muons could be effectively distinguished. Thus, we advocate for the inclusion of a muon tagger outside of the ECAL. There are proposals to make the whole near detector system, including the MPD, movable, allowing measurements to be made both off-and on-axis. The different kinematics of signal and background events, and between new physics models, means that taking data at multiple angles should allow for signal-background separation and model discrimination. We have not investigated this potential here, leaving it for future work. Overall, the inclusion of a gaseous detector in the DUNE near detector complex will provide superb sensitivity to a broad class of new physics models, enhancing the physics program of DUNE beyond neutrino oscillation measurements.

A Meson production details and beam energy comparison
In this appendix, we expand on the simulations discussed in section 3, providing more detail regarding the meson production rates extracted from Pythia8. We also compare results for an 80 GeV proton beam with those for a 120 GeV proton beam, both of which are under consideration by the DUNE collaboration. Table 6 lists the number of mesons produced per proton-on-target under several assumptions, all using the ''SoftQCD:all'' flag within Pythia8 for the production. We provide the associated Particle ID for clarity. We perform four separate simulations; two each with 80 GeV protons and 120 GeV protons. Because the DUNE-LBNF target is expected to be graphite, we perform simulations assuming the protons are hitting either protons or neutrons, however the resulting output is similar in both cases. To obtain an expected number of mesons per proton-on-target, for our simulations, we average the pp and pn results, assuming an isoscalar target. This leads to the results in tables 4 and 5.
From table 6 we see that light meson production rates are around 10% larger for a 120 GeV beam than for an 80 GeV beam, while D/D s meson production differs, roughly, by JHEP02(2020)174 m + (P ) M ν M ℓ Figure 32. Feynman diagrams contributing to the process m + → + νV .
a factor of two. We also explored the energy distributions of outgoing mesons in these two samples and find them to be largely similar. Combining these two differences, we expect that the sensitivity to different types of new physics with an 80 GeV beam would be mildly less powerful than that with a 120 GeV beam but the differences will be minor. Finally, we note that the DUNE collaboration is also exploring a "high-energy" beam configuration, in which the Standard Model neutrino flux consists of higher energy neutrinos. This configuration is obtained by changing the focusing magnets, not the proton beam energy or target composition, meaning that table 6 would be unchanged. All results we present regarding new physics particles being produced in the decay of neutral or shortlived mesons would be unchanged. We leave studies regarding the results of high-energy beam operation with new physics particles coming from long-lived, charged mesons (mostly heavy neutral leptons coming from π ± and K ± decays -see section 7) to future work.

B Derivation of charged meson decays into leptophilic vector bosons
In section 5, many of our results for the sensitivity to the leptophilic vector bosons associated with the L α − L β current relied on production via charged meson decay, m ± → ± νV . In this appendix, we provide a derivation for the partial width of such a decay (see ref. [149] for a version of this result for the K ± → µ ± ν µ V decay channel).
In figure 32 we show the two Feynman diagrams that contribute to this decay, in which V is emitted by either the final-state neutrino ν α or the charged lepton + . We refer to the matrix elements of these two diagrams as M ν and M , respectively. These matrix elements may be expressed as where G F = 1.16 × 10 −5 GeV 2 is the Fermi constant, g V is the new vector boson's gauge coupling, F m is the decay constant of the meson m, and V qq is the relevant quark mixing matrix element for this decay. We label the initial meson momentum as P and the final state momenta as p 1 (charged lepton + ), p 2 (neutrino ν α ), and p 3 (vector boson V ). After squaring the total matrix element, we express kinematical factors in terms of invariants s ij ≡ (p i + p j ) 2 , and eliminate s 13 by using the identity s 12 + s 13 + s 23 = JHEP02(2020)174 m 2 m + m 2 + m 2 ν + M 2 V . Taking the limit m ν → 0, the matrix element squared is As observed in ref. [149], this process, like the two body decay m + → + ν, is helicitysuppressed and vanishes when m → 0. For scalar emission, that is not the case [25,28,149].