Axion-like Particle Searches at DarkQuest

Axion-like particles (ALPs) interacting with the Standard Model can be abundantly produced in proton beam fixed-target experiments. Looking for their displaced decays is therefore an effective search strategy for ALPs with a mass in the MeV to GeV range. Focusing on the benchmark models where the ALP interacts dominantly with photons or gluons, we show that the proposed DarkQuest experiment at Fermilab will be able to test parameter space which has been previously inaccessible. We pay particular attention to the self-consistency of gluon-coupled ALP production and decay calculations, which has been recently shown to be a problem in many existing predictions. We also apply these results to explore existing constraints in the ALP parameter space.


Introduction
Axions and axion-like particles (ALPs) are a generic feature of many theories of beyond-Standard Model (BSM) physics. Originally introduced to resolve the strong CP problem [1][2][3][4], axions occur as pseudo-Nambu-Goldstone bosons (pNGBs) of spontaneously-broken global symmetries, or as zero-modes of higher-dimensional gauge fields [5][6][7][8]. The pNGB nature of axions ensures that they are technically natural, even if the physics associated with the global symmetry lies at a very high scale. As a result, these particles can be the first messengers of the ultraviolet (UV) that can be accessible at experiments.
The nature of UV physics (i.e., the gauge charges and couplings of the heavy BSM states) dictates the couplings of the ALP to SM particles. For example, the defining coupling of the QCD axion is its interaction with gluons which is generated by new particles with strong interactions. Depending on other charges of these heavy states, couplings to other SM gauge bosons are induced. In this minimal model the axion obtains a mass from the QCD condensate, and, as a result, its mass and coupling are fixed by a single parameter. Axion-like particles are a generalization of such a scenario where the mass and coupling can be varied individually.
The shift symmetry of axions also ensures that they interact with SM particles through dimension five operators suppressed by the high scale associated with the UV physics. Depending on the size of this scale and the mass of the ALPs, they can provide a cosmologically stable dark matter candidate [5,[9][10][11]. Alternatively, if their lifetime is short on cosmological time scales, one can produce and detect them in astrophysical systems or in terrestrial experiments. A wide range of observations has been carried out across decades in mass and coupling resulting in stringent constraints on the axion and axion-like particle parameter space -see, e.g., Refs. [12,13]. In this paper we study the sensitivity of the proposed high-intensity proton beam-dump experiment DarkQuest at Fermilab and show that it can probe previously inaccessible parameter space. Many other dark sector models have also been studied in the context of DarkQuest -see, e.g., Refs. [14][15][16][17][18][19][20][21].
Several collider and fixed target experiments with different beams and targets have searched for ALPs over the last 30+ years. The non-renormalizable nature of the ALP coupling naturally leads to long lifetimes on experimental length scales. As a result, if these particles are produced they can travel a macroscopic distance before decaying, giving rise to a displaced energy deposition in a detector downstream of the target. While certain regions of parameter space have been excluded by searching for this signature, large swaths remain untested. These are typically associated either with masses that are kinematically inaccessible, or with lifetimes that lead to decays in the shielding. The DarkQuest experiment provides a unique opportunity to access precisely these regions of parameter space with a compact set-up and a high-energy proton beam. While the proton beam also enables the study of many different couplings of ALPs, in this work we focus on the their interactions with gluons or photons.
Even though the hadronically-coupled ALP is an extensively studied model, its interactions at low energies (below the QCD phase transition) are still not fully understood due to the interplay with non-perturbative QCD dynamics. In particular, Refs. [22,23] have recently pointed out that many existing calculations of ALP reaction rates miss important contributions that are needed for theoretical consistency (to be made precise in the following sections). For some processes, this leads to rate changes of more than an order of magnitude [23]. Therefore, in our calculations we ensure that these theoretical consistency conditions are satisfied. This amounts to tracking the cancellation of unphysical parameters in physical amplitudes describing ALP production and decay. This paper is organized as follows. In Sec. 2 we describe the proposed DarkQuest experiment and discuss potential backgrounds for the search for ALP decays to photons a → γγ (this decay channel is generally present both in the dominantly photon-and gluoncoupled models). We then discuss ALP interactions below the QCD phase transition in Sec. 3. Following Refs. [22,23] we discuss the cancellation of unphysical parameters in chiral perturbation theory, extending some of their findings to three flavours and other interactions. We then use these results to study ALP production channels in proton beamdump experiments in Sec. 4 for photon-coupled ALPs and in Sec. 5 for gluon-coupled ALPs. 1 We apply these calculations in Sec. 6 to project the sensitivity of DarkQuest and to study existing constraints. Our main results for the photon coupling are summarized in Fig. 15 and for the gluon coupling in Fig. 17. Some constraints on the gluon-coupled model are re-evaluated in Fig. 18. We conclude in Sec. 7.

DarkQuest Experiment
DarkQuest is a proposed upgrade of the existing SeaQuest/SpinQuest experiment at Fermilab [24], which collides a 120 GeV proton beam with a thin target. The set-up of the experiment is shown in Fig. 1. The current spectrometer is optimized for measuring energetic muons emanating from the target. The target is followed by a 5 m magnetized iron block, the FMAG, which serves as a beam dump and also sweeps soft charged particles out of the detector acceptance. The FMAG is followed by a series of detector subsystems for tracking and analyzing charged particles; we will focus on ALP decays to photons and therefore will not make use of these.
The DarkQuest upgrade will add an ECAL at z ≈ 19 m and improve data acquisition systems, enabling sensitivity to electrons and photons. This version of the detector will make possible a multitude of searches for various dark sector particles [14][15][16][17][18][19][20][21]. The high beam energy and relatively compact geometry of the dump make this an ideal experiment to study production and decay of long-lived particles. In addition, the ECAL modules will be repurposed from the PHENIX experiment [25] allowing for rapid and cost-effective construction with a physics run possible as early as 2023. The Fermilab accelerator complex can provide 10 18 protons on target (N POT ) in about two years of running. We will refer to this benchmark as phase 1 of DarkQuest. Future upgrades of the complex in the on-going Proton Improvement Plan will enable even more intense beams [26]. We assess the gains from increased luminosity by considering a second phase with N POT = 10 20 .
Unlike searches for long-lived particle decays to charged particles, photon final states are subject to additional backgrounds which we describe below. These necessitate additional shielding after FMAG, or preshower detectors to enable tracking.

Backgrounds
Certain SM processes can mimic the displaced decay signal of ALPs. Below we consider two such processes: production of short-lived particles in the back of the dump and decays of SM long-lived particles.
Proton-nucleus interactions that take place in the back of the beam dump can produce particles that make it out to the decay volume. Production of π and η ( ) mesons and their decay into γγ can therefore potentially mimic our signal if it occurs deep enough in FMAG. The number of mesons produced in the last interaction length of the 5 m iron FMAG is where the first factor is the attenuation from requiring the beam particle to make it this far in the iron without an inelastic scattering, n M is the number of mesons produced per interaction (n π 0 ≈ 2.5 and n η ≈ 0.3). The photons from the decays of these mesons have a high acceptance rate and would correspond to 10 5 (10 7 ) background events in phase 1 (2). These backgrounds can be mitigated by placing an additional 14 − 17 interaction lengths of shielding before or after FMAG. Alternatively, if the segmentation of the ECAL allows a rough vertexing of the decay photons, events that originate from the end of FMAG can be rejected. Following Ref. [16] we will use z 7 m as our fiducial decay region, allowing up to two meters after FMAG for additional shielding or vertexing.
Standard model long-lived particles can also mimic the ALP displaced decay signal. In particular, neutral kaons K S,L and the Λ baryon can have lab-frame lifetimes comparable to the size of the DarkQuest experiment. These hadrons have decay modes into neutral final states that include π 0 's which can produce photons, e.g., K S → π 0 π 0 , K L → π 0 π 0 π 0 and Λ → nπ 0 . We use PYTHIA to estimate the multiplicity n H and typical boosts γ H of H = K 0 and Λ. The multiplicities are n K ≈ 0.36 (this number includes both K 0 and K 0 and is consistent with measurements of the combined K ± multiplicity at NA49 [27]), n Λ ≈ 0.13 (consistent with the recent measurements from NA61/SHINE [28]), while for mean boosts we find γ K ≈ 18 and γ Λ ≈ 17. We use these to estimate the number of hadrons that traverse the FMAG and decay into final states containing photons beyond: In computing the above we assumed that the hadrons are produced in the first interaction length of FMAG, and accounted for the attenuation of their flux from inelastic interactions and decays in the dump. We see that displaced K L and Λ decays can produce significant backgrounds for the ALP search. As above, the simplest solution is additional shielding beyond FMAG. The same amount of shielding would reduce the rates for these displaced decays to negligible levels in phase 1 and 2. Additionally event selections based on the number of photons detected in the ECAL can be used to eliminate processes with multiple π 0 's such as K L → π 0 π 0 π 0 .

ALP Interactions
The general ALP Lagrangian depends on the ultraviolet (UV) completion of the model [22,23]; it contains direct couplings to quarks, leptons and gauge bosons. Moreover, even if the UV completion only generates one of these couplings at a high scale, renormalization group (RG) evolution will generate a myriad of other interactions.
In this work we focus on two simplified ALP models in which the dominant coupling just above the QCD scale is either to photons or to gluons. 2 Following the notation of Ref. [22] we write the ALP Lagrangian as where m a,0 is the bare ALP mass (without including QCD effects), f is the ALP decay constant (with dimensions of energy) and c γγ and c GG are dimensionless constants, which are typically O(1). 3 For the QCD axion m a,0 = 0 and the mass comes entirely from QCD (we will not make this assumption). The benchmark models with c γγ = 0 or c GG = 0 have also been explored in the CERN PBC study [32]. When discussing the model with dominantly-photon interaction we will also use the notation g aγ = c γγ α/(πf ).
The ALP-gluon coupling leads to interactions with hadrons below the QCD phase transition. These interactions can be derived within chiral perturbation theory (χPT), as we describe in the following section. As first pointed in Ref. [22], this procedure is fraught with the possibility of missing quantitatively important terms. We therefore pay particular attention to ensure that final results are physical (in the sense to be made precise below).

Hadronic Interactions at Low Energies
We now discuss the interactions of ALPs in the c γγ = 0, c GG = 0 model. The ALPs that can be discovered at DarkQuest have masses below about a GeV, where the appropriate description of ALP couplings in Eq. 3.1 is in terms hadrons rather than gluons due to confinement. In order to derive this description it is convenient to perform a field redefinition on the quarks

2)
2 This simplifying assumption ensures that we do not generate other important interactions through RG evolution. Such couplings can induce many additional signals and constraints -see Ref. [29]. 3 Parametrically larger coefficients can be obtained in non-minimal constructions such as clockwork [30], or through mixing of multiple ALPs [31].
where κ is a matrix in flavour space (we will consider both two and three-flavour limits). The rotated Lagrangian contains [22,29] where the hatted quantities include the effects of the chiral rotation: The first term in Eq. 3.4 is the original gluon coupling combined with a contribution that arises due to the chiral anomaly; the second term is the induced coupling to photons that also arises due to the chiral anomaly; the third term comes from the quark kinetic term, and the last term is the modified quark mass term. At this point the matrix κ is arbitrary, but we can see that if we choose tr κ = 1, (3.6) the gluon term is eliminated. One also usually chooses κ to be diagonal in the same basis as m q , i.e., κ = diag(κ u , κ d ).
The key point emphasized in Refs. [22,23] is that κ is completely arbitrary up to the condition in Eq. 3.6; as a result, physical observables cannot depend on κ. Cancellation of these unphysical parameters in Feynman diagram calculations therefore provides a non-trivial check of the results. This was test was carried out for several processes involving ALPs in [22,23].
We now derive the interactions with mesons and nucleons starting from Eq. 3.4.

