Discovering Long-lived Particles at DUNE

Long-lived particles (LLPs) arise in many theories beyond the Standard Model. These may be copiously produced from meson decays (or through their mixing with the LLP) at neutrino facilities and leave a visible decay signal in nearby neutrino detectors. We compute the expected sensitivity of the DUNE liquid argon (LAr) and gaseous argon (GAr) near detectors (ND) to light LLP decays. In doing so, we determine the expected backgrounds for both detectors, which have been largely overlooked in the literature, taking into account their angular and energy resolution. We show that searches for LLP decays into muon pairs, or into three pions, would be extremely clean. Conversely, decays into two photons would be affected by large backgrounds from neutrino interactions for both near detectors; finally, the reduced signal efficiency for e + e − pairs leads to a reduced sensitivity for ND-LAr. Our results are first presented in a model-independent way, as a function of the mass of the new state and its lifetime. We also provide detailed calculations for several phenomenological models with axion-like particles (coupled to gluons, to electroweak bosons, or to quark currents). Some of our results may also be of interest for other neutrino facilities using a similar detector technology (e.g. MicroBooNE, SBND, ICARUS, or the T2K Near Detector).

Long-lived particles (LLPs) arise in many theories beyond the Standard Model.These may be copiously produced from meson decays (or through their mixing with the LLP) at neutrino facilities and leave a visible decay signal in nearby neutrino detectors.We compute the expected sensitivity of the DUNE liquid argon (LAr) and gaseous argon (GAr) near detectors (ND) to light LLP decays.In doing so, we determine the expected backgrounds for both detectors, which have been largely overlooked in the literature, taking into account their angular and energy resolution.We show that searches for LLP decays into muon pairs, or into three pions, would be extremely clean.Conversely, decays into two photons would be affected by large backgrounds from neutrino interactions for both near detectors; finally, the reduced signal efficiency for e + e − pairs leads to a reduced sensitivity for ND-LAr.Our results are first presented in a model-independent way, as a function of the mass of the new state and its lifetime.We also provide detailed calculations for several phenomenological models with axion-like particles (coupled to gluons, to electroweak bosons, or to quark currents).Some of our results may also be of interest for other neutrino facilities using a similar detector technology (e.g.MicroBooNE, SBND, ICARUS, or the T2K Near Detector).

