Heavy Axion Opportunities at the DUNE Near Detector

While the QCD axion is often considered to be necessarily light ($\lesssim$ eV), recent work has opened a viable and interesting parameter space for heavy axions, which solve both the Strong CP and the axion Quality Problems. These well-motivated heavy axions, as well as the generic axion-like-particles, call for explorations in the GeV mass realm at collider and beam dump environments. The primary upcoming neutrino experiment, Deep Underground Neutrino Experiment (DUNE), is simultaneously also a powerful beam dump experiment, enabled by its multipurpose Near Detector (ND) complex. In this study, we show with detailed analyses that the DUNE ND has a unique sensitivity to heavy axions for masses between $20$ MeV and $2$ GeV, complementary to other future experiments.

The QCD axion provides an elegant dynamical solution to the Strong CP Problem of the Standard Model (SM) [1][2][3][4]. In its minimal realization, the QCD axion signatures are dominated by its couplings to the CP-odd field strength operators of SM gauge fields and the SM matter fields. The well-known relation (see e.g. [5]) for the QCD axion, relating its mass m a to its decay constant f a -m a = 5.7 10 9 GeV fa meV -along with a variety of experimental/observational constraints implying f a > 10 9 GeV [6], have driven most QCD axion searches to focus on light, sub-eV masses. However, a series of recent model building efforts [7][8][9][10][11][12][13][14][15][16][17] including earlier work [18][19][20][21][22][23] motivate heavier variants of the QCD axion and, within a class of such models relying on a Z 2 −symmetric mirror SM sector, a testable parameter space has been identified where the axion mass can be around or even larger than the GeV scale [15]. In such scenarios, the heavy axion can be obtained by introducing a new strongly-coupled mirror SU (3) sector that also generates a larger axion potential aligned with the QCD-generated potential [7-9, 15, 22, 23]. As a result the axion continues to solve the Strong CP problem, while being heavy enough (for a given f a ) to be more robust against unwanted UV contributions which would otherwise have given rise to the so-called Quality Problem [24][25][26][27]. Therefore, this predicts a new, less-explored heavy axion solving both the Strong CP and the Quality Problem that deems further exploration. Moreover, a pseudoscalar field is a generic constituent of many beyond-the-Standard-Model (BSM) scenarios [28], as well as String Theoretic constructions [29,30]. Hence, the search for GeV-scale pseudoscalar fields, parametrized under a generic effective field theory, dubbed as "Axion-like Particles" (ALP), is a vital component of the BSM program.
The next-generation Deep Underground Neutrino Experiment (DUNE) [58], via its intense, high energy Long Baseline Neutrino Facility (LBNF) proton beam and high-precision near and far detectors, will facilitate the most accurate measurement of various neutrino properties. At the same time, its Near Detector (ND) site [58], composed of a complex detector facility with a full-fledged particle detection and identification system, can be creatively viewed as a beam dump facility for new physics searches for long-lived particle signatures [59][60][61][62]. As we demonstrate in this work, the heavy axions, well-motivated by the Strong CP and Quality Problems, and generic BSM considerations, would provide a compelling target for DUNE ND. Our proposed search at DUNE ND can cover unique parts of the axion parameter space.
This work is organized as follows. In Sec. II, we briefly review the mechanism solving both the Strong CP and Quality Problems and how that motivates a heavy axion. After discussing some of its properties and the EFT parametrization, in Sec. III, we detail our simulation procedure for axion production in the proton beam dump environment with the DUNE ND. That allows us to investigate in Sec. IV both the signatures of heavy axion decay in DUNE ND and possible background contributions. All of these considerations then enable us to derive the DUNE ND sensitivity to heavy axions in Sec. V. When projecting sensitivities, we consider two benchmark models, both orthogonal to the "photon-dominant" axion(-like-particles) frequently considered in the literature. We assume that either (a) the heavy axion coupling is "gluon-dominated", i.e., the coupling between the axion and the QCD field strength tensor is the largest among its other couplings, or (b) the heavy axion coupling is "co-dominant" and the axion couples equally to the different SM field strength tensors of SU (3), SU (2) and U (1) Y . We then conclude in Sec. VI.
Distance from Target (not to scale) < l a t e x i t s h a 1 _ b a s e 6 4 = " k 3 N p B R 1 F X R s u I 6 q d E y 2 r y N r W 4 n k = " > A A A B 6 X i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S I I Q k l E 0 W P R i 8 c q 9 g P a U D b b S b t 0 s w m 7 G 6 G G / g M v H h T x 6 j / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n a X l l d W 1 9 c J G c X N r e 2 e 3 t L f f 0 H G q G N Z Z L G L V C q h G w S X W D T c C W 4 l C G g U C m 8 H w Z u I 3 H 1 F p H s s H M 0 r Q j 2 h f 8 p A z a q x 0 f / r U L Z X d i j s F W S R e T s q Q o 9 Y t f X V 6 M U s j l I Y J q n X b c x P j Z 1 Q Z z g S O i 5 1 U Y 0 L Z k P a x b a m k E W o / m 1 4 6 J s d W 6 Z E w V r a k I V P 1 9 0 R G I 6 1 H U W A 7 I 2 o G e t 6 b i P 9 5 < l a t e x i t s h a 1 _ b a s e 6 4 = " H m m E X P 4 V x 2 q s v H w 6 r 5 t 8 8 8 7 E h m E h m s i I y x B I T b Q o r m x L c + S 8 v k l a t 6 p 5 X n b u z S v 2 q q K M E R 3 A M p + D C B d T h F h r Q B A K P 8 A y v 8 G Y 9 W S / W u / U x G 1 2 y i p 0 D + A P r 8 w c F 7 p R R < / l a t e x i t > p 120 GeV < l a t e x i t s h a 1 _ b a s e 6 4 = " K k 1 W n u N T Y i u H A y P 9 6 0 1 t e x s i h / Y = " > A A A B 7 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K o k o e i x 6 8 V j B t I U 2 l s 1 2 0 y 7 d b M L u R C i h v 8 G L B 0 W 8 + o O 8 + W / c t j l o 6 4 O B x 3 s z z M w L U y k M u u 6 3 s 7 K 6 t r 6 x W d o q b + / s 7 u 1 X D g 6 b J s k 0 4 z 5 L Z K L b I T V c C s V 9 F C h 5 O 9 W c x q H k r X B 0 O / V b T 1 w b k a g H H K c 8 i O l A i U g w i l b y u 6 l 4 d H u V q l t z Z y D L x C t I F Q o 0 e p W v b j 9 h W c w V M k m N 6 X h u i k F O N Q o m + a T c z Q x P K R v R A e 9 Y q m j M T Z D P j p 2 Q U 6 v 0 S Z R o W w r J T P 0 9 k d P Y m H E c 2 s 6 Y 4 t A s e l P x P 6 + T Y X Q d 5 E K l G X L F 5 o u i T B J M y P R z 0 h e a M 5 R j S y j T w t 5 K 2 J B q y t D m U 7 Y h e I s v L 5 P m e c 2 7 r L n 3 F 9 X 6 T R F H C Y 7 h B M 7 A g y u o w x 0 0 w A c G A p 7 h F d 4 c 5 b w 4 7 8 7 H v H X F K W a O 4 A + c z x 9 0 8 I 5 z < / l a t e x i t > ⇡ 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " s i k l I g q f V g F G D O M m K i X m T p / 5 f u k = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 W P R i 8 c W 7 A e 0 o W y 2 k 3 b t Z h N 2 N 0 I J / Q V e P C j i 1 Z / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 n G q G D Z Z L G L V C a h G w S U 2 D T c C O 4 l C G g U C 2 8 H 4 b u a 3 n 1 B p H s s H M 0 n Q j + h Q 8 p A z a q z U o P 1 y x a 2 6 c 5 B V 4 u W k A j n q / f J X b x C z N E J p m K B a d z 0 3 M X 5 G l e F M 4 L T U S z U m l I 3 p E L u W S h q h 9 r P 5 o V N y Z p U B C W N l S x o y V 3 9 P Z D T S e h I F t j O i Z q S X v Z n 4 n 9 d N T X j j Z 1 w m q U H J F o v C V B A T k 9 n X Z M A V M i M m l l C m u L 2 V s B F V l B m b T c m G 4 C 2 / v E p a F 1 X v q u o 2 L i u 1 2 z y O I p z A K Z y D B 9 d Q g 3 u o Q x M Y I D z D K 7 w 5 j 8 6 L 8 + 5 8 L F o L T j 5 z D H / g f P 4 A x K + M 6 Q = = < / l a t e x i t > a < l a t e x i t s h a 1 _ b a s e 6 4 = " s i k l I g q f V g F G D O M m K i X m T p / 5 f u k = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 W P R i 8 c W 7 A e 0 o W y 2 k 3 b t Z h N 2 N 0 I J / Q V e P C j i 1 Z / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 n G q G D Z Z L G L V C a h G w S U 2 D T c C O 4 l C G g U C 2 8 H 4 b u a 3 n 1 B p H s s H M 0 n Q j + h Q 8 p A z a q z U o P 1 y x a 2 6 c 5 B V 4 u W k A j n q / f J X b x C z N E J p m K B a d z 0 3 M X 5 G l e F M 4 L T U S z U m l I 3 p E L u W S h q h 9 r P 5 o V N y Z p U B C W N l S x o y V 3 9 P Z D T S e h I F t j O i Z q S X v Z n 4 n 9 d N T X j j Z 1 w m q U H J F o v C V B A T k 9 n X Z M A V M i M m l l C m u L 2 V s B F V l B m b T c m G 4 C 2 / v E p a F 1 X v q u o 2 L i u 1 2 z y O I p z A K Z y D B 9 d Q g 3 u o Q x M Y I D z D K 7 w 5 j 8 6 L 8 + 5 8 L F o L T j 5 z D H / g f P 4 A x K + M 6 Q = = < / l a t e x i t > a < l a t e x i t s h a 1 _ b a s e 6 4 = " s i k l I g q f V g F G D O M m K i X m T p / 5 f u k = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 W P R i 8 c W 7 A e 0 o W y 2 k 3 b t Z h N 2 N 0 I J / Q V e P C j i 1 Z / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 n G q G D Z Z L G L V C a h G w S U 2 D T c C O 4 l C G g U C 2 8 H 4 b u a 3 n 1 B p H s s H M 0 n Q j + h Q 8 p A z a q z U o P 1 y x a 2 6 c 5 B V 4 u W k A j n q / f J X b x C z N E J p m K B a d z 0 3 M X 5 G l e F M 4 L T U S z U m l I 3 p E L u W S h q h 9 r P 5 o V N y Z p U B C W N l S x o y V 3 9 P Z D T S e h I F t j O i Z q S X v Z n 4 n 9 d N T X j j Z 1 w m q U H J F o v C V B A T k 9 n X Z M A V M i M m l l C m u L 2 V s B F V l B m b T c m G 4 C 2 / v E p a F 1 X v q u o 2 L i u 1 2 z y O I p z A K Z y D B 9 d Q g 3 u o Q x M Y I D z D K 7 w 5 j 8 6 L 8 + 5 8 L F o L T j 5 z D H / g f P 4 A x K + M 6 Q = = < / l a t e x i t > a ⌘/⌘ 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " J A E 3 R y A Q W H B v A O + 2 D K n j 8 s h z 5 F c = " > A A A B + X i c b V D L S s N A F J 3 U V 6 2 v q E s 3 g 0 V w V Z M q 6 L L o x m U F + 4 A m l s n 0 p h 0 6 k 4 S Z S a G E / o k b F 4 q 4 9 U / c + T d O 2 i y 0 9 c C 9 H M 6 5 l 7 l z g o Q z p R 3 n 2 y q t r W 9 s b p W 3 K z u 7 e / s H 9 u F R W 8 W p p N C i M Y 9 l N y A K O I u g p Z n m 0 E 0 k E B F w 6 A T j u 9 z v T E A q F k e P e p q A L 8 g w Y i G j R B u p b 9 s e a H K R t y c v k U x A 3 6 4 6 N W c O v E r c g l R R g W b f / v I G M U 0 F R J p y o l T P d R L t Z 0 R q R j n M K l 6 q I C F 0 T I b Q M z Q i A p S f z S + f 4 T O j D H A Y S 1 O R x n P 1 9 0 Z G h F J T E Z h J Q f R I L X u 5 + J / X S 3 V 4 4 2 c s S l I N E V 0 8 F K Y c 6 x j n M e A B k 0 A 1 n x p C q G T m V k x H R B K q T V g V E 4 K 7 / O V V 0 q 7 X 3 M t a / e G q 2 r g t 4 i i j E 3 S K z p G L r l E D 3 a M m a i G K J u g Z v a I 3 K 7 N e r H f r Y z F a s o q d Y / Q H 1 u c P Z U W T g Q = = < / l a t e x i t > FIG. 1. Simplified schematic representation of the LBNF beam, DUNE target, and DUNE Near Detector complex for our study. Protons with 120 GeV energy strike a target that produces copious amounts of SM mesons and, potentially, axions a. Some fraction of these particles travel in the direction of the Near Detector complex (right) where the liquid argon time projection chamber (green) and multi-purpose detector (red indicating the gaseous argon time-projection chamber and blue indicating the electromagnetic calorimeter/magnet) are situated. The axions can decay in those detectors providing striking signals. The crosses on the π 0 and η/η lines indicate that these particles can directly mix with the axions a. Figure adapted and modified from Ref. [60].