Interactions with Mesons and Photons
We focus on the two flavour theory for simplicity first. We will have to study three-flavour case in order to incorporate interactions with η and η mesons. The second line of Eq. 3.4 becomes the following leading-order chiral Lagrangian: The first line is the ALP kinetic and mass terms, the second term is the meson kinetic and mass terms (m † q (a) is the quark mass matrix with the two exponential factors in Eq. 3.4 but without the γ 5 ); the third line comes from matching the axial current in the quark theory to the axial current in the meson theory. 4 In the two-flavour theory (with a decoupled η ), the matrix meson field is and we use f π ≈ 93 MeV [33,34]. The three-flavour model is considered in Appendix A. Note that our normalization conventions differ slightly from Ref. [22,23]: our f π is a factor of √ 2 smaller than theirs, which accounts for various differences between their formulas and ours.
The Lagrangian in Eq. 3.7 contains off-diagonal terms involving the ALP and π 0 (and η ( ) in the three-flavour theory): This mixing can be accounted for either perturbatively for each process involving ALPs, or by performing a complete diagonalization of the meson-ALP kinetic and mass terms at the Lagrangian level. The former provides a simple explicit illustration of how unphysical κ dependence cancels; the latter is more straightforward to use in the three-flavour case discussed in Appendix A.
In the perturbative approach, a generic ALP production process can be schematically represented as There are two contributions: one from direct ALP production and a second one from π 0 mixing with an ALP; both are needed to demonstrate the cancellation of the unphysical κ parameters. The mixing contribution to the processes requires the following Feynman rule which follows from Eq. 3.10 where the first term arises from the kinetic mixing, and the second term is the mass mixing.
An immediate application of this is to derive the physical ALP-photon amplitude Ref. [22]: parametrizes isospin breaking due the quark mass matrix. In the last line of Eq. 3.14 we used Eq. 3.5 and Eq. 3.6; note that the result is now independent of κ, so it gives the physical coupling of ALPs to photons. The equivalent three-flavour result that accounts for mixing with η and η mesons is given in Eq. B.1.
Now we use the chiral Lagrangian to derive amplitudes relevant for low-mass ALP decays to photons and mesons, as well as ALP production in meson decays. Some of these results are known (in particular, see Refs. [22,23,29]) in the two-flavour χPT. We extend these calculations to three flavours when relevant. For ALPs near the ρ mass and above, virtual vector mesons can play an important role in determining amplitudes of certain decay channels [35,36], which are not included in our calculation. It would be interesting to extend their calculations while keeping track of κ in pseudoscalar and vector meson interactions. This can be done within the hidden local symmetry framework [37], for example.
The amplitudes that are relevant for sub-GeV ALP decays are • Three-pion final states that follow from the chiral Lagrangian In the above expressions we have decoupled the singlet meson, but still worked in three flavours. In our numerical results we retain η − η mixing with their physical masses. Note that the two-flavour result of Ref. [22] is obtained by taking m K → ∞ in the two amplitudes above; in this limit the 3π 0 amplitude vanishes if m a is small, but this is not so in the three-flavour model.
In the mass window where a → ηππ is allowed it is clearly not appropriate to decouple the singlet meson. For these and other three-flavour calculations we work in the simplified η − η mixing scheme (in which we approximate sin θ ηη ≈ −1/3) unless stated otherwise; this produces compact expressions for the amplitudes. 5 The a → ηπ 0 π 0 amplitude is To leading order in isospin violation the amplitude to charged pions is the same The ALP mass scale for which these amplitudes are relevant is in the regime where standard χPT becomes unreliable. In order to improve the accuracy of our estimates, we follow Ref. [35] and use "k-factors" which rescale the resulting partial widths. The k-factors are obtained by comparing our calculations of the η → ηππ amplitudes to measurements [38].
The ALP decay calculations are summarized in Fig. 2 in which we show the ALP decay length and branching fractions to various final states. We note that for m a 1 GeV the partial width to γγ remains appreciable in most of the parameter space, making this a good decay channel to search for even in the dominantly-gluon-coupled ALP model; other decays result in more complex final states but can be used to recover some signal rate if a → γγ is suppressed.
• The π ± and kaon rare decays arise from electroweak interactions. We focus here on kaon decays because the ALP mass range relevant for π ± is well-covered experimentally (the amplitude for π ± → νµ ± a is given in Ref. [29]). Following the notation of Ref. [41], the EW Lagrangian relevant for s → d transitions after the chiral rotation, Linearising the right-handed chiral rotation R allows one to extract terms relevant for various kaon decays. Both of the ALP terms appearing in R and D µ are necessary for κ-dependence to cancel. We find the following amplitudes: where K ≈ 2.23 × 10 −3 is the CP-violating kaon mixing parameter. The charged kaon decay was already presented in Ref. [29] in the m η → ∞ limit. Here we give the amplitude in the same limit, but use the full result in our numerics: (3.31) When the ALP mass is close the masses of the neutral mesons, mixing is enhanced; these regions are highlighted by the vertical gray bands. In the right panel, the lines labelled a → 3π and a → ηππ combine partial widths to charged and neutral pions.
These amplitudes are given in the limit of decoupled singlet meson and to leading order in isospin violation for brevity. Their calculation, as well as subleading isospinviolating terms are discussed in Appendix C. Note that the leading-isospin amplitudes satisfy |A(K S → π 0 a)| ≈ |A(K ± → π ± a)|.