I. INTRODUCTION
While the evidence pointing to the existence of physics beyond the Standard Model (BSM) is overwhelming, the new physics has so far eluded its discovery in direct-detection experiments and colliders.Nevertheless, until very recently our experimental strategy was mostly focused on unveiling the existence of new states with masses at or above the electroweak scale.Interestingly, current constraints are successfully evaded by a plethora of BSM physics models containing light, feebly interacting states that offer viable solutions to most open problems in the Standard Model (SM) and avoid large corrections to the Higgs mass.Popular examples of this kind of models include those with heavy neutral leptons (HNL), which could explain the observed pattern of neutrino masses and mixing as well as the observed baryon asymmetry of the universe [1][2][3]; or models with a rich dark sector, which offer novel candidates for dark matter and may be connected to the SM through renormalizable portals at low energies (see, e.g., Refs.[4][5][6][7]).Thanks to their weak interactions with the SM, models of this sort typically include new particles that are long-lived and decay to SM states with a significant branching ratio.
The existence of light, feebly interacting, unstable states can be probed in multiple ways, from their impact on cosmological observables and astrophysical objects, to direct or indirect signals in laboratory experiments and colliders.In the case of unstable particles with masses in the O(0.1)-O (10) GeV range, searches at fixed-target experiments typically offer the best constraints (see, for example, Refs.[8,9] for recent reviews).The key in this case is that, once produced, a long-lived particle (LLP) may propagate over tens or even hundreds of meters before decaying into visible final states in nearby detectors.
In recent years, the experimental search for LLPs has received considerable attention from the neutrino community.Accelerator-based neutrino experiments are entering a precision era with the primary goal of discovering CP violation in the lepton sector, but they also offer all the necessary ingredients to conduct sensitive fixed-target searches: namely, high-intensity proton beams producing large fluxes of mesons, as well as versatile near detectors.The current generation of accelerator-based neutrino oscillation experiments set already some of the leading constraints for certain LLP models (including dark scalars [10][11][12], axion-like particles [13][14][15], or HNLs [12,[16][17][18][19], among others).These will surely be improved over the next decade by the two upcoming new-generation long-baseline neutrino oscillation experiments, Hyper-Kamiokande [20] and the Deep Underground Neutrino Experiment (DUNE) [21], currently under construction.
In this work, we highlight the unique opportunity offered in this regard by DUNE, which will be exposed to the LBNF neutrino beam [22].In addition to its very high intensity, which leads to a considerably higher flux of pions and kaons with respect to other facilities, the high proton energy available at LBNF (120 GeV) will allow the production of a significant flux of heavier resonances such as D mesons.This will enable the DUNE experiment to provide leading constraints on LLPs with masses above the kaon mass, a window that is otherwise challenging to explore for laboratory experiments.
In order to make our results as model-independent as possible, we will present them as a function of the mass of the new state and its lifetime, in line with our past works [13,[23][24][25] (see also Refs.[26][27][28] for related works that also follow a model-independent approach).However, in order to put the expected sensitivities of DUNE in context and to ease their comparison to present limits, it is useful to consider a specific model.Thus, we will also consider models with pseudoscalar particles, which arise naturally as pseudo-Nambu-Goldstone bosons of a spontaneous global symmetry breaking, and are therefore ubiquitous in extensions of the SM.These particles are often referred to as axion-like particles (ALPs), since the best-motivated example is the QCD axion [29,30].Generic models with light unstable pseudoscalars may lead to a wide set of new physics signals in neutrino detectors depending on their couplings to the SM, including ALP decays into pair of electrons, muons, photons, or multiple mesons.The DUNE experiment will offer the possibility to study many of these, thanks to its highly-capable suite of near detectors (ND) [31], which will include both a liquid argon (LAr) and a gaseous argon (GAr) time projection chamber (TPC).In particular, while many of the constraints on ALPs rely on their coupling to photons and electrons, ALP couplings to muons are subject to fewer limits.For a wide class of models, these can be particularly relevant in the mass window above 2 m µ where the decay channel a → µµ dominates.Modern facilities operating with TPC detectors offer excellent opportunities to provide leading constraints for these scenarios, as pointed out in Refs.[13,32].Moreover, as the imaging capabilities of TPCs allow the possibility of studying complex final states with multiple particles (a key advantage with respect to other detector technologies used in neutrino experiments), in this work we will also consider ALP decays leading to multiple pions (for similar studies, see Refs.[33,34]).
The multiple possibilities offered by DUNE to search for LLP decays have been pointed out in the literature before (e.g., in Refs.[26,[32][33][34][35][36][37][38][39][40][41][42][43][44]).In most of these studies the effect of the background has been neglected, arguing that it could potentially be reduced to a negligible level by means of appropriate selection cuts.However, it should be kept in mind that the near detectors available at neutrino facilities are exposed to a high-intensity neutrino flux.The whole DUNE ND suite, for instance, will register more than 100 million neutrino interactions per year [31].Even though the kinematic features of neutrino-nucleus scattering events are different than those expected from an LLP decay, reducing the expected background to a negligible level is a daunting task, as it was shown, for instance, in Refs.[41,43] for the HNL scenario.The final states expected in the case of HNL decays through mixing are, however, different from those expected for other LLPs.In this work, we compute the expected backgrounds for generic final states involving two photons, two leptons, and multiple pions, with no missing energy.Thus, our background study is, a priori, applicable to BSM extensions including light vectors, scalars, or pseudoscalars.
In principle, considering event rates from LLP-argon interactions at the ND, in addition to LLP decays, may provide additional sensitivity to certain BSM scenarios.This is specially so in regions of the parameter space corresponding to very long lifetimes, since the probability for an LLP to decay inside the detector becomes heavily suppressed.As pointed out in Refs.[35,37,45], the inclusion of scattering events for ALPs could be useful to close the so-called cosmological triangle; however, for ALPs above the MeV scale, the corresponding bounds are considerably worse than those obtained from decay searches (e.g., see Fig. 3 in Ref. [35]).In our case, the study of scattering signals would demand a more involved simulation and analysis of the relevant backgrounds, falling outside the scope of the present work.
This article is structured as follows.Section II describes the computation of the decay signal of a generic LLP.It also summarizes the main details of the ALP benchmark models considered in this work, where we assume that the ALP is coupled to the SM with different sets of effective operators.Current constraints on these benchmark models are then summarized in Section III (as well as in Appendix A).Next, we discuss our evaluation of the backgrounds in Section IV, before presenting our results in Section V. We summarize and conclude in Section VI.

II. PRODUCTION AND DECAY OF LONG-LIVED PARTICLES
The computation of the number of decays of long-lived particles (LLP) in detectors located at a certain distance can be carried out as follows.We start from a Monte Carlo simulation of the production of parent mesons (technical details regarding these simulations are provided in Sec.IV).The number of mesons can be written in terms of the number of protons on target (PoT) and the meson yield per PoT (Y M ) as where d 2 ρ M /(dE M dΩ M ) is the probability that a given meson is produced with energy in the interval [E M , E M + dE M ] and with a trajectory defined by a solid angle in [Ω M , Ω M + dΩ M ].Using Eq. ( 1), the expected number of particles a produced in the decay M → a + . .., with a branching ratio BR(M → a), can be computed as where d 2 ρ M →a /(dE a dΩ a ) stands for the differential probability that a meson will produce an LLP a with energy and trajectory defined by E a and Ω a .As we will see, in certain scenarios the LLP may also be produced directly through their mixing with SM mesons.In this case, the LLP differential flux can be approximated by the meson flux in Eq. ( 1), rescaled with the corresponding mixing angle accordingly (which depends on the masses of the two particles involved).
At this point, it is important to note that the meson fluxes depend on the production point (inside the target and in the decay volume region).Similarly, the detector acceptance depends on the production point of the LLP, since this determines whether its trajectory crosses the detector.In what follows, in order to simplify our notation, we remove the dependence with the meson and the LLP production points; we stress, however, that in our numerical calculations this dependence has been fully accounted for.
Once it has been produced, the LLP may live long enough to propagate to the detector and decay inside, with a probability that depends on its decay length boosted to the lab frame: L a = c τ a γ a β a , where τ a is the lifetime of the particle at rest, while β a and γ a are the boost factors.Such probability reads: where ℓ det is the propagation distance before the particle enters the detector, and ∆ℓ det is the length of the intersection between the trajectory of the particle and the detector (which in most cases approximately coincides with the detector length along the beam axis).Note that both quantities depend on the solid angle Ω a , as well as on the production point of the LLP.Eventually, the total number of LLP decays inside the detector into a given decay channel ch is obtained after integration over the LLP variables, and multiplying by the corresponding branching ratio and the detector efficiency for that channel, ϵ ch : where the produced LLP flux and the decay probability are given in Eqs. ( 2) and (3), respectively, and the integral in solid angle is performed taking into account only those trajectories within the angular acceptance of the detector, Ω det .Even though Eq. ( 4) is exact (and it is, indeed, what we use in our computation of the number of signal events), it might not be very illuminating.We can derive a simpler expression noting that for the experimental setup considered here: (i) the distance to the detector is much longer than the size of the target (where most of the mesons are produced); (ii) the detector size is much smaller than the distance traveled by the LLP before reaching the detector; (iii) the angular acceptance of the detector is small, which mostly selects particles traveling along the beam axis.These allow us to assume that all LLPs are approximately produced at the same point and to neglect the dependence of ℓ det and ∆ℓ det with the trajectory of the LLP.Under these approximations, we obtain: where dε M det /dE a is the angular acceptance of the detector for an LLP with energy E a .Since the detector acceptance relies on the boost of the LLP in the lab frame, it will depend on the energy and momentum distribution of the parent meson M , on the mass of the LLP, and on whether the LLP is produced in a two-body or a three-body decay.As illustration, the differential detector acceptance is shown in Fig. 1 as a function of the energy of the LLP, for K and D two-body decays and for two representative values of the mass of the LLP.As can be seen from this figure, the accepted flux ranges from a few GeV to tens of GeV, with some dependence on the parent meson and the LLP mass.
Up to this point, the discussion is model-independent and applies to any LLP produced from meson decays at a neutrino beamline.The final sensitivity will be fully determined by the mass of the particle and its lifetime, and it will be optimal for values of cτ a which maximize the decay probability in Eq. (3).In order to compute the expected number of events, one just needs to know the parent meson and whether the LLP is produced through a two-body or a three-body decay.Nevertheless, in a given model the production branching ratio and the decay width will typically be partially correlated and depend on a set of common parameters, which will also enter the decay branching ratio into a given channel.The key requirement for these searches to succeed is that the lifetime of the LLP should be long enough so that it reaches the detector before decaying.This pushes into the weakly interacting limit, which demands small couplings and, therefore, typically small production branching ratios.Consequently, in order to reach the maximum sensitivity region, a trade-off must be found in terms of the production branching ratios and lifetime, which is not always easy to do and depends severely on the model.The second aspect to consider is that specific benchmarks will translate into preferred decay channels, which ultimately determine the experimental strategy to follow in order to maximize the signal-to-background ratio.
In the rest of this section, we will focus on three benchmark models in the context of ALPs.This will allow us to write specific expressions, in terms of the model parameters, for the ALP production and decay branching ratios, as well as its lifetime.As we will see, the three models considered would lead to very distinct phenomenology in an experimental search.

A. Benchmark Model I: Gluon Dominance scenario
We start by considering an anomalous coupling between ALPs and the gluon field, leading to a higher-dimensional effective operator.This has been referred to in the literature as the gluon dominance scenario [46], and is one of the main targets for the vast suite of experiments searching for long-lived particles (see, e.g., Refs.[9,47,48]).The sensitivity of DUNE to ALPs in this scenario has been studied previously in Refs.[33,34] assuming that backgrounds could be reduced to a negligible level for the most relevant decay channels.Here, we revisit their limits taking into consideration the expected background levels, and slight improvements (outlined below) on the calculation of the ALP decay widths.Following the same normalization convention as in Ref. [46], the relevant ALP interaction Lagrangian reads: where G b µν is the gluon field strength, G bµν ≡ 1 2 ϵ µνρσ G bρσ , with ϵ 0123 = 1.Also, α s ≡ g 2 s /(4π), and g s stands for the strong coupling constant.
In this scenario, the ALPs would be mostly produced directly through their mixing with neutral pseudoscalar mesons (π 0 , η, η ′ ) below 1 GeV, and from gluon fusion for masses above this value.Here, we use the calculation from Ref. [33].As for their decays, for m a < 3m π the only available decay channel is a → γγ, since CP conservation forbids the decay into two pions, whereas for higher masses, hadronic decay modes (with three or more mesons in the final state) will dominate.ALP interactions with pseudoscalar mesons may be described using chiral perturbation theory including ALPs (aχPT), but a proper treatment of the interactions between multiple mesons requires the inclusion of Vector Meson Dominance (VMD) terms in the Lagrangian [49].Throughout this work, we follow Refs.[50,51], where a data-driven method (similar to the one developed in Ref. [52]) is used to successfully describe ALP interactions within this framework up to relatively high ALP masses (m a ≲ 3 GeV).The effective aχPT Lagrangian is then matched onto the perturbative QCD (pQCD) Lagrangian. Figure 2 shows the decay widths and branching ratios for the most relevant decay channels in this scenario, computed following Ref.[51] (see also Ref. [50]).The sharp features observed in the decay width stem from the large mixing between the ALP and the neutral pseudoscalar mesons whenever m a ∼ m π 0 , m η , m η ′ , indicated by the dotted vertical lines.

B. Benchmark Model II: ALPs coupled through electroweak operators
Next, we consider an effective Lagrangian that includes ALP couplings to electroweak operators through effective operators arising at a high-energy scale Λ (which we set to Λ = f a = 1 TeV): where B and W I stand for the EW vector bosons, ϕ is the Higgs doublet and It is worth mentioning that O ϕ can be traded by a set of (flavorconserving) fermionic operators by means of a hypercharge rotation, as shown in Ref. [53].This finally leads to: where sum over f extends to quarks and charged leptons, with c f f = c ϕ for down-type quarks and charged leptons, and c f f = −c ϕ for up-type quarks.
Main ALP decay widths (left) and branching ratios (right) for the gluon dominance scenario in Eq. ( 6), as a function of its mass.These have been computed following Ref.[51] (see also Ref. [50]).
Although the operators included in Eq. ( 8) are flavor conserving, the O ϕ and O W operators can induce flavor-changing neutral current (FCNC) processes at one loop, as detailed, for instance, in Refs.[54][55][56].For this scenario, the mechanism of production we consider will be kaon decays to ALPs (K → πa), in line with our past work in Ref. [13].The corresponding width can be computed using the aχPT Lagrangian [53,57]: where f 0 is the scalar form factor1 and In Ref. [13] we computed the running of the couplings from Λ = 1 TeV down to µ = 2 GeV, following Ref.[59].This leads to the following matching condition for the effective coupling entering the decay width: We note that a similar effective coupling is generated in this class of models for the decay B → Ka; however, the production of B mesons is insufficient at DUNE and will not be considered here.ALPs produced from kaon decays will have a kinematical threshold at m a < m K − m π ∼ 355 MeV.In this mass window, the ALP can only decay into photons or light charged leptons pairs (e and µ).The decay width for the di-photon decay channel is given by where the effective coupling at low scales is given at one loop by [55,60] Here, we have written c i ≡ c i (Λ), B 0 and B 2 are loop functions (which can be found, for example, in Appendix B of Ref. [13]), τ W = 4m 2 W /m 2 a and α is the fine-structure constant.Finally, the decay width into dilepton pairs is given by where c ℓℓ has been computed at one loop, and at low energies (µ ∼ 2 GeV) reads [55,60] By direct observation of Eqs. ( 11)-( 15), we note that: • Γ(a → ℓ + ℓ − ) is suppressed by a factor of m 2 ℓ /m 2 a with respect to Γ(a → γγ).Therefore, for similar values of c γγ and c ℓℓ , we expect Γ(a → γγ) ≫ Γ(a → ℓ + ℓ − ).
• For ALPs masses m a > 2m µ , the decay channel a → µ + µ − will completely dominate over the decay channel a → e + e − , since Γ(a • As can be seen in Eq. (11), O B has a subdominant effect in the ALP production.However, it can have a significant effect in the ALP decay width, affecting its lifetime.For simplicity, hereafter we set the coefficient of this operator to zero; however, we refer the interested reader to Refs.[13,55] for a related discussion on this issue.
The corresponding decay widths are shown in Fig. 3 for an ALP coupled predominantly through O W (left panel) and O ϕ (right panel).Finally, we conclude this subsection stressing that the production of ALPs from decays of heavier mesons (e.g., D + → π + a) using couplings to electroweak operators would require to generate the effective flavor-violating coupling [k Q ] cu via loop.However, proceeding analogously as for the kaon case, it is easy to see that the corresponding coupling for D decays will be severely suppressed, since Therefore, no sensitivity to these operators is expected at DUNE for m a > m K .

C. Benchmark Model III: Charming ALPs
Here, we consider a different set of operators that could induce ALP production from D decays, which would allow to probe even heavier ALPs.This is the case, for example, of effective operators coupling the ALP to a quark current with off-diagonal couplings in flavor space (see, e.g., the discussion in Ref. [61] and references therein).The inclusion of couplings to quarks m a (GeV) implies, however, that the ALP will have a significant decay width into hadronic final states.Consequently, its lifetime and branching ratios will show a significant dependence on the assumed Wilson coefficients and it becomes necessary to choose a particular texture in flavor space.As an illustrative example, we will consider ALPs that couple to right-handed up-quark currents as in Ref. [62], dubbed as charming ALPs: where c u R is a Hermitian matrix in flavor space, and i, j = 1, 2, 3 are quark family indices.The flavor-violating coupling [c u R ] 12 in Eq. ( 16) induces D + → π + a, with a decay width given by [62,63] Γ(D where λ is defined in Eq. ( 10) and f Dπ 0 m 2 a is the scalar form factor evaluated at a momentum transfer q 2 = m 2 a , which we take from the lattice computation in Ref. [64].For ALPs produced in D-meson decays, we can probe a wide mass window up to m a = m D − m π ≲ 1.7 GeV.As in the gluon dominance scenario, in this case the dominant decay modes will involve multiple mesons in the final state, which requires the use of the aχPT for masses below approximately 2.5-3 GeV.We note that the calculation in Ref. [51] describes ALP interactions with mesons induced by purely diagonal couplings in flavor space.However, we can use it for the purposes of this work since the only off-diagonal couplings considered involve c and t quarks, which are irrelevant for decays in this mass region.In particular, we stress that as long as the ALP couples to up-quarks only, its decay widths are mostly determined by the [c u R ] 11 coupling in Eq. ( 16), while the remaining couplings only enter at subleading order for some channels (such as a → γγ).
In summary, for the purposes of this work, the phenomenology of this class of models is fully determined by [c u R ] 12 (which controls the production in D decays) and by [c u R ] 11 (which controls the decay of the ALP).Hence, it becomes interesting to consider the expected size of these two couplings for well-motivated UV completions leading to the set of operators in Eq. ( 16).The authors of Ref. [62] studied several possible UV completions of this scenario.Here, we highlight two of them, where the ALP emerges as a pseudo-Nambu-Goldstone boson (pNGB) of a spontaneously-broken symmetry: • The first one is a Froggatt-Nielsen (FN) model, where the ALP corresponds to what is usually called a flavon or a familion [65][66][67][68][69].The model includes a complex scalar field S, whose radial component is identified with the ALP, and a new U (1) flavor symmetry.
The inclusion of higher-dimensional operators at a high-energy scale Λ > f a generates the operators in Eq. ( 16) once the scalar acquires a vacuum expectation value, ⟨S⟩ = f a .As shown in Ref. [62], an appropriate choice of the charges under the new symmetry is able to generate the correct up-quark masses (thus alleviating the flavor problem) and would lead to where off-diagonal entries are controlled by ϵ = f a /Λ ∼ m c /m t .
• The second possibility is a dark-QCD (dQCD) model with a confining dark sector that contains n d dark flavors transforming under SU (n d ).A dark confining sector would offer plausible candidates for dark matter among the dark hadrons and therefore is also wellmotivated from the theoretical point of view.In addition, the ALP may be identified with one of the CP-odd states in the dark meson spectrum (a dark pion, see e.g.Ref. [51]).The operators in Eq. ( 16) may be generated, for example, through a Yukawa interaction between SM quarks and dark quarks with a heavy scalar field, which is integrated out of the theory at low energies.Adopting the same assumptions and parameter values as in Ref. [62], this leads to a texture for c u R in the u − c sector that is very similar to the one in Eq. ( 18), up to an overall rescaling factor.
In summary, both of these UV completions would lead to [c u R ] 12 /[c u R ] 11 ∼ O(0.01)-O(0.03),with small differences depending on the particular scenario being considered.Fixing a texture is convenient since it allows to compute the corresponding ALP decay widths for all relevant decay channels, and sets the correlation between the ALP production rates and its lifetime.The corresponding branching ratios for the texture in Eq. ( 18) are shown in Fig. 4, where we have set the matching scale between the aχPT and the pQCD Lagrangians at 2.9 GeV following Ref.[51].We see that for m a < 3m π the ALP decays exclusively to a → γγ.Conversely for m a > 3m π the decay a → π + π − π 0 rapidly dominates, as shown in Fig. 4, with the exception of the regions neighbouring the mass values of m η and m ′ η .The resulting branching ratios for a dQCD-inspired model would be very similar to the ones shown here for the FN-inspired model, although the corresponding ALP lifetime would be longer.This would affect the statistics expected at the detector (inducing a change in sensitivity), but the results would otherwise be qualitatively very similar.Therefore in what follows for concreteness we will adopt the texture in Eq. ( 18) to show our results for this benchmark scenario.
We finalize this section by stressing that, while we have based our discussion on theoretically well-motivated examples, our results for this benchmark scenario will apply to a wider class of models, as long as the ALP does not couple to down quarks directly2 and both  O u R FIG. 4. Main decay widths (left) and branching ratios (right) for the Charming ALP scenario, as a function of its mass.These have been computed following Ref.[51], for the couplings in Eq. ( 18).

III. PREVIOUS CONSTRAINTS
This section contains a brief summary of the most relevant constraints for the two benchmark models under consideration.Most of the bounds from current and past experiments have been previously computed in the literature, and are just recasted here for the scenarios of interest.In the case of an ALP coupled through electroweak effective operators, a comprehensive discussion on the applicable bounds can be found e.g., in Refs.[48,63,70] (see also Ref. [13]).For an ALP coupled to right-handed up-quarks, we follow the discussion in Ref. [62], where the authors derived the most relevant bounds on this scenario.However in this case we have rederived some of the constraints here, finding significant differences as explained in more detail below.
A. Visible decay searches.
At fixed-target experiments, the ALP can be produced through various processes, depending on its couplings to SM particles: from meson decays, through Primakoff scattering, through its mixing with the neutral pseudoscalar mesons, or via proton or electron Bremsstrahlung.Here we distinguish three different cases: • For an ALP coupled to gluons through O G in Eq. ( 6), constraints of this type are obtained from K ± → π ± γγ measurements at NA62 [71], E949 [72], K L → π 0 γγ at NA48/2 [73] and KTeV [74], and from B → Ka, a → γγ searches at BaBar [75].We take these limits from Ref. [48].Additional limits come from searches for LLP decays into two photons at CHARM [76] (derived in Ref. [34]).
• For an ALP coupled through the O W operator in Eq. 8, significant constraints are obtained from E137 searches [77] (see e.g.Refs.[78][79][80]), which we take from Ref. [78], as well as from K → πγγ (taken from Ref. [48]) and from B → Ka(a → γγ) [75].Additional bounds can be obtained from NA64, for an ALP produced in the forward direction through the Primakoff effect in the vicinity of a nucleus, recasting their search for ALPs decaying into two photons [81].Finally, ALPs could be produced at LEP through an off-shell photon (for example, via e + e − → γ * → γa) or in photon fusion (e + e − → e + e − a) and decay to two photons.The corresponding bounds are taken from Ref. [82].
• For an ALP coupled through the O ϕ operator in Eq. ( 8), relevant bounds are obtained recasting CHARM bounds on LLPs decaying into lepton pairs [76,83].Here we use the revised bounds obtained in Ref. [84] (see also Refs.[80,85]).Significant constraints are also obtained from LHCb [86,87], for ALPs produced from B-meson decays (via B → Ka) and decaying within the detector into µ + µ − .We take these from Refs.[55,84].Finally, in our previous work [13] we obtained new bounds from a recast of a MicroBooNE search for e + e − pairs from a long-lived particle pointing to the NuMI absorber [11].Here we also add to these the corresponding recast for a similar search into µ + µ − pairs [12], following the same methodology as in Ref. [13].We also note that the ArgoNeuT experiment has recently obtained a bound on axion-like particles for heavier masses (up to 700 MeV) decaying into µ + µ − pairs [15].These cannot however be easily recasted to the scenario considered here without a proper simulation of the meson fluxes in the NuMI target (in particular, from η meson production), which lies beyond the scope of this work.
• For an ALP coupled to right-handed quarks, Eq. ( 16), relevant constraints can be obtained for CHARM [76], using the null results from a search for a → γγ as outlined in Ref. [62].
We rederive such constraint here, finding significant differences which we attribute to the different treatment of the decay width of the ALP as well as to a different simulation of the D meson production.Specifically, we obtain the D-meson fluxes from Pythia (v8.3.07)[88] for a beam energy of 400 GeV, and compute the signal acceptance of the detector using the same methodology described in Sec.II, for a detector with a decay volume of 3 × 3 × 35 m 3 located at L det = 480 m from the target.Our computation of the lifetime of the ALP, as well as the branching ratio for the a → γγ decay, follows Ref. [51].
Finally, we note that data from atmospheric neutrino oscillation experiments may provide relevant constraints as well.For example, in Refs.[23,24] a model-independent analysis of SK multi-ring data was able to set a constraint on HNL production BR(K → N ) × BR(N → e − like) ≲ 5 × 10 −9 for values of cτ ∼ O(km).A similar analysis would also constrain the scenarios considered here, for ALP decays into γγ or e + e − .

B. Invisible decay searches.
For very light masses, or for sufficiently small couplings, the ALPs may exit the detector without decaying (hence the term "invisible decay").Results for K → ν ν or B → Kν ν may then be reinterpreted as bounds on the branching ratios for K → πa and B → Ka [54,55,80].For heavier masses, while no dedicated search for D → πa exists, an indirect constraint may be obtained from the reinterpretation of D + → τ + ν, τ + → π + ν measurements, as pointed out in Refs.[61,62].
The strongest limits on K → πa are obtained from the NA62 experiment [89,90].However, the very competitive bounds from E787 & E949 [91] are comparable (and even dominate) in the region close to the pion mass.For B → Ka decays, we use the constraint from Belle on B + → K + ν ν [92].Experimental limits on invisible ALP decays also arise from precision measurements of the pion momentum in K → πX [93].Here we draw from our previous work in Ref. [13] where we derived the corresponding constraints for an ALP coupled through the electroweak operators in Eq. (8).For the gluon dominance scenario, we extract the relevant bounds from Ref. [48].A priori, similar bounds may also be derived for the Charming ALP scenario.However it can be seen that the leading contribution is obtained when the c-quark is running in the loop.In particular, following Refs.[13,59], we find Using this effective coupling, we re-derive here the corresponding constraints from K → πa, following the procedure described in Ref. [13]; however, while in Ref. [13] the running was taken from Λ = 1 TeV, here we take into account the relationship between f a and Λ, see Eq. (18).In other words, for a given value of f a the running of the coupling constants is carried out from Λ = f a /ϵ down to µ = 2 GeV.
For the charming ALP scenario, bounds on D → πa can be derived using CLEO [94] or BESIII [95] measurements of the D → τ ν decay.Their data is provided as a function of the missing mass squared in both cases, which can be used to derive a limit on BR(D → πa) as a function of m 2 a , as was done in Ref. [62] for CLEO.Here we decide to use BESIII data instead, as for CLEO we are not able to reproduce their SM fit with the information provided in Ref. [94].Specifically, we take the BESIII data from Fig. 3 in Ref. [95].We have checked that we are able to reproduce their result for the determination of D → τ ν in the SM within a reasonable degree of accuracy.Our best-fit curve leads to χ 2 min /dof ∼ 29/28, indicating a good compatibility with the data with a p-value of 0.4.We then use the fit to derive a constraint on the BR(D → πa) as a function of the mass of the ALP, see App.A for details.
Finally we note that additional bounds arise from searches for mono-photon signals at colliders [54,96,97].The best limits of this class are obtained using BaBar [98,99] or LEP [100] data.Nevertheless, the resulting bounds are milder than the rest of the limits considered in this work and will not be shown here.
C. Bounds from D − D mixing.
Neutral meson particle-antiparticle pairs are allowed by the SM interactions to both oscillate and decay to lighter particles.Their oscillations are parametrized by a 2 × 2 unitary mixing matrix M .The mixing for D 0 − D0 is very suppressed in the SM, leaving room to probe new physics that might induce such mixing.These bounds are relevant for the charming ALP scenario, since the non-vanishing [c u R ] 21 coupling in Eq. ( 18) leads to non-standard contributions in the form of effective operators at low energies involving four quarks.These contribute to the off-diagonal entry in the mixing matrix, M 12 .Here we take the bounds derived in Ref. [62] using the mixing constraints on M 12 from Ref. [101].

D. Astrophysical bounds.
Within the mass range of interest here, the main constraints are derived from supernovae data.These can be further classified into three subcategories: (1) constraints obtained from the requirement that the energy loss induced by ALP emission does not exceed that of neutrino emission [102][103][104][105][106]; (2) bounds coming from the visible photon signal resulting from the ALP burst [107,108]; and (3) limits obtained from the observation of low-luminosity core-collapse supernovae, which constrain the total energy deposition in the progenitor star from radiative ALP decays [109].Additional limits may be derived from X-ray observations of neutron star mergers associated to gravitational waves [110].
Here we take the limits on c W from Ref. [48], while for the charming ALP scenario we take the limits from Ref. [62].Finally, while bounds derived on ALPs coupled to muons or electrons do exist in the literature [104,109,[111][112][113][114], to the best of our knowledge none of these analyses may be easily recasted for our O ϕ operator (which induces couplings to electrons, muons and, via loop, to photons) within our mass range of interest.Consequently, no bounds from supernovae are included in this case.We stress, nonetheless, that astrophysical bounds are typically relevant for couplings much smaller than those within the expected sensitivity reach for DUNE, and are only shown here for completeness.
IV. TARGETING SIGNALS OF LLP DECAYS AT THE DUNE NEAR DETECTORS DUNE [21] is an upcoming long-baseline neutrino oscillation experiment under construction in the United States.Once built, DUNE will consist of two state-of-the-art neutrino detectors exposed to the world's most intense neutrino beam.The so-called near detector (ND) will record neutrino interactions near the source of the beam, at the Fermi National Accelerator Laboratory (Fermilab), in Illinois.The much larger far detector (FD) will be built at a depth of 1.5 km at the Sanford Underground Research Facility (SURF), in South Dakota, 1300 km away from Fermilab.DUNE has a very rich scientific program that includes the precision study of neutrino mixing [115,116], astroparticle physics [117] and searches for BSM phenomena [118].
In this section, we discuss the prospects for detecting the decay of LLPs at the DUNE ND [31], which will be built in a shallow underground hall located 574 m downstream from the neutrino beam origin.In the first phase of DUNE, the ND will consist of a liquid-argon Time Projection Chamber (TPC), ND-LAr, followed by a downstream muon spectrometer (the so-called temporary muon spectrometer, or TMS).In the so-called Phase II of the experiment, TMS will be replaced with ND-GAr, a magnetized, high-pressure gaseous argon TPC surrounded by a calorimeter.We do not consider in our study the magnetized beam monitor known as SAND, which will be present in both phases of the experiment.
The search for LLP decays at the DUNE ND will suffer from a significant background from neutrino interactions, given the intensity of the LBNF neutrino beam.Considering this, ND-GAr appears to be the ideal detector for this kind of searches, given its large volume (proportional to the number of signal events) and low mass (proportional to the number of background events).However, ND-GAr will only be available as part of the upgrades contemplated for Phase II [119].Therefore, in this work we have also estimated the expected sensitivities for the ND-LAr detector, as it will start collecting data much sooner.Optimal selection cuts for ALP signals at ND-LAr and ND-GAr are generally different, due to their different features.Thus, our background analysis is performed for each detector separately.

A. Simulation
We have used large simulation data sets and a basic event selection to estimate the signal efficiency and background rejection in the DUNE ND for the main LLP decay channels of our benchmark models.Each simulation event represents a single LLP decay or neutrino-argon interaction in the active volume of the detectors, thus ignoring possible pile-up effects (i.e. the cross-contamination of different interactions occurring in the same TPC event) since the ND design will be optimized to make them negligible [31].For simplicity, detector effects are not simulated, but we do take them into account in our study with the introduction of the typical detection thresholds and resolutions expected from ND-LAr and ND-GAr.
The first step in our simulation involves the generation of the LLP fluxes at the ND starting from the decays of their parent mesons.For kaon decays, we make use of publicly-available histogram files [120,121] that contain the position and three-momentum distributions of the decay of mesons in the LBNF beamline as obtained with G4LBNF [120], the official Geant4 simulation of the LBNF beamline from primary proton beam to hadron absorber.This simulation code includes kaon production in the target, along the 194 m-long decay pipe, and in the absorber at the end of the decay pipe, leading to a kaon production yield Y K = 0.75 mesons/PoT.For the D mesons, we use the same distributions as in Ref. [42].They were obtained using Pythia (version 8.2.44) [122] to create a pool of events for proton collisions at various momenta, followed by a Geant4 simulation to predict proton inelastic interactions between the 120 GeV primary proton beam and the target.Doing this, the D-meson production yield is Y D = 1.2 • 10 −5 mesons/PoT.As for the luminosity considered, we take a total nominal exposure of 1.1 × 10 22 PoT, corresponding to 10 years of DUNE operation [115].
The parent mesons will decay to LLPs (as detailed in Sec.II), which will then propagate towards the DUNE ND, located 574 m away from the target.In our computation, we approximate the ND-LAr as a rectangular cuboid 7 m wide, 3 m high and 5 m deep, with an assumed fiducial volume that excludes 50 cm from the sides and upstream end and 150 cm from the downstream end [31].The active volume of the HPgTPC is a cylinder with a radius of 260 cm and a length of 500 cm; for the purpose of computing event rates, we define a fiducial volume by excluding the outer 37 cm of radius and 30 cm on each end of the cylinder [31].
We assume that the dominant background source in this search will be neutrino-argon interactions in the active volume of the TPCs.We estimate that other possible background sources (such as neutrino-rock interactions or cosmic muons) will be negligible in comparison, as the resulting events will not be aligned with the direction of the beam, in general.Using the GENIE neutrino Monte Carlo generator (version 3.2.0)[123] and the public DUNE flux histogram files [120,121], we have produced 2 × 10 7 ν µ -Ar interactions.
Below, we discuss the event selection we have devised for the four main ALP decay channels relevant for the models discussed in Sec.II.Table I summarizes our results.We have verified that our estimates do not change significantly for different LLP masses (m a ) in the relevant range for each channel.In this regard, it is worth noting that the reconstructed invariant mass of the LLP would provide a handle for the discrimination of signal and background that we have not exploited in our study.

B. Event selection: µ + µ − decay channel
A priori, the most important background source for the di-muon decay channel is ν µ chargedcurrent events with charged pions, as it is relatively easy to confuse muon and pion tracks due to their similar stopping power (dE/dx) in argon.About 38% of the ν µ -Ar interactions have a charged muon and a charged pion above threshold in the final state, with an expected rate of the order of 6×10 5 events per ton-year at the DUNE ND.Actual di-muon events from charged-current charm production only represent less than one percent of the total background events.
We start our event selection requiring candidate signal events to have only two µ-like tracks (i.e., µ ± or π ± ) above threshold.This allows the rejection of background events with hadronic activity near the interaction vertex.We consider a proton detection threshold of 40 MeV for the ND-LAr and of 5 MeV for the ND-GAr [31], obtaining similar rejection factors for both detectors: about 0.3% of the initial events meet the above criterion.From this point on, each detector requires slightly different considerations, discussed next.
In ND-GAr, the TPC combined with the electromagnetic calorimeter (ECAL) and the muon identification system that surround it will provide superb µ/π separation capabilities, reaching 100% purity in the identification of muons for a wide range of momenta [31].Moreover, the magnetic field in ND-GAr will allow the measurement of the charge sign of muons.Thanks to these capabilities, we can reduce the initial sample of background events by a factor 6 × 10 −6 for, essentially, perfect signal efficiency.Lastly, we can further reduce the background sample taking into account the particular kinematics of signal events (see Fig. 5, top row): (i) The reconstructed LLP transverse momentum should be low; that is, the LLP trajectory points back in the direction TABLE I. Signal efficiencies and background event rates for the different decay channels, before and after event selection according to the cuts discussed in the main text.Results are shown separately for the two DUNE near detectors considered.Background event rates are provided per year, and for the total fiducial volume considered for each detector.We highlight in bold type the large backgrounds expected for some of the decay channels, as well as the reduced LAr ND signal efficiencies for most decay channels considered.

Selection cut
Signal of the target.(ii) The muons in signal events are highly boosted, and thus the angle between them should be small.We assume that a momentum resolution of the order of 5% or better and angular resolution of the order of a few degrees can be achieved in ND-GAr for momenta up to 10 GeV/c [31].Overall, these cuts let us achieve a background rejection in excess of 10 7 (see Tab. I), resulting in a background-free search in 10 years of data taking.
In the case of ND-LAr, the detector will not be able to fully contain high-energy muons or measure lepton charge, but the downstream spectrometers (TMS in the first phase of DUNE, and ND-GAr in the second one) will measure the charge sign and three-momentum of the muons that enter them.Events with muon kinetic energies below 1 GeV will be contained within ND-LAr, while events with higher energy muons traveling within 20 degrees of the beam direction will exit ND-LAr and enter the spectrometer [31].TMS will only be able to measure muons up to ∼ 6 GeV/c before they range out, corresponding to 40% of our LLP decays. 3We will assume as well that the combination of the dE/dx measurement in ND-LAr plus the µ/π separation capabilities of the TMS -pions will interact inelastically in the steel layers of TMS with high probability, while muons will behave as minimum ionizing particle -will be enough to reach essentially perfect purity in the identification of muons, such as in the ND-GAr.Finally, the two kinematical cuts described above are applied, achieving a background-free search in 10 years of data taking.As a point of comparison, the analysis described in Ref. [124] for the identification of di-muon neutrino trident events achieved a background suppresion of 6 orders of magnitude using exclusively kinematical cuts in ND-LAr.All histograms are area normalized.Top row: transverse momentum (left panel) and angle between tracks (right) for µ + µ − candidate ALP-decay events in the DUNE ND-GAr.Central row, left: angle with respect to the beam direction of e + e − events.Central row, right: angle between the two photon showers for γγ candidate events in ND-LAr.Bottom row: transverse momentum (left panel) and angle between the charged-pion tracks of π + π − π 0 candidate events.

C. Event selection: e + e − decay channel
The dominant background source for this decay channel is ν µ neutral-current single-pion (NCπ 0 ) events.They can be mistaken as an e + e − signal if only one of the photons from the π 0 → γγ decay -with a branching ratio close to 99% [125] -converts in the detector, and no visible hadronic activity occurs at the vertex of the neutrino interaction.Less than 1% of the ν µ -Ar interactions result in a single π 0 , corresponding to an expected rate of about 1.5 × 10 4 events per ton-year at the DUNE ND.
In LAr, the attenuation length for gamma rays in the energy range from 100 MeV to 10 GeV is of the order of 20 cm [126].This limits considerably the probability that only one of the two photons from a π 0 decay interacts in ND-LAr.Through Monte Carlo simulation, we have estimated that probability to be 1.3%.Conversely, in high-pressure argon gas, the attenuation length for gammas is well above 10 m.Therefore, in most cases, both gammas from the π 0 decay will escape the TPC, interacting in the ECAL; only in about 1% of the events the decay will result in a single-gamma conversion in the GAr, with no associated activity in the ECAL [31,36].
Lastly, pair-conversion events in both ND-LAr and ND-GAr can be suppressed requiring the reconstructed direction of the LLP to be aligned with the neutrino beam (see Fig. 5, middle left panel).Background events, in contrast, will follow a nearly isotropical distribution, as the neutral pions will be coming from neutrino-nucleus interactions.A rather conservative cut of 100 mrad (about 5 • ) for both ND-GAr and ND-LAr suppresses the background events by an order of magnitude.
For ND-LAr, as far as we are aware, there is no estimation yet of the reconstruction efficiency of high-energy electron showers.Therefore, taking into account the significant depth of the showers (more than 2 m for an electron of 1 GeV) and the busy DUNE ND environment, we conservatively assume an efficiency of 10%.

D. Event selection: γγ decay channel
The most important background for this channel is, like in the previous case, the decay of neutral pions resulting from ν µ neutral-current interactions.
In the case of ND-GAr, we assume that these events can only be reconstructed with good efficiency and high purity if both gamma rays are detected in the ECAL [31].Nonetheless, this calorimeter will not be able to contain completely the electromagnetic showers from high-energy photons up to 10 GeV, as it will only be 8-12 radiation-lengths deep [31].The impact of this shower leakage on the energy and angular resolutions can only be fully understood with a detailed simulation, beyond the scope of this paper.We will only assume here that background events can be suppressed taking into account that signal events should point back in the direction of the neutrino beam.As the pointing accuracy of photon showers reconstructed in the ND-GAr ECAL is low (the resolution for photons of 1 GeV is about 10 • ), a selection cut on the opening angle between the two photons has limited impact and, hence, we do not consider it.
In ND-LAr, full reconstruction of both gammas requires the containment of a significant fraction of the two electromagnetic showers within the fiducial volume of the detector.Following Ref. [127], we estimate the reconstruction efficiency to be 5% in our energy range of interest.As ND-LAr should be able to reconstruct with accuracy of a few degrees the direction of the showers, the kinematical cut on the LLP direction results more effective than in the previous case, and we can also discriminate between signal and background applying a cut on the angle between the two photons (see Fig. 5, middle right panel).
In this channel, the invariant mass of the di-photon system could be used effectively for background rejection, as the distribution for background events would peak around the π 0 mass (conversely, this limits as well the sensitivity to LLPs with mass in that energy region).As this selection cut requires a detailed understanding of the invariant mass resolution of the detectors, we are not applying it in this analysis.
The most important background source for this decay channel is ν µ neutral-current events with multi-pion production, with an expected rate of the order of 39000 events per ton-year.The same arguments given in previous sections for the reconstruction of µ/π tracks and the two photon showers from the π 0 decay apply here.Background events can be suppressed with selection cuts on the transverse momentum and angle between the two charged-pion tracks (see Fig. 5, bottom row).

V. RESULTS
To compute the expected sensitivity to the Wilson coefficients in Eq. ( 8), or to the value of f a in Eqs. ( 6) and ( 16), we perform, for each detection channel, an unbinned Gaussian χ 2 analysis that takes into account the expected background event rates as outlined in Sec.IV (see Table I).As the backgrounds stem mainly from neutrino neutral-current interactions in the detector, we include one nuisance parameter (ξ) in order to account for systematic uncertainties affecting their overall normalization.This is done with the pull-method approach [128] taking a prior uncertainty on the background, σ bg = 20%, to account for the large uncertainties coming from the corresponding cross section (see, e.g., Refs.[127,129]).Thus, for a given decay channel ch with an expected non-zero background rate N bg,ch , we define our χ 2 simply as where T ch is the total expected event rate including both signal and background events, while O ch corresponds to the assumed observed events, and {Θ} stands for the parameters of the model.We take the observed events as expectation in the absence of a BSM signal, O ch = N bg,ch , and the associated statistical error as σ = √ O ch .The predicted event rates read: Note that the χ 2 definition above, in terms of total event rates, will typically lead to conservative results: an improvement in sensitivity may be obtained for a binned analysis that takes into account the different distributions of the signal versus the background in the kinematic variables of interest, which we leave for future work.
Our sensitivity regions are obtained taking the corresponding χ 2 cut at a given confidence level (C.L.), for 1 degree of freedom (d.o.f.).Thus, they can be interpreted as the upper limit that DUNE would be able to set on a given parameter, for an ALP with mass m a , in the absence of a BSM physics signal.Finally, for the channels without SM background, we follow the Feldman-Cousins prescription [130] and require N dec,ch > 2.44 for limits at 90% C.L.

A. Model-independent sensitivity limits
Using the χ 2 analysis just outlined, we first perform a model-independent sensitivity analysis.If we assume that the production branching ratio and the lifetime of the ALPs are independent, the number of decays approximately depends as in Eq. ( 5).This allows to derive a sensitivity limit on the product of the production and decay branching ratios as a function of cτ a /m a , shown in Fig. 6.Mild differences are obtained for different masses, however, induced by the dependence of the detector acceptance on m a (which affects the boost of the particles to the lab frame).Here we follow the same approach as in Ref. [25] and provide our limits as bands, where the width indicates the variation in the obtained limit when the mass of the LLP is varied between 10 MeV (upper edge of each band) and up to the production threshold in each case (lower edge of the band).In Fig. 6 we show two sets of bands, depending on the parent meson: kaons (blue) and D mesons (red).Moreover, for each parent meson, we show results for a → γγ using the ND-GAr, computed taking the corresponding background event rates from Table I, as well as the limiting sensitivity in the background-free case (which would only be applicable for decays into µ + µ − ).Note, however, that the upper edge of each band corresponds to m a = 10 MeV, for which the decays into µ + µ − or multi-pion final states are not kinematically accessible.The best sensitivity for a → µ + µ − searches is indicated by the dashed lines in each case, corresponding to m a ∼ 2m µ .
In order to compare to current limits in the literature it is convenient to pick a specific mass.This exercise is done in Fig. 7, where the show the DUNE sensitivity limits compared to previous constraints, for three representative masses: m a = 50 MeV (left) and 300 MeV (center), where the strongest limits at DUNE would be obtained for LLP produced in K decays; and m a = 1.2 GeV (right), which could be probed at DUNE only if the LLP is produced from D decays.
As shown in the figure, the sensitivity limits at DUNE depend heavily on the final state the LLP decays into, due to the different backgrounds expected.In particular, for ALP decays into γγ the analysis would be affected by large backgrounds, leading to a limit that is considerably worse than in the background-free scenario.This should be taken into consideration when comparing 10 −4 10 −3 10 −2 10 −1 10 0 10 1 10 2 10 3 10 4 cτ a (m) The different lines for DUNE (solid, dashed, dotted) are obtained for different background assumptions, depending on the final state being considered (a → γγ, a → ℓ + ℓ − , a → π 0 π + π − ), see Sec.IV.In the case of the right panel, the limit for a → ee would coincide with the one shown for a → 3π, as backgrounds would be similar for both final states, while the limit for a → µµ would be slightly better since it would be background free.Present constraints are indicated by the shaded areas and correspond to 90% C.L., with the exception of CHARM and LHCb (shown at 95% C.L.).
our results with similar studies in the literature obtained neglecting the effect of backgrounds.However, we see that DUNE is expected to improve over present constraints for masses below the kaon mass even in the less favorable case where the LLP decays into photons.In the case of decays into lepton pairs, DUNE is expected to reach values of BR(M → a) × BR(a → ℓ + ℓ − ) ∼ 10 −17 for optimal values of the LLP lifetime 4 , many orders of magnitude below current constraints.In case of a heavier LLP (right panel), the expected limits suffer from a significant reduction in the production of the parent meson, and we find that the limit for an ALP that decays purely into photons is not as strong as the one obtained from CHARM data (salmon shaded regions, with dotted edges).However, we see that DUNE will be able to considerably improve over current limits if the ALP decays predominantly into lepton pairs or into 3π, which would be background-free at the ND-GAr.

B. Sensitivity limits for specific scenarios
The computation of model-independent sensitivity limits (Figs. 6 and 7) is useful since allows our results to be easily recasted to other scenarios.However, one should bear in mind that in particular models, correlations may arise between the production and decay of the LLP if they depend on the same set of model parameters.This changes the relative importance of different sets of constraints, as these may be optimal for different values of the lifetime of the LLP.Moreover, while in Sec.V A the constraints shown are obtained for searches for invisible or visible LLP decays, in specific scenarios additional bounds arise (e.g. from SN1987A or meson mixing, see Sec.III).Thus, in the rest of this section we evaluate the sensitivity of DUNE for the benchmark scenarios considered in Sec.II, as illustrative examples.

Gluon dominance
As outlined in Sec.II, the main production mechanism in this scenario is through ALP mixing with neutral pseudoscalar mesons (for ALP masses below ∼ 1 GeV) and gluon fusion (for higher masses).Once produced, if the ALPs are sufficiently long-lived they will reach the DUNE detectors before decaying into either a pair of photons or multi-pion final states, see Fig. 2.
Here we revisit the results previously computed in Refs.[33,34] for this scenario, in light of our background estimates in Sec.IV and the refined computation of the relevant decay widths in Ref. [51].Our results are shown in Fig. 8, where we show separately the expected sensitivity regions for a search for a → γγ (red lines) and a → 3π (blue lines).Moreover, due to the different background rejection capabilities, we show separately the results for ND-LAr (dashed lines) and ND-GAr (solid lines).As can be seen from the figure, DUNE is expected to improve over current limits (shown by the shaded gray regions) only for ALP masses large enough to have a significant branching ratio to multi-pion final states, whereas for lighter ALP masses the sensitivity is affected by the large backgrounds expected for the di-photon channel.Also, note that even though the expected background for a → γγ is higher for ND-GAr than for ND-LAr, its much higher signal efficiency leads to a better performance thanks to the higher signal-to-background ratio.The vertical dotted lines indicate the masses of the neutral pseudoscalar mesons, where the large mixing with the ALP leads to sharp features in our regions (see also Fig. 2).

ALPs coupled through EW operators
As outlined in Sec.II, the main production mechanism at DUNE for this scenario is through kaon decays, K → πa.Once produced, the ALP may decay into three different final states: e + e − , µ + µ − and γγ, while hadronic modes are not allowed in this mass range (ALP decays a → ππ and a → π 0 γ are forbidden by CP and C, respectively, and the decay into three pions is kinematically not allowed in this mass window).
Our resulting sensitivity to ALPs coupled through EW operators is shown in Fig. 9. Let us first discuss the results shown in the left panel, obtained for an ALP coupled predominantly through the O ϕ operator.As outlined in Sec.II, in this case the ALP tends to decay preferentially to leptons.Once the muon decay channel is open, the ALP lifetime is considerably reduced and the sensitivity to a → e + e − is diminished in the region of large couplings.In the region of small couplings, on the other hand, the dependence on the total lifetime of the ALP approximately cancels and the result depends exclusively on the decay width for a given channel, a → ee, µµ.Overall, we see that DUNE has an excellent opportunity to considerably improve over current constraints for an ALP decaying into µ + µ − , by more than one order of magnitude, for both ND-LAr and ND-GAr.In the case of a → e + e − , an improvement is only expected for ND-GAr, while the results for ND-LAr are severely affected by the reduced signal efficiency, see Tab.I. Finally, the results for a → γγ are not shown in this case since the reduced branching ratio, combined with the much larger backgrounds expected, yields sensitivities that are not competitive in light of present bounds.
The right panel in Fig. 9 shows the results for an ALP coupled predominantly through the O W operator.While in this case the ALP tends to decay predominantly into photons, this channel is affected by much larger backgrounds.Therefore, although the branching ratio into leptons is suppressed for this operator (see Eq. ( 15)), the final sensitivities obtained for a → ℓℓ are similar and even surpass those obtained for a → γγ for large masses, when the decay channel a → µµ is opened.In particular, in the high-mass region we see that DUNE is expected to improve by almost an order of magnitude over current limits by the E137 experiment.

Charming ALPs
In this section, we present the sensitivity projections for the charming ALPs model described in Section II C, using the DUNE ND with a total exposure of N PoT = 1.1 × 10 22 .The sensitivity projections are shown in Fig. 4 separately for the final states with dominant branching ratios (a → γγ, a → π + π − π 0 ).The same way as in the previous scenarios considered, our sensitivity contours also exhibit here the typical shape of a visible decay search.For large values of the couplings, the ALPs become very short-lived and decay before reaching the detector, leading to a loss in sensitivity induced by the exponential term in the decay probability.On the other hand, small couplings are suppressed by both the production rate in Eq. ( 17) and the fact that the ALPs become too long-lived.
For decays into photon pairs, we include both production from kaon decays and from D decays.Conversely, hadronic decay channels are only available for m a ≳ 3m π and therefore only D → πa is considered in this case.In the case of a → γγ it is worth mentioning that in the region below for masses below m a < m K − m π the bound for this model is still dominated by the contribution from K → πa, since the suppressed production branching ratio for this scenario (see Eq. ( 19)) is ).Solid (dashed) lines correspond to ND-GAr (ND-LAr).Shaded gray regions indicate current bounds, see Sec.III.To facilitate the comparison with previous literature we provide on the right axes the corresponding limits to the so-called photon-dominance and fermion-dominance scenarios [46,47] (see Ref. [13] for the mapping to our Wilson coefficients).
partly compensated by the large kaon flux available.Again in this case we see that the resulting sensitivities for this decay channel are not competitive with current constraints, with similar results for the two near detectors.
As highlighted in Sec.II C (see Fig. 4) for low masses the branching ratio is dominated by the decay a → γγ whereas hadronic decay channels (driven mainly by a → π + π − π 0 ) rapidly take over for masses m a ≃ 3m π .Since the π 0 decays promptly to two photons, the final signal topology for those channels will be a → π + π − 2γ.Other relevant decay modes (not included here for simplicity) would be those involving η or η ′ mesons (e.g. a → π + π − η), which also decay promptly into two photons and would lead to a similar signature.From Fig. 4 we see that a general improvement is expected over present bounds, specially for ALP masses m a > 800 MeV.Note that for ALP masses in the vicinity of m η and m η ′ the decay width increases very rapidly, which translates into a sudden loss of sensitivity.

VI. SUMMARY AND CONCLUSIONS
The DUNE experiment will be exposed to the LBNF high-intensity neutrino beam.Besides the production of pions, kaons and other light mesons (such as η and η ′ ), the high proton energy at LBNF will gnerate of a significant flux of heavier mesons, such as D and D s .This provides a unique opportunity to study the production of exotic particles feebly coupled to the SM, with masses between O(10) MeV and O(2) GeV.If unstable, these could decay into SM particles and be searched for in the DUNE near detectors, which include both a liquid argon TPC (ND-LAr) and a gaseous argon TPC (ND-GAr).
The potential of DUNE to search for long-lived particles (LLPs) has been studied previously in  10.Sensitivity projections for our Benchmark Model III, as described in Sec.II C. Sensitivities are shown separately for a search into di-photon final states and for searches into final states with three pions, which are the dominant decay modes for this scenario.Solid (dashed) lines correspond to ND-GAr (ND-LAr).Shaded gray regions are disfavored by current constraints (see also Fig. 11).Vertical lines indicate the masses of the neutral pseudoscalar mesons, which explain the sharp features in the sensitivity regions.
the literature.However, most of those studies have neglected the impact of backgrounds, arguing that they could be reduced to a negligible level through appropriate selection cuts.In this sense it is worth noting that most studies in the literature consider the ND-GAr detector, which offers a superb environment for these searches thanks to the lower backgrounds expected from neutrino beam interactions.However, ND-GAr is expected to start taking data only during Phase II of the experiment.In light of this, it becomes pressing to assess the capabilities of ND-LAr for this kind of searches, which are a keystone for the BSM program at DUNE.
In this work, we have first computed the background rates according to the expected angular and energy resolution for both ND-GAr and ND-LAr detectors, for decays into pairs of photons, charged leptons or into three pions, with no missing energy.Our results (summarized in Table I) show that these may be reduced to a negligible level for channels involving pairs of muons in the final state; however, the background is significant for LLPs decaying into two photons, which dramatically reduces the expected final sensitivity.Our results also indicate a reduced ND-LAr signal efficiency for decays into e + e − or into three pions.In general, we note that our background study would be applicable to a wide range of BSM models with LLPs, including not only light pseudoscalars (which is the focus of this work), but also light scalars or vector bosons, as considered, for example, in the phenomenological studies in Refs.[26,32,35,36,38,40].
We have then derived sensitivity regions for generic LLP scenarios using the decay channels outlined above, within a model-independent approach.Our results (shown in Figs. 6 and 7, for the ND-GAr) are provided in terms of the production branching ratio and lifetime of the LLP, and may be easily recasted to specific scenarios, including both pseudoscalar and scalar particles.Overall, we find that DUNE has the potential to significantly improve over current constraints for LLP decays into muon pairs, where backgrounds may be reduced to a negligible level while keeping relatively high signal efficiencies.Conversely, we find that searches for decays into two-photons are significantly harder, due to the large backgrounds expected from neutrino interactions, contrary to expectation.
In order to put our sensitivities in context and to compare the potential of DUNE to present bounds, we then study three benchmark models involving axion-like particles (ALPs) coupling to the SM through effective operators: (i) the so-called gluon dominance scenario; (ii) a scenario where the ALP is coupled to electroweak operators; and (iii) the so-called charming ALP scenario, where the ALP is coupled to right-handed up-quark currents with a non-trivial flavor structure.Our choice of scenarios intends to demonstrate the potential of DUNE to constrain generic ALP models inducing significant couplings not only to photons (as commonly studied in the literature) but to charged leptons and mesons as well, over a wide range of masses.Overall we find that both ND-LAr and ND-GAr will be able to improve significantly over current limits in certain regions of parameter space, for decay channels involving muon pairs or multiple pions.This is particularly so in the case of ALPs produced from D decays, where the parameter space is largely unconstrained by laboratory experiments due to the intrinsic difficulties associated to the production of D mesons.However, we note that the sensitivities in the region where the ALP decays preferably into photon pairs would be affected by the large backgrounds expected, and improving over current limits will be challenging for DUNE.Finally, we stress that if the ALPs decay preferraly into e + e − pairs, ND-GAr will be needed in order to improve over current constraints, while the capabilities of ND-LAr are somewhat limited by the assumed signal efficiency.
In summary we have shown that, thanks to the high-intensity, high-energy proton beam available at LBNF, and in combination with the excellent capabilities of the argon gas TPC near detector, DUNE has the potential to improve significantly over current constraints for a wide landscape of BSM models with light unstable long-lived particles.We stress that similar searches using the liquid Argon TPC near detector may be possible.However, the higher backgrounds and lower efficiencies expected translate into significant reductions in sensitivity in most cases.In this regard, the availability of a gas TPC at DUNE will be a key asset in order to ensure the leadership of DUNE on LLP searches at neutrino facilities.
We emphasize that this work uses only publicly available information from the DUNE collaboration, and that the conclusions are our own and do not represent necessarily the official views of the DUNE Collaboration.
Note added: After this paper was finished, the authors of Ref. [131] pointed out the importance of effects from ALP mixing with neutral pseudoscalar mesons due to the coupling between ALPs and quarks.Although the focus of Ref. [131] is the fermion universal case, such effects would also be relevant for the benchmark models considered here.Mixing effects may lead to an increase in the expected sensitivity for DUNE in certain regions of parameter space, as new production mechanisms become available; however, previous constraints would also have to be reevaluated accordingly.We leave a detailed study of mixing effects between ALPs and neutral pseudoscalars (including a proper reevaluation of past constraints) for future work.

ACKNOWLEDGMENTS
We warmly thank Pilar Hernández for her involvement in the early stages of this work, and Carlos Pena for illuminating discussions.We also thank Adrián Carmona, Felix Kahlhoefer, Luca Merlo, Laura Molina-Bueno, Christiane Scherb and Edoardo Vitagliano for useful discussions and comments.This work has received partial support from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant procedure, we first check that the resulting χ 2 gives a good fit to the extracted data within the error bars provided in Ref. [95], considering the number of degrees of freedom in the fit (n dof ).At the best-fit, we obtain χ 2 min = 29 for n dof = 28 .
This leads to a p-value of 41%, indicating good compatibility with the data.Indeed, our best-fit curves match very well those in Ref. [95], as they should.Next, we have checked that we approximately reproduce their fit to the D + → τ + ν.Specifically, we obtain a best-fit for the number of events which corresponds to an overall precision of about 30%.This reproduces reasonably well the result of the collaboration (137 ± 27 events) albeit with a larger error bar in our case.Thus, our limits may be taken as conservative, keeping in mind that a dedicated analysis (including all data obtained in the full energy range, and a more sophisticated implementation of systematic uncertainties) would probably lead to a better result.
Given the good agreement between our result and that obtained in Ref. [95] for the SM case, we then proceed to add a contribution from D → πa events as a function of m 2 a (which can be can be directly identified with MM 2 ).To this end, we approximate the energy resolution of the detector in MM 2 by fitting the distribution of D → µν to a Gaussian.We believe this should be approximately correct, given that the signal in this case should be a delta function centered at zero.We obtain the best-fit for a Gaussian with mean σ ≃ 0.026 GeV 2 and a small bias in its central value, µ = 0.016 GeV 2 (in what follows, we assume both parameters are independent of MM 2 in the energy range under consideration).The differential distribution of the expected number of events from D → πa in each sample s is computed for each value of the ALP mass, as: where N D stands for the total number of D mesons produced, while ϵ s is the efficiency associated to each sample, ϵ π is the pion detection efficiency (which we set to ϵ π = 90% according to Ref. [95]), and P exit is the probability for an ALP to exit the BESIII detector before decaying.
The probability P exit is estimated assuming that the parent D mesons are produced in pairs, in collisions with a center-of-mass energy √ s = 3.77 GeV, and taking the detector radius ∼ 3 m, following Ref.[132].Finally, the number of D mesons is estimated from the observed number of events and the branching ratio reported by the collaboration in Ref. [95], taking N ev = N D BR(D → τ ν)BR(τ → π)ϵ π , BR(D → τ ν) exp = 1.2 × 10 −3 and N ev = 137.The contribution of the number of D → πa events to each bin in MM 2 is obtained integrating Eq. (A4) within the limits of each bin.These are subsequently added to the total number of predicted signal events, as an additional contribution to the χ 2 in Eq. (A1).A limit can then be obtained on the production branching ratio BR(D → πa), after marginalization over the same set of nuisance parameters as in the SM case.The resulting limit is shown in Fig. 11 at 90% CL, for a χ 2 with n dof = 28.18)) taking fa/Λ = ϵ = mc/mt.The different regions indicate current bounds from NA62 [89,133], E787 & E949 [91], CHARM [134], and BESIII [95], all of which have been computed here, see Sec.III for details.On the other hand, bounds from SN1987a and D − D mixing are directly taken from Ref. [62].

FIG. 1 .
FIG.1.Detector acceptance as a function of the energy of the ALPs.On the left panel, the ALPs have a mass of 300 MeV and originate from the decay K + → π + a.On the right panel, the ALPs have a mass of 1 GeV and are produced via D + → π + a.

c ϕ = 1 FIG. 3 .
FIG. 3. Main decay widths for the ALP considered in our benchmark model II (BM-II), as a function of its mass ma.Results for an ALP coupled predominantly through OW are shown in the left panel, while the right panel shows similar results for O ϕ .

1
FIG. 5. Distributions for signal (blue line) and background (grey solid histogram) of the kinematic quantities used in our event selection.The vertical dashed red line indicates the selection cut value used.All histograms are area normalized.Top row: transverse momentum (left panel) and angle between tracks (right) for µ + µ − candidate ALP-decay events in the DUNE ND-GAr.Central row, left: angle with respect to the beam direction of e + e − events.Central row, right: angle between the two photon showers for γγ candidate events in ND-LAr.Bottom row: transverse momentum (left panel) and angle between the charged-pion tracks of π + π − π 0 candidate events.

5 DUNE ND-GAr 1 . 1 •
FIG.6.Expected sensitivity to LLPs in the model-independent scenario, assuming their branching ratios and lifetime are completely uncorrelated.In the absence of a new physics signal, DUNE is expected to disfavor the region above each line at 90% confidence level (C.L.).The width of the bands indicates the variation in our results when the mass of the LLP is varied between 10 MeV (upper edge of each band) up to the production threshold in each case; see text for details.The dashed lines correspond to ma ∼ 2mµ, and therefore indicate the best sensitivity for searches of a → µ + µ − .Blue (red) bands show the results for LLP produced from kaon (D-meson) two-body decays.For each parent meson, the different bands are obtained assuming different background rates; see main text for details.

4 cτFIG. 7 .
FIG.7.Expected sensitivity to LLPs as a function of their lifetime.The different panels show the results for ma = 50 MeV (left), 300 MeV (center) and 1.2 GeV (right), as indicated.Colored regions are disfavored by present constraints from searches for invisible and visible LLP decays, see Sec.III.The different lines for DUNE (solid, dashed, dotted) are obtained for different background assumptions, depending on the final state being considered (a → γγ, a → ℓ + ℓ − , a → π 0 π + π − ), see Sec.IV.In the case of the right panel, the limit for a → ee would coincide with the one shown for a → 3π, as backgrounds would be similar for both final states, while the limit for a → µµ would be slightly better since it would be background free.Present constraints are indicated by the shaded areas and correspond to 90% C.L., with the exception of CHARM and LHCb (shown at 95% C.L.).

DUNE 1 . 1 •
FIG.8.DUNE sensitivity projections for the gluon dominance scenario, Eq. (6).Shaded regions are disfavored by present experiments.Red (blue) lines show the expected DUNE sensitivity (at 90% CL) for a search into a → γγ and a → 3π final states.Results are shown separately for the ND-LAr (dashed lines) and ND-GAr (solid lines).

1 a 1 •
FIG.9.DUNE sensitivity projections (one Wilson coefficient switched on at a time) for c ϕ (left panel) and cW (right panel) as a function of ma, assuming fa = 1 TeV.Our sensitivity lines are shown at 90% CL, individually for each final state topology as indicated by the labels (a → γγ, a → e + e − and a → µ + µ − ).Solid (dashed) lines correspond to ND-GAr (ND-LAr).Shaded gray regions indicate current bounds, see Sec.III.To facilitate the comparison with previous literature we provide on the right axes the corresponding limits to the so-called photon-dominance and fermion-dominance scenarios[46,47] (see Ref.[13] for the mapping to our Wilson coefficients).

FIG. 11 .
FIG. 11.Summary of present bounds on applicable to the Charming ALP scenario, as a function of ma, for the FN-inspired model (see Eq. (18)) taking fa/Λ = ϵ = mc/mt.The different regions indicate current bounds from NA62[89,133], E787 & E949[91], CHARM[134], and BESIII[95], all of which have been computed here, see Sec.III for details.On the other hand, bounds from SN1987a and D − D mixing are directly taken from Ref.[62].