II. PROPERTIES OF A HEAVY AXION
This section details the properties of a heavy axion considered in this work. First, in section II A, we introduce such a heavy axion that solves the Strong CP Problem. After discussing the axion Quality Problem, in Section II B, we demonstrate how such a heavy axion can solve both these problems simultaneously and thereby be theoretically well-motivated. After reviewing the theoretical constraints on such a heavy axion, in Section II C, we describe the axion properties under benchmark choices of the axion effective field theory (EFT) to study experimental prospects in the later sections.
Before proceeding, we wish to clarify some definitions and conventions that will be used throughout this work. The (dimensionful) coupling f a is the axion decay constant that enters as a part of axion coupling to the CP-odd QCD field strength tensor, responsible for solving the Strong CP problem. We will also occasionally use the quantity f G = 4π 2 f a , which is adopted in some phenomenological studies of heavy axions in the literature. Finally, since production rates and decay widths of a are typically inversely-proportional to f a and f G , we will often use the coupling g agg = f −1 G so that these rates/widths are proportional to positive powers of g agg , which also allows for more transparent comparisons between our projections and those in the literature.

A. The Strong CP and the Quality Problem
As defined, the QCD axion has a gluonic coupling (for a recent review see [63]) through which it solves the Strong CP Problem dynamically. This can be seen explicitly by considering the QCD-generated axion potential [3,5,64] As the axion dynamically relaxes to its minima at a = −θf a , it makes the effectiveθ parameter vanish -solving the Strong CP Problem and explaining the smallness of the (as-yet unobserved) neutron electric dipole moment [65]. The above QCD-axion potential gives rise to the well-known relation [5], However, a variety of terrestrial, astrophysical and cosmological constraints (see e.g. [6] and references therein) require f a > 10 9 GeV for the QCD axion and therefore, the relation in Eq. (3) precludes observing an otherwise phenomenologically interesting, accelerator-observable parameter space where m a ∼ 10 MeV -100 GeV. Hence, it is interesting to ask whether there exist models solving the Strong CP Problem which can occupy this mass regime. At the same time, from an ultraviolet (UV) perspective, various axion models often suffer from the so-called "Quality Problem" [24][25][26][27]. To see this, we recall that an axion can be realized as the Goldstone boson of a spontaneously broken Peccei-Quinn (PQ) U (1) PQ symmetry. In the far UV, quantum gravitational effects are expected to break all global symmetries [66,67]. Therefore, at low energies, the U (1) PQ can, at best, survive as some accidental symmetry. 1 To illustrate the severity of the Quality Problem, we consider a Planck-suppressed U (1) PQ -breaking operator, In the above, the axion a arises as the Goldstone mode of the (composite) field Φ ∼ f a e ia/fa at low energies, and the resulting minima-structure of the axion potential need not align with the QCD-generated potential. Therefore, unless the UV contribution in Eq. (4) is small compared to the QCD-generated potential in Eq. (2), the axion will relax to a minima dictated by Eq. (4) where generally a = −θf a and the axion solution to the Strong CP Problem will be spoiled. More concretely, given the constraint f a > 10 9 GeV for the minimal QCD axion, we see that unless we forbid all operators of the type in Eq. (4) up to N = 9, the axion-solution to the Strong CP Problem no longer works. The question of why the UV theory should respect U (1) PQ to such high quality is the Quality Problem.