Interactions with Nucleons
Next we consider the nucleon-ALP interaction. In the low energy limit (such as that relevant for DM axion detection or astrophysical processes) it is enough to consider twoflavour χPT -these results appear in many places, including Refs. [23,42]. In particular, Ref. [23] has demonstrated the κ-independence of this interaction. We first repeat their analysis to illustrate this cancellation of unphysical parameters and then extend it to the three-flavour model. We will find that mixing with η and η becomes very important near the GeV-scale. The two-and three-flavour results for the effective proton interaction g pa are compared in Fig. 3.

Two-flavour Model
The leading pion-nucleon Lagrangian is [33,34] where Ψ = (p, n) is the nucleon doublet. The restrictions to "isovector" or "isoscalar" parts will be explained shortly. The nucleons transform non-linearly under SU (2) L × SU (2) R but simply under the isospin subgroup. The vielbein u µ | vec is given by where ξ 2 = Σ, r µ and l µ are right-and left-handed external currents, and is the axial current. 7 In the above we have isolated contributions from isovector and isoscalar axial currents; the isoscalar pieces receive contributions from η and ALP interactions, but in the two-flavour limit the η is decoupled, so Decomposing the coupling matrixĉ qq into isovector and isoscalar components gives: The leading ALP-proton interactions are therefore Note that the isoscalar piece is proportional to tr κ making them physical (i.e., no meson mixing contributions are required to cancel off unphysical κ dependence).
As with the physical ALP photon amplitude, the amplitude for p → pa involves mixing with π 0 as well as a direct coupling [23]; their sum is κ-independent: where we made use of Eqs. 3.6 and 3.13.
Three-flavour model Production of ALPs near the GeV scale is sensitive to mixing with heavier mesons, so it is useful to include η, η couplings to baryons, which necessitates we study the nucleon octet B. In order to include the singlet η we construct a U (3) L × U (3) R invariant Lagrangian following [43,44] We will see that linear combinations of the coupling constants D, F and D s are related to g A and g 0 from the two-flavour formalism. For example, if we forget about the ALP the leading interactions are Comparing this with the two-flavour result, Eq. 3.37, we therefore expect that g A ≈ D + F , and g 0 = D + 2D s + F . Refs. [45,46] found that D ≈ 0.80, F ≈ 0.46 and Ref. [44] fit D s in the range [−0.6, −0.2] (but they do not provide a best fit value). These values are consistent with three-flavour extractions of g A = 1.283 and g 0 = 0.384 from lattice QCD [47,48] if we take D s ≈ −0.41. The additional ALP term that results from Eq. 3.39 is In rotating to the physical basis, the physical ALP-proton coupling receives contributions from both Eq. 3.40 and 3.41; we write the result as where the quantities ·a are defined in Appendix A. By inserting explicit expressions for these, one can check for the cancellation of κ dependence. One can also take various limits, like decoupling the third flavour, or taking the mass of the singlet m η 1 → ∞ to recover Eq. 3.38. As before, working in the "simplified η − η mixing" scheme with sin θ ηη = −1/3 gives The result is reassuringly κ-independent and retains various enhancements in mixing with π, η and η when the ALP mass is near the corresponding meson mass. We compare g pa in two-and three-flavour theories in Fig. 3 (the two-flavour g pa is just the quantity in the parentheses of Eq. 3.38). We see important differences at large ALP masses, where the coupling is suppressed for m a > m η (although at these large masses one should start including higher meson resonances). At low masses the two results are not quite equal because m π /m η ( ) = 0.
In the discussion above we took the hadrons as point-like. In realistic calculations we must introduce form-factors to account for their extended size and substructure, which we will discuss in Sec. 5.

ALP Production and Signals From the Photon Coupling
In this section we focus on the photon-coupled ALP model and set c GG = 0. We mostly follow the notation of Ref. [49] and write the interaction as where The interaction in Eq. 4.1 enables the decay of ALPs into photons with width and a variety of production channels in proton fixed target experiments. The most important of these are Primakoff production γA → aA (with γ a secondary photon from, e.g., meson decay) and quasi-elastic "photon fusion" processes pA → apA, which are illustrated in Fig. 4. We discuss these in more detail below. While the large flux of secondary photons typically dominates the production of ALPs, photon fusion can produce ALPs with a  Z × pp (LUXqed) Figure 5. Cross-sections for various photon-coupled ALP production channels in collisions of a 120 GeV proton beam with an iron target. Primakoff (γA) and coherent photon fusion (pA) processes are used to estimate the sensitivity of DarkQuest. The latter channel is restricted to proton momentum transfers of Q 2 ≤ 1 GeV 2 for computational simplicity. Cross sections for processes with higher momentum transfers are obtained using the LUXqed photon PDF, but end up being subdominant in the parameter space accessible at DarkQuest.
slightly higher boost, enabling sensitivity to shorter lifetimes. For both processes the production is coherent over the nucleus for ALP masses GeV. Incoherent processes (such as analogues of deep inelastic scattering) can in principle produce more massive ALPs, but their cross-sections are typically very suppressed as we discuss below. In Fig. 5 we compare the cross-sections for different ALP production processes in the DarkQuest beam and target configuration. We describe their calculation next.

Primakoff Production: γA → aA
Proton-nucleus collisions produce a large multiplicity of neutral mesons such as π 0 and η which decay to photons, giving rise an intense secondary photon flux. These photons can undergo a Primakoff-type interaction with a nucleus in the dump producing an ALP as shown in the left panel of Fig. 4. This mechanism has been identified as the dominant production mode for photon-coupled ALPs in proton beam dump [16,18] and (primary) photon beam experiments [50] across a broad range of parameter space. Below we describe how we simulate this production channel at the DarkQuest experiment.
First, we need to estimate the rate, spectrum and angular distributions of photons produced in proton-nucleus collisions. Following Refs. [16,18] we use PYTHIA 8.240 to model meson production and their subsequent decay into photons using the SoftQCD flag. Ref. [18] has validated the meson multiplicity and angular distributions at several energies relevant for proton beam dump experiments, finding reasonable agreement of simulated transverse momentum and small longitudinal momentum distributions with data. However PYTHIA underestimates meson production at larger p z /p max z . These highly-boosted mesons produce the most energetic photons, and the most boosted ALPs (see Eq. 4.9 below). Since our simulation underestimates this region of phase space, we underestimate the sensitivity to shorter lifetimes and heavier ALPs. We leave an improved analysis utilizing more realistic distributions (such as those obtained from Ref. [51]) for future work. Since PYTHIA simulates pp collisions, we rescale the total cross-section to the experimentally fitted pA value where A is the nuclear mass number and we used the fit from Fig. 1 of Ref. [52] (based on data of Ref. [53]). 8 The PYTHIA simulation produces an average of approximately 6 photons per pp interaction with a mean energy of ∼ 3 GeV; the typical transverse momentum of these photons is given by the QCD scale ∼ 0.2 GeV.
The secondary photons produced in pA collisions can interact the nuclei in the dump to produce ALPs. The differential cross-section for Primakoff production follows from the interaction in Eq. 4.1 which is in agreement with the expression in Ref. [50] up to a different parametrization of the photon-nucleus interaction (this result is also equivalent to the expressions in Refs. [18,49,54] in the limit of m a , |t| E a , m A ). Here s and t are Mandelstam variables with s = (p γ + p A ) 2 , t = (p γ − p a ) 2 , with t 0,1 being the kinematic boundaries of t given in Eq. 4.7 below; Z and m N are the target nucleus charge and mass. The electromagnetic form-factor F is modelled using the Helm parametrization [55] wheres = 0.9 fm and This parametrization is also used in Refs. [18,49]. Finally, note that the kinematically- and For a given incident photon, the ALP energy is given by For small ALP masses, the distribution in Eq. 4.4 peaks at t ∼ −m 4 a /(2E 2 γ ), so from Eq. 4.9 it is clear that typically ALPs inherit most of the incident photon energy.
We use the secondary photons from PYTHIA to produce a sample of ALP events by drawing from the differential distribution in Eq. 4.4. These ALPs are displaced and decayed into photons. We then estimate the total event yield from this mechanism via where the outer sum is over the secondary photons and the inner sum is over the ALP momenta; s γA is the Mandelstam variable for the γA system; the function C implements ALP decay to γγ, and computes the probability weight associated with the various experimental cuts on the location of the vertex and properties of the photons which are described in Sec. 6 for DarkQuest; N mc,pA is the number of pA collisions simulated in PYTHIA and N mc,γA is the number of ALP events generated per secondary photon. We will only consider secondary photon production in the first interaction length of FMAG (T (p) = 16.77 cm) and ALP production in the first radiation length (T (γ) = 1.757 cm); this simplification allows us to neglect attenuation of the initial proton and secondary photon beams. For ease of discussion below, we will separate the event yield into a total cross section σ Prim and a dimensionless acceptance A Prim and A Prim follows from these two equations.
Representative kinematic distributions for ALP decay products for the Primakoff production mechanism are shown in Figs. 6 and 7, while the total cross-section is compared to other processes in Fig. 5.
The general features of ALP production and decay via the Primakoff process are evident in the distributions in the left column of Fig. 6. Since the γA interaction peaks at small momentum transfers (see Eq. 4.4 and Eq. 4.9), it follows that ALPs and their decay products roughly inherit the properties of the secondary photons that produce them. As a result, the typical energies and angles are similar for different ALP masses up to threshold effects at small photon energies.