B. Addressing the Quality Problem
In Ref. [15] a model addressing the Quality Problem was constructed, generalizing on previous work [7-9, 22, 23], in which the Strong CP Problem is solved through the presence of a Z 2 -symmetric mirror sector (containing the primed fields). The axion coupling is then given by, This Z 2 symmetry is softly broken by the only relevant operator in the SM and the mirror sector, the Higgs masses. Consequently, the mirror Higgs VEV H can "naturally" be much larger than H . In this case, the mirror quarks decouple at higher energies without impacting the RG running of the mirror QCD at lower energies. This results in the mirror confinement scale Λ QCD being much larger than Λ QCD . Therefore, the QCD-axion potential receives a parametrically larger contribution from the mirror sector with [15], At the same time, thanks to the Z 2 symmetry, this enhanced contribution still aligns with SM-QCD generated potential in Eq. (2). Thus the same axion solves the Strong CP Problem in both the sectors, and, since the axion potential is parametrically enhanced due to the presence of the mirror sector, it is less susceptible to the Quality Problem. We now briefly summarize the theoretical constraints on our model in Fig. 2, while refering the reader to Ref. [15] for more detailed explanations.
Theoretical constraints on the axion parameter space for the class of models considered in this work that solve both the Strong CP and the Quality Problems, adapted from Ref. [15]. The white region is the theoretically allowed/motivated region. See the text for explanations of different labels. The parameters f G and f a are related by f G = 4π 2 f a . Fig. 2 presents the theoretically-motivated region of parameter space for this heavy, high-quality axion model, as a function of the axion mass m a and the axion decay constant f a ≡ f G /(4π 2 ). Since our mechanism makes the axion only heavier, it can not populate the region labeled "QCD Axion" where it would be lighter than the QCD axion. In the region labeled "f a < Λ QCD ", the axion EFT breaks down because m a > f a in that region. In the region labelled "PQ Quality Problem," the axion suffers the Quality Problem discussed in Eq. (4) due to operators with N ≥ 6. Finally, in the region in the bottom right denoted "H Quality Problem", the mirror VEV H spoils the Strong CP solution via, In particular, to avoid the Quality Problem from Eq. (7), we require H < 10 14 GeV, so there is only a maximal amount by which the axion can be made heavier in this scenario. By inspection, Fig. 2 encourages us to focus on heavy axions in the keV-TeV mass range with f G between 10−10 9 GeV. A natural question that emerges is how much of the open theoretical parameter space can be covered by existing and upcoming experiments. To discuss this, we detail a phenomenological discussion of the heavy axion properties next.

C. Heavy Axion EFT, Mixing and Lifetime
A robust consequence of the above mentioned class of models is the defining GG coupling of the axion. For this purpose, we consider an effective Lagrangian, with α i = g 2 i /(4π) given in terms of SM gauge couplings, and α 1 = 5/3α Y , in terms of the hypercharge gauge coupling. To illustrate the significance of a non-zero c 3 , we will focus on two scenarios which are complementary to the case of photon or electroweak-dominance, c 2 , c 1 c 3 , frequently assumed in the literature due to its testability. In more detail, we will focus on the cases of, • Gluon dominance: c 3 = 1, c 1 , c 2 = 0; Both of the above cases are motivated from the generic Axion considerations, as well as UV considerations. In fact, these two cases match well respectively to the KSVZ [68,69] and DFSZ [20,70] scenario of the minimal axion theory.
These choices have an important effect on the phenomenology of such heavy axions, since for m a 1 GeV, the axions predominantly decay into hadronic final states, as opposed to diphoton final states on which a significant number of searches rely. Consequently, interesting parts of the axion parameter space open up, as we will see below. Simultaneously, a non-negligible c 3 gives rise to important axion production channels at the LHC and various proton beam dump experiments.
Below the scale of electroweak symmetry breaking, the EFT in Eq. (8) gives rise to an axion photon coupling, with [42,48,71] To obtain c γ for m a Λ QCD , we have assumed the η − η mixing angle sin θ ηη = −1/3 as in [48] while noting the significant uncertainty θ ηη −(10 • − 20 • ) (see Sec. 15 of Ref. [6] and references therein). Some results for more general mixing angles can be found in Refs. [54,72]. The factor of 1.92 in Eq. (10) can be obtained after including higher order corrections [5] on the leading order contribution due to quark masses 2 3 4m d +mu mu+m d ≈ 2. Importantly, we see that even in the absence of the a tree-level c 1 , c 2 , the anomaly and axion-pseudoscalar meson mixing introduce a signigicant axion-photon coupling below Λ QCD .
The phenomenology will be largely dictated by the lifetime of the axion. For m a < 3m π , the axion decays exclusively in diphoton final states with a width, whereas above that threshold the hadronic decay modes open up [48] and quickly become dominant except regions near resonant mixing. We show the resulting lifetime of the axion in Fig. 3 which is used in the following to determine the reach of the DUNE ND. The blue (red) line in Fig. 3 assumes the Gluon Dominance (Codominance) scenario, and both lines assume f G = 1 PeV. This lifetime is proportional to f 2 G . We note that the two scenarios are nearly identical for m a 1 GeV but below 1 GeV the above distinctions are quite important. For m a below 100 MeV, we can see that in the Codominance scenario, the axion has a larger lifetime as the corresponding c γ is smaller for the particular choices of c i = 1. Between 100 MeV and 1 GeV, two effects are important. One is the near resonance mixing with the SM mesons, which determines the dips in this lifetime plot. The other one is the cancellation between different meson mixings and the direct contributions to c γ in Eq. 10 for which we get peaks in the lifetime. We also note here, the mixing expansion in this equation at the meson pole regime should be regulated by the unitarity of the mixing matrix, which we neglected here in the equations but implemented effectively in the next section in our numerical computation.

III. SIMULATION DETAILS: AXION PRODUCTION AND DUNE NEAR DETECTOR
This section details the simulations we perform and how we determine the DUNE experimental sensitivity to heavy axions. In Section III A we explain the approach we employ to calculate the axion production, both from meson mixing and from gluon-gluon fusion. Section III B explains how we include the DUNE ND complex in these simulations, including both the liquid argon near detector and the gaseous argon multi-purpose detector.