Photon Fusion
Axion-like particles can also be produced directly in the fusion of two virtual photons from the incoming proton and the nucleus as illustrated in the right panel of Fig. 4. Several methods for evaluating this production rate have been considered before [49,56]. Here we implement a version of the equivalent photon approximation in which we treat the incoming proton as a beam of nearly-on-shell photons [49], enabling the following factorization of the production cross-section [57]: where dσ γA /dt 2 is the differential cross-section in Eq. 4.4 (we have relabelled t → t 2 to emphasize that t 2 is the Mandelstam parameter of the 2 → 2 sub-process γ * A → aA, not of the actual 2 → 3 reaction); dn γ is the (proton) virtual photon spectrum given by [57] where Q 2 = −q 2 and q = (ω, q) denotes the four-momentum of the virtual photon, x = ω/E p is the fraction of proton beam energy the photon carries away. The momentum transfer Q 2 can be related to the magnitude of the transverse momentum q t of the virtual photon: where we took the limit m p /E p 1. We see from this relation that Q 2 min ≈ m 2 p /(1 − x). The form-factors D and C encode the non-point-like nature of the proton; these are fit to electron-proton scattering data and are given by where with µ 2 p = 7.78, Q 2 0 = 0.71 GeV 2 . These "dipole" form-factors are valid at the ∼ 10% level for momentum transfers Q 2 1 GeV 2 [58].
The factorization in Eq. 4.13 is enabled by several approximations discussed in general in Ref. [57] and scrutinized in the context of ALP production in Ref. [56]. For example, the use of the two body sub-process cross-section in Eq. 4.4 assumes that Q 2 m 2 a , an approximation that breaks down for small ALP masses. Another simplification used above was taking m p /E p 1. The quality of these approximations depends on the ALP mass. For the DarkQuest beam energy and the mass range of interest our calculation typically overestimates the exact photon fusion rate by a factor of at most ∼ 1.5 − 2 [56] (our method is termed "photon absorption" in this paper). We will show that photon fusion is sub-dominant to Primakoff production in rate by about an order of magnitude, so we will not pursue refinements of these approximations.
We generate events based on the distribution of Eq. 4.13 using vegas [59,60] to perform importance sampling on x, q 2 and t 2 ; these are transformed into an ALP four momentum which is then decayed into γγ. We use these events to compute the total event yield, cross section and acceptance in analogy to Eqs. 4.10 and 4.11, 4.12 for the Primakoff process.
Various kinematic distributions for ALP decay products for the photon fusion production mechanism are shown in Figs. 6 and 7. The acceptance function and total cross section are shown in Figs. 5 and 12. Compared to the Primakoff process, heavier ALPs produced in photon fusion tend to be more boosted. This results in decay photons that are more forward and energetic as shown in the right column of Fig. 6. In Fig. 7 we show that this effect leads to more collimated photons for larger ALP masses; however these are still typically well separated enough.

Other Channels
The production channels discussed above feature coherent scattering of protons or secondary photons off the iron nuclei in the dump, such that their rates are enhanced by Z 2 . Production of heavier ALPs leads to larger momentum transfers and loss of this coherence. As a result processes like γp → a + X have a much smaller rate; moreover since the initial secondary photon is the same as in coherent γA reactions, this channel does not give any new kinematic reach. Naively, larger masses should be accessible in pA or pp interactions, but we have restricted ourselves to momentum transfers Q 2 ≤ 1 GeV 2 in Sec. 4.2 so that we could use simple dipole expressions for the proton form factors. At larger momentum transfers one must consider inelastic contributions and generation of photons from parton evolution of the proton constituents. These effects have been included in the LUXqed photon parton distribution function (PDF) [61,62]. We use LUXqed to estimate the cross-sections for pA → aA + X with Q 2 ≥ 10 GeV (the minimum momentum transfer allowed in that PDF) and pp → a + X. The former is calculated using MadGraph [63], while the latter is computed directly using LUXqed and the narrow width approximation via where f γ (x, Q 2 ) is the photon PDFs evaluated using lhapdf [64]. The results are shown in Fig. 5, where we scaled the pp cross section by Z = 26. We see that the ALP production rates are indeed dominated by Primakoff (γA) and coherent photon fusion processes (pA with Q 2 ≤ 1 GeV 2 ) for m a below a few GeV. We will see that the couplings g aγ required for ALP events to fall into the DarkQuest acceptance are 10 −6 GeV −1 for m a 1GeV; for such tiny g aγ we expect less than one ALP to be produced even in phase 2 from hard pA and pp interactions. We therefore neglect these processes in our analysis.

ALP Production and Signals From the Gluon Coupling
Gluon-coupled ALPs feature a wide array of possible production and decay channels. In addition to the Primakoff and photon fusion ALP production (which are enabled by the induced ALP-photon coupling discussed in Sec. 3.1.1), rare meson decay and other hadronic processes are allowed. We will show that the latter dominate because they are not suppressed by the electromagnetic coupling. Two representative processes are shown schematically in Fig. 8.
It is not obvious how to best model ALP production in hadronic interactions; various approaches have been used including emission in a parton shower [65] or hadronization [65,66]. While using the gluon coupling directly in a parton shower is manifestly κ-invariant, allowing the production of ALPs in hadronization by mixing with mesons can lead to κdependent results. This is because this is usually done by replacing neutral mesons from a Monte Carlo simulation of pp collisions by an ALP and reweighing the cross-section by the η (′) π π a A p a Figure 8. Representative production mechanisms of gluon-coupled axion-like particles in proton beam-dump experiments. In the left panel, an ALP is produced in a rare meson decay (other production channels of this kind include rare pion and kaon decays). The right panel shows the bremsstrahlung of an ALP from an incoming proton beam that undergoes a scattering off a nucleus.  Figure 9. Cross-sections of various gluon-coupled ALP production processes in collisions of a 120 GeV proton beam on an iron target. For rare meson decays, the cross-section is the meson production cross-section times its branching fraction into ALPs. The line labelled A × gg refers to the gluon-gluon fusion process in proton-nucleon collisions scaled by the atomic number of the target.
corresponding (κ-dependent) mixing. We therefore consider a simpler alternative to clearly track κ dependence: we instead compute ALP bremsstrahlung from a proton, following the recent results of Ref. [67]. It would be interesting to compare these different methods while ensuring cancellation of all unphysical parameters.
In Fig. 9 we show the cross-sections for the DarkQuest beam and target configuration. We see that below a GeV the ALP production rate will be dominated by rare meson decays and proton bremsstrahlung. 9 In this plot, as in all other rate calculations we assume that the meson production occurs in the first interaction length of the iron dump.   Figure 11. Distributions of angular separation between the photons in a → γγ for two gluoncoupled ALP production mechanisms (rare η decays and proton bremsstrahlung) for m a = 0.5 GeV. The dotted line indicates the approximate minimum photon separation assuming the decay happens at z = 8 m and the detector is at z = 19 m. The sharp edge in the bremsstrahlung histogram is due to the limited range of validity of the "quasi-real" approximation (see Sec. 5.2) which prevents the ALP from carrying away too much of the beam momentum (and therefore cuts off the very collimated part of the phase space).