A. Axion Production Details
Meson Mixing: To determine the axion production due to meson mixing, we simulate the Long Baseline Neutrino Facility (LBNF) beam as a 120 GeV proton beam colliding with a fixed target using Pythia8 with the "SoftQCD:all = on" option. For all of our simulations, we assume the total number of protons-on-target (POT) is N POT = 1.47 × 10 22 over the course of 10 years. 2 We find that approximately 2.89 π 0 , 0.33 η, and 0.03 η are produced per POT at this beam energy.
As alluded to in Eq. (10), axions mix with the SM pseudoscalar mesons through the GG coupling. Here we summarize the mixing angles [42,48,72], π =π phys + θ aπ a phys + · · · ≈ π phys + 1 6 η =η phys + θ aη a phys + · · · ≈ η phys + 1 η =η phys + θ aη a phys + · · · ≈ η phys + where f π ≈ 93 MeV. In these equations, the ellipses contain π − η and π − η mixing terms, which subdominantly contribute to ALP production considered below. Using Eqs. (12), (13) and (14), the number of axions produced from ALP-meson mixing is obtained as 3 , The axion flux with axion mass below 1 GeV at the DUNE ND then depends on Eq. (15), the ND cross-sectional area, and the acceptance fraction. The acceptance fraction is defined as the fraction of produced axions that are travelling in the direction of the DUNE ND upon production, folding in the production angular depenendence in the beam dump environment. We discuss the details of the DUNE ND in our simulations in Section III B, and we note here that the acceptance fraction for the different meson-mixing production mechanisms is O(10 −2 ). Fig. 4 displays the expected flux at the DUNE ND for axion production from meson mixing as a function of the mass m a . We show the separate contributions from π 0 , η, and η as different colors, which allows us to see the different axion-meson mixing dominating for different regions of m a . We dash the contributions for m a 1 GeV, where we expect gluon-gluon fusion to serve as a better description of axion production in this environment. This flux is shown for f G = 1 PeV and scales with f −2 G . Gluon-Gluon Fusion: Above O(GeV) masses, the direct production mode from gluon-gluon fusion could potentially dominate the contribution to the heavy axion flux at beam dump facilities. There, given the momentum exchange is above the GeV scale, the parton distribution function description is valid. The operator ∝ aGG determines the production rate.
We evaluate the production cross section convoluted with the leading order parton distribution function NNPDF [73,74] with our calculation available at this link, following: where µ F and µ R are the factorization and renormalization scale, and s is the center of mass energy, approximately (15 GeV) 2 . The factorization and renormalization scale is set at the heavy axion mass µ 2 R = µ 2 F = m 2 a . In Fig. 5(left), we show the inclusive production rate for a decay constant f G = 1 PeV. We used the 1-loop renormalization-group-equation running of the strong coupling constant α S (Q 2 ) embedded in NNPDF. In Fig. 5(right), we show the same rate in a linear scale, effectively zooming in around m a ≈ 1 GeV. We also show how the cross section varies when we consider factorization scales between m 2 a , 1/2m 2 a (dashed lines) and 2m 2 a (dotted). In both panels, we see the scale uncertainty is sizeable, especially below the GeV scale. In fact, this regime is where PDFs have large uncertainties and scheme dependence. We observe with our current NNPDF choice the cross section only starts to dominate when m a 0.9 GeV, and so our results are not subject to the large uncertainty in the low mass regime. Another interesting phenomena is the cross-over of different scale choices, this is a result of the PDF evolution and the α 2 S (Q 2 ) running from the production cross section. In Fig. 6, we show the normalized energy distribution from the gluon-gluon fusion production for a benchmark heavy axion with a mass of 1.2 GeV. The red curves show the distribution at production, which is universal for all axion decay constants. For m a √ s , where √ s 15 GeV for the LBNF-DUNE beam under consideration here, the distribution from gluon-gluon fusion are similar. In the same plot, we also show the differential distribution in orange, yellow and cyan for the axion accepted by the DUNE ND (accounting for decays of axions en route) for axion decay constant f G of 50, 100, and 200 PeV, respectively. The shift in distribution is mainly driven by the necessary boost for a heavy axion to arrive at the DUNE ND before decaying. More boost is needed for a smaller decay constant, and hence the distribution shifts towards higher energies. We discuss the detection considerations and details in Section IV.
A few other axion production modes could be significant. These include bremsstrahlung effects from the proton-proton collinear emission with resummation. The result will depend on how one treats the finite mass effect from the axion and its derivative coupling. Similarly, there can be collinear emissions from the quarks and gluons in the collision, involving model dependent axion-quark couplings. Through the axion coupling to quarks, possible flavor-changing decays from mesons will also contribute to axion production [75]. Last but not least, the proton-proton collision at beam dumps will create secondary collisions from the remnants of the first collision, enlarging the number of mesons produced and hence enriching the flux for heavy axions below O(GeV). All these effects could help improve the axion flux and, therefore, the DUNE ND sensitivities to heavy axions. We leave a detailed analysis of these different contributions with more model dependence for future studies. TABLE I. Signals of ALP decay a → γγ and a → hadrons in the liquid argon near detector (ArgonCube) and the gaseous argon detector (MPD). We also list the dominant source of backgrounds in the detectors for each of these searches, and some properties that distinguish between the signals and backgrounds.

Signature
Liquid We are interested in signatures of axion decay in the DUNE ND Complex. Specifically, we consider such signatures inside the liquid argon time-projection-chamber detector ArgonCube and the gaseous argon time-projection-chamber Multi-Purpose Detector (MPD). ArgonCube is situated at a distance of 574 m from the DUNE proton target and has a total active volume of 7 m wide, 3 m high, and 5 m long. Fiducialization reduces this to a fiducial mass of roughly 67.2 t [58,62]. The MPD is situated directly downstream of ArgonCube (designed to be a spectrometer of muons and other particles that do not stop in ArgonCube) with a cylindrical volume that is roughly 5 m in diameter and 5 m in height. This corresponds to an active mass of 1 ton. The MPD is situated inside an electromagnetic calorimeter and a magnetic field, allowing for precision measurement (and charge and particle identification) of the particles traveling through its fiducial volume. Fig. 1 provides a schematic drawing of the DUNE target and Near Detector Complex (note that many elements are removed from this figure for simplicity, including the magnetic focusing horns and a significant amount of earth between the decay volume near the target and the detector hall).
Given the dimensions of the detectors and the DUNE target/detector distance, we find that O(10 −2 ) of the axions produced via meson mixing will travel in the direction of the near detector complex. We include axion decays inside both the liquid and gaseous detectors in our simulations, corresponding to a total decay length of roughly 10 m. In Section IV, we discuss the various experimental signatures of this heavy axion decay and how we can reduce associated backgrounds in the two detectors.