Rare Meson Decays: M → a + X
ALPs can be abundantly produced in rare decays of π ± , η, η , K L,S and K ± mesons. We do not consider the π ± production channel here only because that parameter space is well covered by existing searches. As in Sec. 4, we use PYTHIA 8.240 to model meson production rate and kinematics. We find that the number η ( ) and K mesons produced per pp interaction is n η ≈ 0.30, n η ≈ 0.034 (5.1) and The η multiplicity agrees with measurements of Refs. [68,69] (albeit at different beam energies); the η multiplicity is consistent with scaling the η rate by sin 2 θ ηη as a naive estimate. The kaon multiplicities qualitatively match the results of Ref. [27] (see their Fig.  130).
In order to avoid dealing with the attenuation of the proton beam we will focus on collisions in the first nuclear collision length of the target. Given a total number of protons on target, N POT , the total number of mesons produced is then where S is a symmetry factor (= 1 if the final state mesons are distinguishable, and = 2 if they are not); β i are with β(x, y) = 1 − 2(x + y) + (x − y) 2 . (5.8) The partial widths for the two-body decays of kaons are simply where |p a | is the magnitude of the ALP three-momentum We simulate ALP production by sampling the integrands in the above expressions, reconstructing the full ALP four-vector and boosting it to the lab frame specified by the meson four-vector from PYTHIA. The ALPs are displaced and decayed to photons. A typical kinematic distribution of the daughter photons from η -produced ALPs is shown in the left panel of Fig. 10 (other mesons give similar distributions). The angular photon separation is shown in Fig. 11.

Proton Bremsstrahlung: pA → a + X
In this section we consider a model for ALP bremsstrahlung off protons based on the results of Ref. [67] (see also Ref. [70]); a similar calculation of gluon-coupled pseudoscalar production geared towards LHC/Forward Physics Facility energies can be found in Ref. [71].
A key advantage of this process is that it allows us to simply use the physical, κ-independent ALP-proton interactions obtained in Sec. 3.1.2. The idea is to express the rate for pp → X + a in terms of the cross-section for a well measured process like pp → X times a splitting function that encodes ALP radiation (we will generalize the discussion to protonnucleus collisions at the end of this subsection). The pp cross-section is dominated by low-momentum transfer reactions involving non-perturbative physics that can be modelled using Pomeron and meson exchange. Since we are considering relatively light ALPs, pp → X + a should still be dominated by similar processes. The pp processes can be divided into elastic non-diffractive (both protons remain intact), single-diffractive (SD, one proton is dissociated, producing a multiparticle hadronic state X) or double-diffractive (DD, both protons dissociate). It was emphasized in Ref. [67] that if the underlying reaction is elastic or SD, the radiation of any light state from initial and final state protons interferes, leading to a severe cancellations in the emission rate. 10 As a result, we expect the production rate to be dominated by non-elastic, non-single-diffractive events (NSD) in which the initial state radiation of the ALP from the beam proton (emission from the target is suppressed -see Ref. [67]) does not interfere with other possible contributions.
Under certain kinematic conditions (roughly small momentum transfers, i.e. large beam energy, small ALP mass and forward production) ALP bremsstrahlung can be factorized from the underlying pp scattering, an approximation termed "quasi-real" in Ref. [67] because the slightly-off-shell intermediate proton is approximated as on-shell. We summarize the main results here, with a detailed derivation given in Appendix E. We start with the following interaction L ⊃ g pa (∂ µ a)pγ µ γ 5 p (5.11) where g pa is computed in Sec. 3.1.2. Under certain conditions the spin-summed and averaged matrix element can be written as 12) where z is momentum fraction of the initial beam carried away by the ALP, p T is the ALP transverse momentum, H = p 2 T + z 2 m 2 p + (1 − z)m 2 a is related to the invariant mass of the intermediate (off-shell) proton and A(pp → X) is the amplitude for the underlying hadronic interaction. The differential cross-section can then be written as where the ALP splitting function w a is and σ(s ) is the cross-section of the underlying hadronic process pp → X, evaluated at a slightly different center-of-mass energy with s ≈ 2m p p p (1−z) with p p the beam momentum. This expression is not yet complete as we need to dress it with form-factors to account for the non-point-like nature of the beam particle; we discuss this in the following subsection. The factorization of ALP radiation in Eq. 5.13 requires several assumptions that limit its range of validity; these assumptions are detailed in Ref. [67] and discussed in Appendix E. Their physical content is to require the ALP, beam and recoil proton to be ultrarelativistic. This limits the range of z, p T and m a one can consider in this approximation. These conditions are achieved when the beam energy is much larger than other energy scales in the process; this is a good assumption for the DarkQuest configuration with a 120 GeV beam and a forward detector (limiting the range of relevant p T ) searching for ALPs with m a GeV.
The above discussion followed closely Ref. [67] and focused on pp collisions. However, it is easy to translate it to proton-nucleus collisions which are of more direct relevance for DarkQuest. The only modification is to replace the pp NSD cross-section σ(s ) by an equivalent quantity for pA scattering. In order to avoid the initial/final-state radiation interference discussed in [67] we focus on processes where the beam proton, or both beam proton and target nucleus are disrupted. This corresponds to any inelastic scattering process except for the ones where only the target nucleus is excited (termed target single diffractive, TSD, in [52]). We use fits to data for inelastic and TSD cross-sections from Ref. [52]. The resulting cross-section is At the relatively low √ s 15 GeV of interest, these cross-sections are very weakly dependent on s (see, e.g., Ref. [72]) and we treat them as constant.
We simulate ALP production and compute the rate by sampling Eq. 5.13 (with the cross-section in Eq. 5.15 and the form-factors discussed in the next subsection) using vegas [59,60] and reconstructing kinematics using Eq. E.1. The resulting bremsstrahlung cross-section is compared other processes in Fig. 9. The ALPs are displaced and decayed to photons. A typical kinematic distribution of the daughter photons is shown in the right panel of Fig. 10. The angular photon separation is shown in Fig. 11. We note that bremsstrahlung is able to produce heavier ALPs than rare meson decays, and these ALPs tend to be much more boosted. As a result, the decay photons are more collimated.

Form Factors
The previous calculation assumed that the protons and mesons are point-like. In order to account for their finite size and internal structure we include form factors following Ref. [67]. This means that each ALP-proton vertex should be multiplied by a scalar function of momenta; this can depend on any Lorentz invariant quantity relevant for the vertex, but momentum conservation implies that without loss of generality we can take where p and p − k are the proton momenta before and after ALP radiation, and k is the ALP momentum. In most discussions of proton form factors one assumes that the protons are on-shell, in which case the only non-trivial argument, k, is space-like. We are, however, interested in the situation where p is slightly off-shell and k 2 = m a > 0. There is no data available for such a time-like form-factor, so we must construct a model. Following Ref. [67], we take the incoming proton to be on-shell, and write the total form factor as We will refer to F 1 as the time-like form factor; F pp * is inserted to control the factorization of the ALP bremsstrahlung cross-section into the form of Eq. 5.13, so it is chosen to be equal to 1 when p − k is on-shell and falls off as it becomes off-shell. As in Ref. [67] we use where Λ is an unknown parameter O(GeV). We will take Λ = 1 GeV. This form-factor has been successfully used to fit a variety of experimental data involving nucleon-nucleon-meson interactions [73,74] with Λ ≈ 1 GeV; while these datasets have √ s 2 GeV, the typical intermediate off-shell nucleon invariant mass squared is just s for their 2 → 2 reactions. For the bremsstrahlung process, (p − k) 2 − m 2 p = −H/z where H is given below Eq. 5.12; the form-factor therefore suppresses soft and high p T ALP emission, both regimes where the "quasi-real" factorization is spoiled [67]. It would be interesting to fully validate this approach in a kinematic region similar to DarkQuest, with √ s ∼ 15 GeV; for example, the same approach can be used to study the forward emission of single π 0 or η. Unfortunately we are not aware of such a data set.
We now turn to the time-like form factor F 1 . There are few studies on axial vector time-like form factors and no data in the relevant kinematic region. Ref. [75] and references therein describe how to analytically continue a space-like axial vector form-factor to the time-like region. In order to do this one needs to know the singularity structure of the form-factor; this is not captured by the usual "dipole" form factors that are fit or data, or to lattice results (see, e.g., [47,48]). Thus, Ref. [75] uses a more physical "two-component" model in the space-like region and analytically continues that. This model is The first factor is meant to describe the "intrinsic structure" of the nucleon consisting of three valence quarks. The second factor describes the contribution of resonances. We are specializing to the axial vector isovector current coupling for now; the lightest meson with the right quantum numbers I G (J P C ) = 1 − (1 ++ ) is the a 1 (1260) with m a 1 = 1.23 GeV and Γ a 1 = 400 MeV. γ and α are parameters, while m a 1 is fixed. γ is fit to electromagnetic scattering data with the result γ = 0.515 GeV −2 , leaving only α to fit to axial data in the space-like region. Ref. [75] obtained α = 0.95 for a joint fit to multiple data sets with −k 2 2 GeV 2 . The authors then proceed to extend this two-component model into the time-like domain; since there are singularities along k 2 > 0, phases arise, requiring the introduction of an additional parameter, δ: (5.20) δ = 0.397 was found in Ref. [76]. We emphasize that this form-factor is an extrapolation into the time-like region, where there is no data; some proposed measurements are discussed in Ref. [75].
The previous discussion focused on the isovector axial form-factor. However, in the physical basis the ALP-proton coupling has both axial isovector and isoscalar pieces (this is easily seen in Eq. 3.37). We have not found construction of a time-like form-factor for the axial isoscalar coupling similar to Eq. 5.20, so we develop one here. We will use two lucky coincidences. First, in the space-like region the lattice results which provide dipole fits for the isovector [47] and isoscalar [48] couplings find roughly the same "axial mass" (fitting parameter to the dipole form-factor which determines its Q 2 dependence) within uncertainties (m A = 1.169 (72) (27) GeV vs m A = 1.261 ± 0.188 GeV). The second coincidence is that the relevant meson for axial isoscalar coupling with I G (J P C ) = 0 − (0 +− ) is the h 1 with a mass m h 1 = 1.166 GeV and Γ h 1 = 0.375 GeV. Luckily these numbers are essentially the same as for the a 1 and the isovector case! The remaining parameters are α and δ. δ is associated with the intrinsic form-factor of the proton (i.e., the valence quark distribution), so we can again use δ = 0.397; α would need to be re-fit to space-like data. However, the lattice results above produce similar Q 2 dependence in the isoscalar and isovector form-factors, so we expect that the isoscalar α would be similar to the isovector one; therefore we use α ≈ 0.95.
To summarize, we argued that a reasonable first approximation for the time-like formfactors is to use the same momentum dependent form factor for both axial isovector and axial isoscalar couplings of the ALP to protons (the normalizations of these form factors are still different, and captured by the values of g A and g 0 , or, equivalently, by D, F and D s ). In any case, the time-like form factor is evaluated at k 2 = m 2 a and can be factored out from the overall rate. So if better models become available, projections or experimental constraints on the ALP coupling can be simply rescaled.