IV. EXPERIMENTAL SIGNATURES OF HEAVY AXION DECAY
For the ALP masses we expect to be sensitive to at the DUNE ND, two classes of a decays are of interest: a → γγ and a → hadrons. For m a 1 GeV, the latter (a → hadrons) consists mostly of a → πππ and ππγ. Here, we highlight the characteristics of these respective signals in the DUNE ND complex focusing on both ArgonCube and MPD.
In Table I we list some defining characteristics of the two signals, a → γγ and a → hadrons, both in ArgonCube and MPD. We also list the types of neutrino-scattering backgrounds that contribute to these searches, and some properties of the backgrounds that allow for separating our signal from these events. For both signatures, and in both detectors, the a decay will be very forward-the large boost factor of a and the decay kinematics require this. Meanwhile, background events, such as those from neutral-current (NC) π 0 production (and subsequent π 0 → γγ decay) will be more isotropic, and may also have some measurable nuclear recoil that would not be present in the signal. In our signal a → γγ, since it is a fully visible decay, the invariant mass m 2 γγ = m 2 a , whereas the background events should reconstruct m 2 π 0 in this case. Additionally, the π 0 are being produced by neutrino scattering and will MeV, produced via meson mixing (with π 0 ) decaying in the DUNE Near Detectors, with respect to the total diphoton energy E γγ and opening angle between the photons ∆θ γγ for the decay channel a → γγ. The left panel assumes a is long-lived relative to the distance between the DUNE target and detector, whereas the right panel assumes it is short-lived. One-dimensional distributions for each of these observables are shown on top of/to the right of the two-dimensional distribution panels.
have at most E π 0 5 GeV (a conservative estimate). Lastly, the large boost of the a will result in small opening angles in the γγ final state, whereas the (less-boosted) π 0 from NC production will have larger opening angles.
We illustrate a subset of these distinctions in Fig. 7, where we show two different signal event distributions for a 200 MeV axion produced via meson mixing and decaying in the DUNE ND with the signal a → γγ. The main panels in Fig. 7 (left) and (right) display the distribution of these signal events as a function of the total diphoton energy E γγ as well as the angle between the two outgoing photons in the lab frame, ∆θ γγ . On the top (right) of the main panels, one-dimensional histograms display the distributions of E γγ (∆θ γγ ) independently. Both panels of Fig. 7 correspond to an axion with m a = 200 MeV produced from meson mixing, predominantly with π 0 . The distinction between these two panels is in the lifetime of a, i.e., whether it is long-lived (left) or short-lived (right), relative to the target-detector distance of 574 m. In the left panel, we assume a is long lived and cτ 574 m. Here, the probability of a given a to decay in the detector is proportional to (γcτ ) −1 which scales as m a /E a and favors lower-energy a from the production distribution. In contrast, the right panel assumes a is short-lived, cτ 574 m. In this scenario, only the high-energy a have γ high enough that their time-dilated lifetime is on the order of the target-detector distance and can survive that journey. This results in high-energy a being favored, which also implies very small opening angles in the diphoton system.
A recent study explored the capability of the gaseous argon MPD to search for decays of dark sector particles, including dark photons and dark Higgs bosons that can decay fully visibly into the final state e + e − [60]. This background channel has a decent degree of overlap with the a → γγ channel we are interested in because its dominant background is from the NCπ 0 production. 4 The searches in Ref. [60] involved lower-energy new-physics particles than those compared here, so the high energy of E γγ provides an additional mechanism to separate our signal from the NCπ 0 backgrounds. Therefore, we expect that a nearly background-free search for a → γγ is possible and will proceed under that assumption.
If we shift our focus to the hadronic final states, a → hadrons, the signal characteristics are not too different. Especially if we consider the a decay as the process a → hadrons, we can characterize the final state in terms of the total hadronic energy E had. and an opening angle ∆θ had. , which is a proxy for the total jet size of all of the final state hadrons. Signal distributions of this variety are shown in Fig. 8, where we now assume that m a = 1.2 GeV and that a is produced via the gluon-gluon fusion process discussed above. As with Fig. 7, we display the event distributions with respect to E had. and ∆θ had. . 5 The left (right) panel assumes that a is long-(short-)lived relative to the distance between the DUNE target and ND. This explains why lower energies are favored in the left panel and higher energies in the right one. We note that the one-dimensional E had. distributions in Fig. 8 on top of each panel nearly match the shapes of the histograms in Fig. 6 for the "at Production" and f G = 50 PeV choices. In Fig. 9, we show the number of expected signal events as a function of the axion decay constant f a for various axion mass points for the gluon-gluon fusion sample. For each mass point, we draw three curves of the same color to indicate the numerical uncertainties of our study. For low decay constant f a , the production rate is high, but the detection probability is exponentially suppressed by the detector distance D over the lab-frame decay probability, exp(−D/βγcτ ). For high decay constant f a , the lifetime is long, and the expected signal number is suppressed by 1/f 2 a for production rate and L/(βγcτ ) for the detection probability. Here L is the effective detector length in the line-of-flight for the axion. Due to the large spread of the axion boost factors at production, the transition between these two limits spreads over the decay constant over a decade or so. We can see, as anticipated, the expected number of signal events decreases for increasing mass due to the rate suppression. More importantly, a larger signal mass means a smaller boost, a shorter lifetime. To reach the DUNE ND, it requires larger f a to overcome the arrival flux suppression exp(−D/βγcτ ). Similar suppression exists for resonant mixing, as we show in yellow color when the axion mass is nearly degenerate with η meson. Nevertheless, thanks to the large flux at DUNE, the DUNE ND will be able to probe the high f a regime uniquely, as shown in the next section.
Backgrounds in this channel are mainly from scattering events that produce many final-state pions, etc., including charged-current scattering that produces a single charged muon and one or more pions (CC1µ2π). Deeply-inelastic-scattering (DIS) events, where the argon nucleus is completely broken up, can also result in events that would mimic this signal. However, as all of these events are generated by neutrino scattering, their total energy will (as in the NCπ 0 case) be less than roughly E had 5 GeV. Our signal events, as demonstrated by Fig. 8, will have hadronic energies 30 GeV, and even higher if the axion is short-lived. As with the a → γγ final state, the direction of these events is very forward-going, whereas the background will be more isotropic, and the opening angle is much smaller in the signal distributions than the backgrounds. With all of these features, we expect the a → hadrons search channel to be background free, like the a → γγ channel.
Before proceeding, we also wish to discuss one unique strength of the search at DUNE ND: combining searches for decaying heavy axions in both the liquid and gas detectors into one combined analysis. The background contributions discussed above are from beam neutrinos scattering in one of the detectors. These background rates scale with detector mass, and so the expected background contributions in the liquid detector are a factor of over 50 higher than in the gaseous detector. Meanwhile, the signal rate of decaying axions is more-or-less proportional to the volume of the detector and, therefore, will be roughly equal in the two detectors. A combined analysis, where the expected signal-to-background ratio can be robustly predicted from one detector to the other, can improve the overall DUNE ND capability.