Other Channels
Production of heavier ALPs is no longer coherent over the proton or the target nucleus. Moreover, at larger masses some of the assumptions used in the bremsstrahlung calculation can become invalid, as does the use of chiral perturbation theory. We therefore estimate the ALP production cross-section in gluon fusion at larger ALP masses. Since the inclusive ALP production rate should be continuous as a function of mass this calculation also serves as a sanity check for the bremsstrahlung cross-section in mass range m a ∼ 1 − 2 GeV. The inclusive rate for gluon-gluon fusion can be written as where s ≈ 2m p E beam ( √ s ≈ 15 GeV for DarkQuest) and f g (x, Q 2 ) is the gluon PDF. We use the CT18NLO PDFs [77] via lhapdf [64] to evaluate the cross-section which is shown in Fig. 9. We have scaled the cross-section by A ≈ 56 to account for the number of nucleon targets per iron nucleus. This gluon fusion line terminates at m a ≈ 1.3 GeV since the PDFs cannot be evaluated at lower Q values. We see that our estimate of the bremsstrahlung cross-section clearly does not saturate the hadronic cross-section at m a 2 GeV. This is acceptable for our estimate of the DarkQuest sensitivity which does not extend much above a GeV. It might be possible to improve agreement at lower masses by including ALP mixing with heavier pseudoscalar resonances.

Sensitivity Projections and Existing Constraints
In this section we use the production mechanisms from Secs. 4 and 5 and the experimental set-up described in Sec. 2 to estimate the sensitivity of DarkQuest to photon-and gluoncoupled ALPs. In both scenarios we consider only the decays a → γγ; in the gluon-coupled case sensitivity can be improved further by considering hadronic decay modes as well. In Sec. 2.1 we argued that additional shielding is likely needed to suppress SM processes that can mimic the signal; following Ref. [16] we therefore require ALP decays to occur between z = 7 and 8 m after the front of FMAG. We will select events where both photons from a → γγ hit the 2 m×2 m ECAL placed at z = 19 m with E γ ≥ 1 GeV. The ECAL modules have transverse granularity of 5.5 cm, so we will also consider a requirement on minimum photon separation ≥ 5.5 cm [25]. This can be used as an additional handle for rejecting SM processes with more than 2 photons as briefly mentioned in Sec. 2.1. If DarkQuest is further enhanced with a preshower detector in front of the ECAL, an invariant mass measurement may be possible. Finally we will study the sensitivity of two phases of DarkQuest with N POT = 10 18 and 10 20 .
We only consider ALP production in the front of the FMAG. For photon-coupled ALPs this means the first interaction length for γ fusion and the production of secondary photons in the first interaction length followed by the Primakoff process within the first radiation length -see Sec. 4. For gluon-coupled ALPs this means we assume η ( ) production and pA bremsstrahlung occur in the first interaction length. This simplification enables us to consider a mono-energetic 120 GeV proton beam and results in a conservative estimate of the signal yield.

Photon coupling
The effect of the selections outlined above on event acceptance is shown in Fig. 12 for two benchmark masses and the two main photon-coupled ALP production mechanisms. We see that geometric requirements on the decay vertex and decay photons are the dominant sources of signal "loss". At small masses, the requirement of well-separated photons limits the sensitivity to short lifetimes since events more boosted ALPs also give more collimated photons. We combine these acceptances with the total cross sections shown in Fig. 5 to derive the sensitivity of DarkQuest to photon-coupled ALPs. Figures 13, 14 and 15 are the main results of this work for the photon-coupled ALP and show that DarkQuest can significantly improve on existing constraints, shown as shaded gray in these plots (these are briefly described in the following subsection).
As was pointed out in Ref. [16], ALP production via the Primakoff mechanism from secondary photons is the dominant process. However, despite having a smaller crosssection, γ * γ * fusion produces somewhat more boosted ALPs making this channel nearly as powerful in the short lifetime regime. This is evident in Fig. 13 where we show the sensitivity of each channel separately. In this and following figures we define the sensitivity contours as 10 signal events following Ref. [16].
Next we consider the impact of the photon separation selection in the ECAL. The sensitivity with and without this cut is shown in the left panel of Fig. 14. The requirement of observing well-separated photons is not stringent in most of the parameter space. It only starts to penalize the reach for m a 0.05 GeV. In the right panel of that figure we compare the sensitivity of the two proposed phases. Both phase 1 and phase 2 are sensitive to new parameter space. Phase 2 does not substantially improve the reach in the short lifetime/large coupling regime because of the exponentially falling acceptance in that direction. However, phase 2 gains almost a factor of 2 in the mass reach over phase 1.
Finally, we compare the sensitivity of DarkQuest phase 1 to other near-future accelerator efforts in Fig. 15. There we show projections for phase 1 of FASER [78] with 300 fb −1 ; NA62 with 10 18 POT [18]; NA64 with 5 × 10 12 EOT [54]; phase 0 of LUXE-NPOD [79]; Belle II with 20 fb −1 [80]; and a reanalysis of existing PrimEx data [50]. These experiments are either under construction, are already built or even have completed running, so results can be expected on a timescale comparable to the DarkQuest timeline of ∼ 5 years. We see that DarkQuest's unique baseline, intensity and beam energy make it extremely competitive and complimentary to other experiments. In fact, DarkQuest is able to probe masses and couplings that are not accessible at any other experiment already in phase 1. On a longer timescale, the photon-coupled ALP parameter space will be further probed by future phases of FASER [81] and LUXE-NPOD [79], and by Belle II [80] and DUNE [82]. Heavier photon-coupled ALPs can be produced and detected at LEP [83], ATLAS [84,85], CMS [86], 11 and in heavy ion [87] and electron-ion collisions [88].

Existing Constraints
The constraints on the photon-coupled ALPs are not subject to the same kind of theoretical consistency conditions as the gluon-coupled case, so we can make use of existing results. The bounds come from reanalyses of electron beam dump experiments E137 and E141 [80], proton beam dump experiments CHARM and νCAL [18], e + e − colliders LEP [89], BaBar [80] and Belle II [90], and from the photon-beam experiment PrimEx [50]. A breakdown of the various constraints is shown in Ref. [80]. We combined all of these bounds in the gray regions of Figs. 13, 14 and 15.