V. DUNE NEAR DETECTOR SENSITIVITY TO HEAVY AXIONS
Combining all of the ingredients discussed to this point, we are now prepared to estimate the DUNE ND sensitivity to heavy axions.

A. Gluon Dominance Projections
We first focus on the case of gluon dominance discussed in Sec. II C where c 3 = 1, c 1 , c 2 = 0. 6 Combining both the meson mixing and gluon-gluon fusion production modes, we determine the parameter space for which we would expect three or more signal events in ten years of data collection at DUNE ND, the red shaded region in Fig. 10. The blue shaded regions correspond to the same ones shown in Fig. 2 based on theoretical considerations regarding the axion Quality Problem. The brown shaded regions are the same as in Fig. 2 as well. The horizontal shaded region labeled "Colored Particles" is disfavored since from a UV perspective, aGG coupling generically originate after integrating out colored fermions with masses ∼ yf a . Requiring a maximal Yukawa coupling y ∼ 4π along with the LHC constraints on colored states to have masses above 2 TeV [79][80][81][82] gives the above bound.  [72,76] and cosmology [77] are shown in dashed lines given the astrophysical uncertainties and model dependences. The region "Existing Constraints" include the bounds from partially invisible kaon decays from E787 and E949 [72], electron beam dump [31,33], CHARM [78], visible kaon decays [54], B decays [48], LHC dijet searches [43]. We also include projections relying on the proposed displaced track trigger at the HL-LHC [15], FASER [78] and NA62 [72].
We also include a number of existing experimental/observational constraints on this parameter space in grey. For small mass m a 100 MeV, astrophysical and cosmological constraints are relevant -the region labelled "SN1987A" indicates the region of parameter space for which such axions would cool the supernova and carry away too much energy [72,76], 7 whereas the region labelled "BBN + N eff " indicates where such axions would affect light-element abundances and contribute to the number of effective degrees of radiation [77,83]. Both the "SN1987A" bound and "BBN + N eff " are shown via dashed lines because of their associated uncertainties, see e.g. [84] and [77], respectively. For m a between 1 MeV and 1 GeV, a number of searches have been performed in the contexts of both electron and proton beam dumps and corresponding bounds were discussed in [31,33,78] for these types of axions, as well as searches for rare meson decays [48,49,54,72,85]. Furthermore, there are constraints from LHC dijet searches for m a > 50 GeV as obtained in [43].
Other planned experiments with similar timescales as DUNE are capable of performing searches in this region of parameter space. We include some projections of these in Fig. 10 as well. 8 The FASER [78] (dot-dashed cyan) experiment at the LHC will be able to probe a similar mass regime as DUNE but with smaller f G due to its close proximity/high-energy production source. A proposed displaced decay search using the high-luminosity LHC track trigger [15] can probe the region encompassed by the dotdashed brown line at heavier masses than DUNE. Finally, for m a < m π , NA62 has powerful sensitivity (dot-dashed purple) via the search for K + → π + + a where a is undetected [72]. Comparing our DUNE projections against these other future proposals, the complementarity of these different search strategies is obvious -the combination of all of these will allow for considerable reach in the theoretically-motivated parameter space in a way that no individual experiment can accomplish on its own. DUNE will specifically be most powerful in this long-lived region of small g agg or large f G , especially for 20 MeV m a 2 GeV.  [72,76] and cosmology [77] (appropriately recasted for our purpose) are shown in dashed lines given the astrophysical uncertainties and model dependences. The region "Existing Constraints" include the bounds from kaon decays from E787 and E949 [72,85], electron beam dump [31,33], as well as diphoton and dijet searches at the LHC [43]. We also include projections relying on the proposed displaced track trigger at the HL-LHC [15], NA62 [72] and diphoton searches at LHCb [51] and HL-LHC [43]. For the present "Codominance" scenario, we have not reanalyzed the coverage discussed in [54] for NA62, NA48, KOTO and in [85] for Belle-II. They would cover a mostly complementary region roughly for 1/f G > 10 −4 GeV −1 and m a few GeV. The exisiting CHARM data [86] would also cover some part of the parameter space for roughly 100 MeV < m a < 1 GeV (details in the text).