Gluon coupling
The effect of experimental selections outlined above on event acceptance is shown in Fig. 16 for m a = 0.5 GeV and two gluon-coupled ALP production mechanisms (we will see that  DarkQuest can cover new parameter space around this mass; lower masses are more constrained).
We see that geometric requirements on the decay vertex and decay photons are the dominant sources of acceptance loss. At this mass, the requirement of well-separated photons does not limit the sensitivity. It is interesting to note that ALPs arising from pA bremsstrahlung are very forward, so their acceptance is only dictated by the probability to decay within the fiducial volume: additional selections have no impact on the acceptance. We combine these acceptances with the total cross sections shown in Fig. 9 to derive the    [80] and a reanalysis of existing PrimEx data [50]. Shaded regions are existing constraints described in the text. sensitivity of DarkQuest to ALPs. Figure 17 is the main result of this work for the gluon-coupled ALP and shows that DarkQuest can improve on existing constraints shown as shaded gray (these are discussed in the following subsection). In particular DarkQuest can access new parameter space with m a 0.4 GeV, where it is competing against νCAL. This improvement is mainly enabled by DarkQuest's shorter baseline. The ALP branching fraction to photons becomes suppressed in this mass range (see Fig. 2); this means that sensitivity can be somewhat improved by considering other ALP decay channels, such as a → 3π and a → 2π + γ.
In the left panel of Fig. 17 we show the projected sensitivity for the two phases of DarkQuest. The futuristic phase 2 offers higher mass and lower coupling reach in the long-decay length regime (the lower part of the sensitivity contours).

Existing Constraints
Re-interpretation of existing searches is sensitive to the assumed model for ALP production and decay. This is particularly important for the gluon-coupled ALP because of the theoretical consistency issues relating to the cancellation of unphysical parameters. It is therefore beneficial to consider the entire set of existing searches within the same, consistent framework. While a comprehensive recasting of all gluon-coupled ALP constraints is beyond the scope of this work, we have performed simplified re-analyses of several experimental results, which provide the leading constraints in the parameter space relevant for DarkQuest. The complementarity of various experimental probes is shown in Fig. 18; we discuss these results below. For an exhaustive compilation of various bounds, see, e.g., Refs. [29,99].
We consider the following experimental results: • KOTO search for K L → π 0 +invisible [100]: This result was considered in the context of ALPs in Refs. [29,41]. We follow the analysis procedure of Ref. [41] and generate a sample of K L 's from the observed momentum distribution [101]. We then decay the K L 's into ALPs, and compute their probability to decay outside of the detector volume, while requiring the π 0 to be in the detector acceptance. We use these weights to construct an effective invisible branching fraction for K L → π 0 a, which is compared to the mass-dependent upper limit of Fig. 4 of Ref. [100].
• NA62 search for K + → π + + invisible [102] was also considered in Refs. [29,41]. In particular, our result for Γ(K + → π + a) matches that of Ref. [29]; however, we find that their excluded region does not quite match the limits reported by NA62 [102] (the experimental result features a gap in sensitivity for 0.1 GeV m a 0.16 GeV which is absent in the theory paper). Ref. [102] provides branching fraction limits as a function of lifetime, which we interpolate and translate into the ALP parameter space using the total ALP width.
• ATLAS Monojet search [104] provides an interpretation of their results in terms of a model where the ALP is produced via the gluon coupling but decays invisibly. As for the KOTO search, we estimate an effective invisible branching fraction from the kinematics of the produced ALPs. We simulate pp → a + j in MadGraph v3.1.1 [63] at √ s = 13 TeV, finding a median ALP energy of ∼ 700 GeV for the event selections of [104]. For each simulated event we find the probability of the ALP to decay outside of the detector which we take to have a radius of 12.5 meters (this coarse model of the detector acceptance should be improved in future studies). We then rescale the ATLAS bounds for invisibly-decaying ALPs, c GG /f 10 −3 GeV −1 , 13 by the effective invisible branching fraction.
• The old proton beam-dump searches νCAL [107,108] and CHARM [109] have been used to constrain a multitude of long-lived particles, including dark photons [110,111] and inelastic dark matter [19]. We simulated the production of gluon-coupled ALPs for both experiments in processes described in Sec. 5 and decayed them into photons. We follow the same procedure for estimating the acceptance and signal rate as for DarkQuest, taking into account the different experimental geometries and event selections. We find that νCAL has superior reach throughout the relevant parameter space, so we only show the νCAL result.
• In Ref. [112] the authors used experimental results from PIENU [113] and PIBETA [114] experiments to constrain rare pion decays π + → aeν, where the coupling of the ALP is described by a mixing angle with π 0 , sin θ. We recast their results into the c GG /f parameter space by using the following κ-independent proxy for the mixing angle: which we evaluated using the electroweak chiral Lagrangian in the limit of m η → ∞. The resulting constraints are labelled "πeν" and "πβ" in Fig. 18.
• We rescale the GlueX [115] constraint (based on the analysis of Ref. [50]) into our gluon coupling normalization (we do not recompute the limit in the κ-independent formalism).
We also translated the bounds on the photon coupled-ALPs into the c GG /f parameter space using the induced photon coupling in Eq. B.1 and accounting for the reduced branching fraction of a → γγ at larger m a . Only two photon-only results cover new parameter space compared to the hadronic experiments above. These are the Belle II [90] and LEP [89] searches for e + e − → γa(γγ).
The complementarity of various experimental probes is shown in Fig. 18. In addition to direct searches for gluon and the (induced) photon coupling, we highlight parameter space that may feature UV-dependent bounds coming from additional QCD-charged particles below the TeV scale (the parameter space above the lower dashed line) and the chromomagnetic dipole moment of the top quark µ t (the parameter space above the upper dashed line) [29]. The former line is calculated in the minimal KSVZ model [116,117] with a single pair of vector-like quarks with no non-QCD charges; this model generates c GG = 1 and the quark mass is bounded above by 4πf / √ 2. Note that it is not trivial to apply an existing collider search to derive a stringent constraint on f , since the exotic quarks are stable GlueX πβ πeν ATLAS Monojet coloured states µ t Figure 18. Constraints on the gluon-coupled ALPs evaluated in this work. The recasting procedure for each bound is described in Sec. 6.2.1. Note that we are considering the model where only the gluon coupling is present at low scales. If we defined our theory at a high scale, RG evolution would generate a multitude of other interactions leading to additional constraints [29,99]. The dashed lines indicate the presence of new QCD-charged matter at the TeV scale (lower dashed line) and a chromomagnetic dipole moment of the top quark (upper dashed line). due to a conserved U (1) [118,119]. Therefore any collider analysis must make specific assumptions about their decays. We leave such an exploration to future work. Similarly, the constraint from the chromomagnetic moment of the top quark is also sensitive to the direct coupling of the ALP to t, which we have set to zero. We chose not to shade the parameter space highlighted by the dashed lines; we combined the remaining bounds in the gray regions of Fig. 17.
As mentioned above, a consistent UV treatment, such as that adopted in Ref. [23] where the gluon coupling is defined at a high scale, generates a multitude of other couplings that result in additional bounds in the parameter space of Fig. 18 -see, e.g., Refs. [29,99].

Conclusion
The dimension-five ALP couplings are some of the leading candidates for BSM physics interacting with SM particles from the perspective of effective field theory. It is therefore important to continue testing these models with existing and future experiments. In this work we showed that the proposed proton beam-dump DarkQuest, will be able to probe new parameter space in models where the ALP couples dominantly to photons and gluons. We expect that similar results can be attained in more general scenarios, e.g., in which the ALP has interactions with quarks or leptons.
Our analysis contained several new aspects of ALP production and decays. In computing amplitudes involving ALPs in chiral perturbation theory we carefully tracked their dependence on unphysical parameters, ensuring their cancellation in the final result. We also provided new calculations of ALP coupling to nucleons in the three-flavour theory, and used it to estimate ALP emission rate in proton-nucleus bremsstrahlung. It would be interesting to compare this mechanism to other methods used to evaluate ALP production in hadronic interactions, such as emission in a parton shower.
We made several simplifying assumptions in estimating the sensitivity of proton beamdump experiments to ALPs, mainly relating to the production of mesons and photons in beam-target collisions. We took the interactions to happen in the first interaction or radiation length, neglecting the possible production deeper in the target with the attenuated beam. Clearly this yields a conservative estimate of the total yield, so it would be useful to study ALP production in more sophisticated simulations to see whether the true sensitivity is appreciably stronger. We also focused on ALP decays to photons; while in the photon-coupled case this is the only decay channel, gluon-coupled ALPs have a multitude of available channels in the parameter space to which DarkQuest is sensitive. These final states can also be looked for to recover the branching fraction penalty associated with a → γγ-only searches.
Finally, our compilation of bounds on the gluon-coupled ALP in Fig. 18 revealed several gaps in experimental coverage. It is interesting to think of existing or future observations that can be used to test these in a model-agnostic way. All of these occur at fairly large couplings so naively this should be straightforward. One obvious possibility are indirect bounds from searches for QCD-coupled states that are needed to generate the dimension-five ALP gluon coupling. However, in minimal models these particles carry a new baryon number, making them potentially collider-stable. As a result, any collider search is necessarily sensitive to the operator mediating their decay, which requires specifying yet more unknown UV physics. Clearly, it would be more satisfying to close these gaps with observables entirely within the ALP effective theory.