B. Codominance Projections
Here we focus on the scenario where c 1 = c 2 = c 3 and derive the coverage for 3 events for 10 year data taking at DUNE ND, shown in red in Fig. 11. While this coverage is similar to the one in Fig. 10, around m a ∼ 400 MeV, the effective photon coupling c γ , in Eq. (10) becomes small for our choices of c i . As a result, the coverage region shifts upwards exhibiting a "peak"-like feature.
The theoretical constraints from the axion Quality Problem remain the same. The set of experimental/observational constraints are shown in grey. Since the "SN1987A" constraints are dominated by the GG coupling, it remains the same. However, the "BBN+N eff " constraint is dominated by the FF coupling, and hence it gets modified based on Eq. (10) as a result of non-negligible c 1 and c 2 . The constraints from partially invisible kaon decays are also modified due to the non-negligible aWW coupling [49,72,85]. We recast the electron and proton beam dump results from [33] as appropriate for the present case of Codominance. Since the diphoton decay modes are now non-negligible, for higher masses m a > 20 GeV, both the diphoton and dijet searches at the LHC give relevant constraints [43].
Some comments regarding a few omissions in Fig. 11 are in order. We expect some part of the parameter space for 100 MeV < m a < 1 GeV would be covered by the existing CHARM data [86] which we have not derived for the Codominance scenario. Also, we have not derived the constraints and projections from KOTO and NA62/48 from visible kaon decays for this scenario, which would cover some parameter space for 150 MeV < m a < 350 MeV and roughly 1/f G > 10 −4 GeV −1 , mostly complementary to our DUNE ND coverage. In Ref. [54], such constraints were derived for the cases GG and WW -dominance separately. Some complementary coverage, again for roughly 1/f G > 10 −4 GeV −1 and m a few GeV, would also come from B → Ka processes at Belle-II similar to what is discussed in [85] for the case of WW -dominance.
To summarize, similar to the case of "Gluon Dominance" above, we see that DUNE ND would provide a powerful coverage, complementary to other existing and projected constraints, especially for large f G and 30 MeV m a 1 GeV.

VI. CONCLUSION
Recent studies of the Strong CP Problem (and the associated Axion Quality Problem) have led to a renewed interest in heavy axions with masses in the MeV-TeV regime. Meanwhile, a number of upcoming and planned experiments are capable of searching for decays of long-lived particles in a beam dump environment. One of the best example, in terms of the total protons on target (POT) and large, multipurpose detectors, of such an experiment is the Deep Underground Neutrino Experiment (DUNE) with its Near Detector (ND) complex. Combining the intense, high-energy proton beam (with a large number of POT per year) and the fine-grained NDs (both the liquid and the gaseous argon ones, allowing for particle identification and energy resolution) provides an exciting prospect for such searches.
In this paper, we have thoroughly explored the DUNE ND complex's ability to search for heavy axions in the MeV-GeV regime. We have revisited previous considerations of heavy axion production through both neutral, pseudoscalar mixing as well as through gluon-gluon fusion. Motivated by the Strong CP Problem, we have focused on two cases of these heavy axions via an Effective Field Theory treatment -one where the axion's dominant coupling is to the SM gluon field strength tensor, and one where it couples democratically to each of the SM gauge group field strength tensors. This is a different focus than the often-studied photon-dominant scenario for axion-like particle searches in beam dump environments.
The DUNE NDs offer several ways of identifying the decays of these heavy axions in their dominant decay channels, which are, depending on the axion mass, into photon pairs or hadrons. We have identified how these searches can leverage different signal characteristics to fully suppress neutrino-related backgrounds, allowing for very powerful searches of these rare signatures. Comparing to other projections for these classes of heavy axions, DUNE provides complementary sensitivity, specifically to very long-lived axions. Performing this type of search in tandem with other collider-based or meson-decay-based searches will allow us to cover as much of the theoretically-motivated parameter space as possible. There exist many more ways to explore these intriguing heavy axion theories at DUNE, including a large variety of production modes, from bremsstrahlung, meson decays mediated by operators beyond the gluon field strength, meson flavor changing decays, hadronic Primakov processes, as well as the rich decay channels from different Axion EFTs.
Whether or not an axion exists as a solution to the Strong CP Problem, as well as if it is in this heavy-axion category, remains to be seen. Regardless, experiments such as DUNE can perform unique searches for these and other new-physics scenarios without detracting from their overall scientific missions (in this case, neutrino oscillation studies). It is imperative that these searches are performed so that our planned experiments can extract as much scientific knowledge as they can. If such a heavy axion does exist within the reach of DUNE ND, then not only will DUNE revolutionize the field of neutrino physics, it will revolutionize our understanding of axions as well.  [57,78], MATHUSLA [55] and CODEX-b [56] for the Gluon Dominance scenario. The latter two projections on MATHUSLA and CODEX-b are taken from a recent analysis in Ref. [87] where the coverage are subject to correction (see details in the text).

Appendix A: Comparison Against Other Proposed Experiments
To place our projected limits on the gluon dominance scenario discussed in Section V A, here we provide a version of Fig. 10 with a larger set of future, proposed experimental sensitivities. In addition to the projections from NA62, FASER, and the HL-LHC shown and discussed in the main text, we include here projections from FASER2 [57,78], CODEX-b [56], and MATHUSLA [55]. Given the energies/detector locations of these different proposals, we see that, DUNE will still have unique sensitivity at large f G /small g agg as long as 20 MeV m a 2 GeV. These other proposals, specifically CODEX-b and MATHUSLA, provide sensitivity at higher m a in the same region of parameter space as the LHC Track Trigger proposal [15], 1 GeV m a 10 GeV, another interesting regime for heavy axion searches. We note here that the MATHUSLA and CODEX-b projections are taken from a recent analysis in Ref. [87], where new production modes from gluon splitting, gluon-gluon fusion, and meson decays are included. We found a typo in the gluon splitting function and that the reach in the O(few GeV) regimes are somewhat better than one would estimate. 9