A Three-Flavour Mixing
Our goal is to derive the ALP-meson mixing in the three-flavour regime, so we must include π 0 , η and η . Let the chiral basis states be π 0 , η 8 and η 1 . The starting Lagrangian is then where we used ln det Σ = tr 2iΠ/f π . For now we assume that the full chiral symmetry is U (3) L × U (3) R ; this is softly broken by the mass term and by the anomaly as we will discuss below. The pion field that follows from this assumption is The coefficient in front of η 1 is chosen to ensure the correct normalization of its kinetic term; it also means that the corresponding generator of the U (1) A transformation has the same normalization as the non-Abelian ones (in the Peskin & Schroeder convention [124]): The anomaly contribution in Eq. A.1d is the leading term that breaks U (3) L × U (3) R to SU (3) L ×SU (3) R ×U (1) V and is consistent with large-N arguments [125]. Its coefficient is chosen such that if the mass term preserved the U (1) (i.e., it had vanishing trace), it sets the mass of the η 1 state. It is also important to note that the anomaly term above is the one after doing the field redefinition to eliminate the GG coupling of the ALP. If we did not perform this transformation the ALP would appear in the combination (A.4) The relative coefficient here can be fixed by demanding that the ALP is a spurion of the U (1) A transformation. The rotation that eliminates GG coupling is Σ → e 2iα Σ, where α = −κc GG a/f ; one can check that the same transformation eliminates the ALP from the anomaly term.
To leading order in a/f , the ALP-dependent quark mass matrix above is The kinetic mixing termĉ qq receives contributions both from "fundamental" ALP-quark couplings and terms ∝ κ; since we are working in the simplified case in which the ALP only couples to gluons, only the latter terms are present. We will not specialize to this case until very end though, we will only assume thatĉ qq is diagonal: Our goal is to bring the ALP kinetic term into a canonical form and then diagonalize the resulting ALP-meson mass matrix. We can write the ALP-pion Lagrangian as where ϕ = (π 0 , η 8 , η 0 , a) T and In the above we define variables that help to take the isospin limit: (A.10) The isospin limit corresponds to taking δ I → 0. We have also defined tree-level meson masses which agree with, e.g., Ref. [39].
Perturbatively diagonalizing and normalizing Z yields a modified mass matrix  where following Ref. [35], we have defined shorthands for the various elements in this matrix that make explicit dependence on small quantities δ I and = f π /f (the m 2 a entry also gets shifted, but only at O(f 2 π /f 2 ), so we dropped that term). The explicit expressions for the matrix entries are not very illuminating. The remaining meson mixings η 8 − η 1 , π 0 − η, π 0 − η are removed by perturbative pairwise diagonalization. For the η 8 − η 1 system we find the mixing angle tan 2θ ηη = 2M 2 (A. 13) In practice this angle is experimentally determined from the ratio of η and η partial widths into photons (see the Quark Model review in Ref. [38]); the approximation of Ref. [35] uses instead sin θ ηη = −1/3 for simplicity, which we use throughout this work.
The mass matrix after this transformation can be written in terms of new shorthand variables:  (A.14) Next we remove the π 0 − η and π 0 − η mixings. The matrix entries that mix these states are proportional to m 2 π and the isospin breaking parameter δ I , so the mixing angles will scale as m 2 π δ I /m 2 η ∼ 0.02 (0.007) for η and η respectively; we will therefore solve for these mixing angles perturbatively and drop higher order terms in these quantities. We find sin θ π 0 η ≈ M 2 π 0 η δ I where the diagonal entries now represent the physical masses of the mesons (the mixing with the ALP will change this only by O(f π /f )).
We can now finally diagonalize the ALP-meson mixing. As before, we do this by a series of three pair-wise rotations to eliminate the mixing of the ALP with each of the mesons. Solving for the mixing angles to remove off-diagonal terms in the mass matrix we find: sin θ π 0 a ≈ − The terms highlighted in blue here are proportional to δ I and therefore represent the leading isospin-breaking corrections to the mixing angles; note that these are not suppressed by the small ratio m 2 π /m 2 η .

B ALP Decays Into Photons in Three-Flavour χPT
The three-flavour equivalent of the expression Eq. 3.14 is where the coefficients in front of the mixing angles in the first line are ratios of the electromagnetic anomaly for each meson relative to the π 0 . In the second line made use of the explicit expressions forĉ γγ andĉ qq in terms of c GG and κ q (Eq. 3.5), the expressions for the mixing angles in terms of tree-level meson masses, and the meson masses in terms of the quark masses. The third line arises from the leading isospin-violating contributions that appear in the mixing angles in Eqs. A.18c (other isospin violating contributions at the same order in δ I are suppressed by m 2 π /m 2 η ( ) ). Note that we have not expressed the η 8 and η 1 masses in terms of the physical η and η mass since this facilitates taking interesting limits below.
This expression is remarkable in that all of the unphysical κ q dependence has cancelled; the result (in the δ I → 0 limit) is proportional to m 0 , the coefficient in front of the anomaly term in the chiral Lagrangian, Eq. A.1, In addition to the vanishing of the κ q dependence we can perform another consistency check by taking the two-flavour limit to see if we reproduce both terms of Eq. 3.14. In order to do this, we can express m η 1 , m η 8 and tan θ ηη in terms for the quark masses and m 0 ; one can then take the limits m s → ∞ and m 0 → ∞. As long as the mixing angle is self-consistently chosen (i.e., via Eq. A.13), the order of the limits does not matter. In the limit m s , m 0 → ∞, Eq. B.1 gives matching exactly Eq. 3.14.
C Rare Eta Decays η and η decays can be computed from the Lagrangian in Eq. A.1 by using the expressions of Eq. A. 19. For η ( ) → π 0 π 0 a decays the only terms that contribute come from the mass terms; for the decays involving charged pions in the final state the kinetic mixing terms also contribute. For generic η − η mixing angles, the expressions for the amplitudes corresponding to these processes are long, so we will present the limit sin θ ηη = −1/3 which significantly simplifies the expressions; this limit was also used in Refs. [35,126]. The measured mixing angle differs from this by a few percent (see the discussion in Ref. [39]), but this is sufficient accuracy for BSM calculations. The other approximation is to drop sub-leading isospin violating terms ∼ δ I m 2 π /m 2 η -these are terms that are proportional to sin θ πη ( ) that are not multiplied by m 2 η .
Under these assumptions we find the following expressions before plugging in the explicit expressions for the mixing matrix elements (Eq. A.20): At this point it is clear that there must be some non-trivial cancellations among the κ q between the various terms. Plugging in the mixing elements in Eq. A.20 we find at leading order in isospin violating parameter δ I A(η → π 0 π 0 a) = An important observation is that all κ q dependence has disappeared; in the π + π − final states this requires the combination of contributions both from the kinetic and the mass terms of the chiral Lagrangian. Moreover, this cancellation is not specific to the sin θ ηη = −1/3 choice; we have checked that it is true for an arbitrary mixing angle. The isospin breaking contributions arise at O(δ 2 I ) because πa ∝ δ I is also multiplied by δ I in the amplitudes above.

D Rare Kaon Decays
For brevity we present here the amplitudes in the limit of decoupled singlet meson and to leading order in isospin violation, which matches the assumptions in, e.g., Ref. [23]. In our numerics we use the full results that include η − η mixing. The charged kaon decay amplitudes are The corresponding K + amplitude has the opposite sign. The δ I → 0 limit of this expression matches the result of Ref. [23] once the difference in f π conventions is taken into account.
For the rare neutral kaon decays we first compute the amplitude for the neutral stronginteraction eigenstates K 0 and K 0 : (D.2b) The corresponding K 0 amplitude has the opposite sign. We find that the leading isospinviolating contribution is important at small ALP masses, and changes the amplitude by up to ∼ 50%. The amplitudes of the physical weak eigenstates K S,L follow from the definitions K S = K 1 + K K 2 , K L = K 2 − K K 1 , where K 1,2 = (K 0 − K 0 )/ √ 2 are the CP eigenstates and K ≈ 2.23 × 10 −3 is the CP violation parameter:

E Proton Bremsstrahlung
We follow the procedure outlined in Ref. [67] and consider the process p(p) + p(p t ) → a(k) + f (p f ) where the four-momenta labels are given in the parentheses; f represents an unspecified hadronic final state. The intermediate beam proton has momentum p . We