Prospects for detecting axionlike particles at the Coherent CAPTAIN-Mills experiment

We show results from the Coherent CAPTAIN Mills (CCM) 2019 engineering run which begin to constrain regions of parameter space for axion-like particles (ALPs) produced in electromagnetic particle showers in an 800 MeV proton beam dump, and further investigate the sensitivity of ongoing data-taking campaigns for the CCM200 upgraded detector. Based on beam-on background estimates from the engineering run, we make realistic extrapolations for background reduction based on expected shielding improvements, reduced beam width, and analysis-based techniques for background rejection. We obtain reach projections for two classes of signatures; ALPs coupled primarily to photons can be produced in the tungsten target via the Primakoff process, and then produce a gamma-ray signal in the Liquid Argon (LAr) CCM detector either via inverse Primakoff scattering or decay to a photon pair. ALPs with significant electron couplings have several additional production mechanisms (Compton scattering, $e^+e^-$ annihilation, ALP-bremsstrahlung) and detection modes (inverse Compton scattering, external $e^+e^-$ pair conversion, and decay to $e^+e^-$). In some regions, the constraint is marginally better than both astrophysical and terrestrial constraints. With the beginning of a three year run, CCM will be more sensitive to this parameter space by up to an order of magnitude for both ALP-photon and ALP-electron couplings. The CCM experiment will also have sensitivity to well-motivated parameter space of QCD axion models. It is only a recent realization that accelerator-based large volume liquid argon detectors designed for low energy coherent neutrino and dark matter scattering searches are also ideal for probing ALPs in the unexplored $\sim$MeV mass scale.

We show results from the Coherent CAPTAIN Mills (CCM) 2019 engineering run which begin to constrain regions of parameter space for axion-like particles (ALPs) produced in electromagnetic particle showers in an 800 MeV proton beam dump, and further investigate the sensitivity of ongoing data-taking campaigns for the CCM200 upgraded detector. Based on beam-on background estimates from the engineering run, we make realistic extrapolations for background reduction based on expected shielding improvements, reduced beam width, and analysis-based techniques for background rejection. We obtain reach projections for two classes of signatures; ALPs coupled primarily to photons can be produced in the tungsten target via the Primakoff process, and then produce a gamma-ray signal in the Liquid Argon (LAr) CCM detector either via inverse Primakoff scattering or decay to a photon pair. ALPs with significant electron couplings have several additional production mechanisms (Compton scattering, e + e − annihilation, ALP-bremsstrahlung) and detection modes (inverse Compton scattering, external e + e − pair conversion, and decay to e + e − ). In some regions, the constraint is marginally better than both astrophysical and terrestrial constraints. With the beginning of a three year run, CCM will be more sensitive to this parameter space by up to an order of magnitude for both ALP-photon and ALP-electron couplings. The CCM experiment will also have sensitivity to well-motivated parameter space of QCD axion models. It is only a recent realization that accelerator-based large volume liquid argon detectors designed for low energy coherent neutrino and dark matter scattering searches are also ideal for probing ALPs in the unexplored ∼MeV mass scale.
In the 70s and 80s, the QCD axion was extensively arXiv:2112.09979v5 [hep-ph] 26 May 2023 searched for in beam dump, fixed target, and reactor experiments (see, for example, refs. [32,[62][63][64][65][66][67][68][69]). This effort has been revived with the modern Intensity Frontier experimental program to explore more generic ALP signals. Over the past three decades advances in accelerators have enabled modern experiments to gain an order of magnitude or more in instantaneous beam luminosity, while advances in instrumentation have led to significant improvements in detection efficiency, in energy, spatial, and timing resolution, lower detection thresholds, and particle identification (see e.g., refs. [70][71][72][73] for a review). In particular, the high intensity photon flux and associated electromagnetic cascades in reactor and accelerator neutrino experiments offer new opportunities to explore ALP production via its electromagnetic and leptonic couplings [30,31,74].
For ALPs that couple predominantly to photons, the Primakoff and inverse-Primakoff processes can be exploited for ALP production and detection, respectively. These processes, illustrated in Figs. 1a, 1b, are coherently enhanced by a factor of Z 2 , where Z is the atomic number of the target nucleus. The ALP flux can produce electromagnetic signals in the detector via inverse Primakoff scattering or decay to a photon pair within the detector's fiducial volume (Fig. 1d).
For ALPs with significant couplings to electrons, other processes can also contribute to its production, including Compton-like scattering (Fig. 1c), e + e − annihilation (Figs. 1e, 1g), and ALP-bremsstrahlung (Fig. 1i). Such ALPs can be detected via inverse-Compton scattering (Fig. 1f) and e + e − conversion (Fig. 1j), or, if sufficiently short-lived, via decay to e + e − within the detector (Fig. 1h). Detailed descriptions of the amplitudes relevant for these processes are given in Appendix A.
Previous studies have shown that ongoing reactorbased neutrino experiments such as CONUS, CONNIE, MINER etc. [30,31], and upcoming accelerator-based neutrino experiments such as DUNE [74] will be able to probe parameter space for ALPs in the MeV mass range coupling to electrons, photons, and nucleons. This mass range has remained inaccessible to terrestrial and astrophysical observations to date. In this paper, we investigate the sensitivity of the Coherent CAPTAIN-Mills (CCM) experiment to ALPs coupled electromagnetically or electronically and push the sensitivity envelope in the MeV mass region of ALP parameter space. The LANSCE accelerator and proton storage ring provides an 800 MeV, 100 µA, short 290 ns pulse of protons (triangular shape) impinging on a thick tungsten target and producing significant hadronic activity and a high intensity flux of photons and electromagnetic cascades in the O(0.1 − 1000) MeV energy range. ALPs produced through photons and e ± interacting with the tungsten target material could be detected at the 5-ton (fiducial cylinder approximately 1 m height by 2 m diameter) liquid Argon detector located 23 meters away from the target and 90 • from the beam direction, as shown in Fig. 2. ALP's that interact electromagnetically in the CCM liquid Argon will produce a copious shower of scintillation light (∼40,000 photons/MeV) at 128 nm with a prompt 6 nsec and slower 1.6 µsec decay constant [75,76]. The scintillation light is wavelength shifted by Tetraphenyl Butadiene (TPB) surfaces and then detected by an array of PMT's. See [77] for detection details.
In 2019 a six week engineering beam run was performed with the CCM120 detector, named due to it having 120 inward pointing main PMTs. The CCM120 experiment met expectations and performed a sensitive search for sub-GeV dark matter via coherent nuclear scattering with 1.79 × 10 21 Protons On Target (POT) [77,78]. Due to the intense scintillation light production and short 14 cm radiation length in LAr [76], the relatively large CCM detector has good response to electromagnetic signal events in the energy range from 100 keV up to 10's of MeV. This low energy kinematic range, which is sensitive to ∼MeV ALP mass range, has not been previously explored at proton beam dump experiments. CCM's novel sensitivity to this region could probe new parameter space of BSM particle production and yield new insights into the nature of the LSND [79] and MiniBooNE event excesses [80,81]. Another key feature of CCM is that it uses fast beam and detector timing to isolate prompt ultra-relativistic particles, which ALPs in the MeV mass range may be for the energy scale of the Lujan proton beam source. This can distinguish them from the significantly slower neutron backgrounds that arrive approximately 225 ns after the start of the beam pulse (relativistic particles traverse the 23 m distance in 76.6 ns) [77]. Furthermore, the Lujan beam low duty factor of ∼ 10 −5 and extensive shielding are efficient at rejecting steady state backgrounds from cosmic rays, neutron activation, and internal radioactivity from PMTs and 39 Ar. In order to determine the sensitivity reach of CCM's ongoing run, we use the beam-on background distribution determined from the recent CCM120 run [77], with a further expected factor of 100 reduction from extensive improvements in shielding, veto rejection, energy and spatial resolution, particle identification analysis, and reduced beam width.
In § II we discuss the ALP phenomenological models, and their production and detection modes at CCM; in § III we review the treatment of backgrounds and cuts to optimize the signal efficiency at CCM120; in § IV we set limits on the ALP parameter space from current CCM120 data and projected sensitivities for the ongoing 3-year physics run, and in § V we conclude. Additional details on the signal prediction are outlined in Appendix A and details on the optical model and reconstruction in Appendix B.

II. ALP MODELS AND THEIR ENERGY SPECTR A AT CCM
The focus of this analysis is on generic models of ALPs with couplings to photons and electrons. These interactions are parameterized by the following Lagrangian terms: where e is the electron Dirac fermion, and (F µν ) F µν is the electromagnetic (dual-)field strength tensor. For specific ALP models, including the QCD axion and its variants, the couplings g aγ and g ae are not independent. For the purposes of this study, however, we will adopt a simplified model approach by considering two limiting cases: in the first, we set g ae = 0, so that the ALP phenomenology is completely determined by its electromagnetic interactions parameterized by g aγ ; and in the second case, we assume that g ae is sufficiently large to dominate ALP production and detection at CCM-this limit holds when g ae ≫ α m e g aγ . At LANSCE's Lujan source, the 800 MeV proton beam impinging on the tungsten target produces a high intensity photon flux from cascades, neutron capture, pion decays, etc. These processes were modelled using GEANT4 10.7 with the QGSP_BIC_HP library [82], and the resulting photon, electron, and positron spectra are shown in Fig. 3 for the energy range E γ,e ± = 0.1 − 500 MeV. Flux errors are primarily associated with differential neutron production uncertainties and for tungsten estimated to be less than 10% [83]. ALPs produced from the processes shown in Fig. 1 could then propagate to the CCM detector, where they could produce an electromagnetic signal via inverse Primakoff or Comptonlike scattering, diphoton decay, or e + e − conversion or decay. We show in the next section that the background . ALP-induced differential event rates at the CCM detector. Top panel: energy spectrum of 1-γ and 2-γ final states from inverse Primakoff scattering and diphoton decays. Bottom panel: energy spectrum of (e − , γ) final states from inverse Compton scattering and e + e − conversion/decays. In some scenarios the signal abruptly turns on, making it clearly identifiable relative to the backgrounds that generally fall off at higher energy (as described in §. III). spectrum relevant to these signal channels falls off exponentially at energies greater than a few MeV; thus, the ALP signals could be potentially visible due to the harder spectral shapes of γ's and e ± 's from ALP scattering and decays. Monte Carlo simulations of the CCM detector response to gammas and electrons from ALP events indicate that the signal reconstruction efficiency for 5 fiducial tons of LAr is above 50% for events above 100 keV.
In Fig. 4 (top panel) we show the energy spectra of ALP-induced events with one or two photons in the final state resulting from either inverse Primakoff scattering or diphoton decays. Signals from heavier ALPs are characterized by a significantly harder spectrum, which in principle makes them more eas-ily distinguishable from backgrounds. On the other hand, the ALP lifetime decreases with ALP mass at fixed coupling, and therefore heavier ALPs need to be more boosted in order to decay within the detector fiducial volume.
In Fig. 4 (bottom panel) we show the energy spectra of ALP-induced events with an electron or an e + e − pair in the final state, arbitrarily normalized to yield a non-negligible signal-to-background ratio. At lower ALP masses m a ≲ 2 m e the signal is dominated by inverse Compton scattering, whereas for ALP masses above the e + e − threshold, the signal is dominated by a → e + e − decays exhibiting a "lineshape" spectral feature. This is because the dominant ALP production mechanism in this kinematic range is resonant e + e − annihilation to an ALP of energy E a = m 2 a /(2m e ). If the ALP lifetime is such that it decays outside the detector, the ALP signal is instead dominated by external ALP conversion to an e + e − pair in the detector.

III. ALP SEARCHES AT CCM120
FIG. 5. The event rate after the length cut in CCM120 around the measured beam T0 (defined as T=0ns) when the earliest speed of light particles from the beam would arrive in CCM (see [77]). For short nucleon scattering like events (<44ns) the allowed ROI end time is 190ns after T0. For long electromagnetic like events, changes in length cut efficiency began 170ns after T0. For some energy bins efficiency changes began 150 ns after T0, which defined the final analysis ROI of 150 ns.
In 2019 CCM120 carried out an initial engineering run to test the detector's capability to search for FIG. 7. The ALP simulation (see Appendix B) vs. prebeam radius showing the preference for a tighter radius cut. ALP events followed an isotropic pattern throughout the detector, while background events from radioactivity preferentially occurred near the edges. The position resolution is ±10 cm near the center, measured from calibration data.
BSM physics, including ALPs [84]. This search was composed of three steps: defining a beam-related neutron background free region of interest (ROI) in which to search for near speed of light ALPs, determining the signal characteristics of ALP events through analysis cuts to obtain the data, and defining confidence limits from the observed data excess.
The first step was to define a beam-related back- ground free ROI. While CCM is able to use the 10 µs of prebeam data collected to make highly accurate background rate predictions for the steady-state backgrounds such as cosmic rays, radioactivity, long lived neutron activation, etc, such prediction accuracy is not possible for the beam-related fast neutron backgrounds. As such a region free of such backgrounds is desirable. This region is defined as the time region starting with the earliest arrival of speed of light particles from the target created by the incident proton beam T 0 and ending with the so-called neutron beam-turn on. T 0 was determined through the use of direct measurement of speed of light particles from the gamma flash generated in the target with an EJ-301 scintillation detector on a dedicated beam line directly viewing the target. This procedure is described in more detail in ref. [77] and [84]. The end time was determined from changes in the efficiency of the selection cuts and the slope of the rate of high energy signal-like events, as shown in Fig. 5. The beam related neutron free ROI was determined from the data to be 150 ns. This is consistent with MCNP simulations for the earliest fast neutrons (kinetic energy > 100 MeV) arriving through the bulk shielding (5 m steel and 2 m concrete) and neutron time of flight measurements with multiple EJ-301 detectors in the vicinity of CCM behind the bulk shielding. This contribution is determined to be negligible.
The second step was the definition of the event selection cuts. A number of data quality cuts were used in all CCM120 searches, specifically a beam quality cut to ensure good quality beam triggers, a previous event cut to remove triplet events, and a veto cut to remove events from outside the detector. These cuts are described in more detail in ref. [77].
The signal selection cuts were added after these quality cuts and were defined based on calibration data using a 22 Na γ source (Fig. 6) and detector simulations (See Appendix B) of various energy electromagnetic events in CCM. The first such cut was a length cut requiring events >38 ns. Using a dynamic length event builder which defined conditions for the end of the event, the event length worked as a crude particle identification (PID) method. This PID was able to distinguish electromagnetic (electron or photon) events from nuclear recoil (neutrons or coherent scattering) events using the triplet light production of the former. This cut was >99% efficient on ALP signals of all energies, and 66% efficient on prebeam background.
The second cut was a minimum energy cut from the high energy (>1 MeV) limits of the ALP search. The cut required events to have more than 10 PE of visible energy, approximately equivalent to 1 MeV at the center of the detector. This cut was ∼ 95% efficient on ALP signals above 1 MeV from simulations while being 30% efficient on prebeam background, as the majority of the background was <1 MeV. The third cut was a strict position cut requiring radius <80 cm and height between -40 and 40 cm from the detector center (Fig. 7). This cut was ∼ 75% efficient on signal and 34% efficient on background.
Finally, two shape cuts were implemented requiring that the proportion of TPB coated PMTs seeing light was >0.6 and that the single PMT which saw the most light saw less than 20% of the total charge (Fig 8). These two cuts combined were ∼ 90% efficient on signal and 40% efficient on background. Over all reconstruction cuts, CCM120's selection criteria was ∼ 65% efficient on signal with variation according to the energy distribution of the signal (see Table 1), and flat 3.2% efficient on background.
The resulting data and background spectrum are shown in Fig. 9. The number of data events in the signal region is 294590, with a scaled prebeam steady state background prediction of 294614.3 ± 241.7 (syst) ±542.8 (stat) and a subtraction of −24.3 ± 594.2. The error is statistics dominated. The systematic error on the steady state background prediction accounts for the additional variance from the systematic decay in background rates observed on the microsecond scale that are characteristic of radioactive decay from beam thermal neutron activation. After performing an exponential fit to account for this background, the error on the fit prediction is used as the systematic uncertainty on the background prediction. There is no significant excess over the entire energy range, however, axion signals can have strong energy dependence as shown in Fig. 4 and require axion model fits to properly analyze. Since the χ 2 of the data points in Fig. 10 are consistent with no excess at the 2σ level, there is no significant axion signal in the 2019 data, but the data can be used to determine exclusion regions as a function of model parameters. The third step in the CCM120 search was to use the data so selected to define a set of confidence limits. We defined these confidence limits as exclusion regions using a ∆χ 2 analysis and the asymptotic approximation given by Wilks' theorem. Using the ALP induced event rates (Fig. 4) for each ALP model, including mass and coupling, we defined a χ 2 value for each point on the phase space map in the two ALP models, where the χ 2 compared background (bkgnd i ) and ALP prediction (pred i ) to the data (data i ) with respect to the total error σ i (stat+syst) across all reconstructed energy bins i as seen in Eq. (2).
This χ 2 was used to evaluate the goodness of fit of the models shown in Fig. 10. We then defined ∆χ 2 model = χ 2 model − χ 2 best−f it and from Wilk's theorem used the resultant ∆χ 2 values as a test statistic for a χ 2 test with 2 degrees of freedom across the entire ALP parameter space. We could thus exclude phase space with ∆χ 2 > 4.61 to a 90% confidence limit. We also calculated the experimental sensitivity using only background predictions, with an otherwise similar ∆χ 2 search methodology. The confidence limits and sensitivities so generated are shown in Figs. 11 and 12.
For g aγ coupling, the existing beam-dump limits are better than the CCM120 exclusion limit as shown in Fig. 11. For g ae coupling (as shown in FIG. 10. The subtraction reconstructed energy spectrum between data and background prediction for the CCM120 ALP search, compared to predicted event spectra at various masses for a photon (top) or electron (bottom) coupling (signal parameters chosen to show large event rates). Smearing uncertainties at ±1σ are shown by the shaded bands. See Appendix B for details on the simulation and reconstruction. Fig.12), however, the CCM120 data constrains regions of parameter space in the m a ≲ MeV mass range, which so far has only been probed by astrophysical observations. In some regions of the parameter space, the CCM constraint is found to be marginally better than both astrophysical and terrestrial constraints.

IV. REACH PROJECTIONS IN ALP PARAMETER SPACE
The new and upgraded CCM200 detector (with 200 inward pointing PMTs, improved TPB foils, double veto PMTs, and LAr filtration (see ref. [77]) was completed in 2021 and has begun a 3-year  FIG. 11. The expected and actual 90% CLs from CCM120 for the ALP-photon coupling gaγ. Also included is the projection region for CCM200 three year run using background taken from CCM120's spectrum reduced by two orders of magnitude for various conservative improvements (dashed green line) and a background free assumption (extent of shaded green region). QCD axion model parameter space for the KSVZ benchmark scenario spans the region indicated by the arrows [85].
physics run from 2023 to 2025, which expects to collect an integrated luminosity of 2.25 × 10 22 POT. In this section we present projections for CCM200's sensitivity to the two ALP models discussed. These projections were obtained under two background scenarios: (i) a conservative one assuming background levels shown in Fig. 9 with an extra flat factor of 100 suppression for various improvements, and (ii) the background-free hypothesis. The measured beam-on background distribution at CCM120 as a function of visible energy has been studied in [77], and shown in Fig. 9 for the ALP cuts (from ref. [84]). We conservatively assume that, with the CCM200 improvements in shielding, particle identification, and energy and spatial resolution, an additional factor of ∼100 in background rejection is possible. To first order this reduction is assumed to be flat in energy and has been demonstrated by data studies and simulations. The ultimate goal is to go beyond this conservative estimate towards zero background with further improvements in added shielding, improved energy resolution with LAr filtration, reduced beam width by a factor of two or more, and analyses improvements in background rejection. The 90% C. L. projection reach of CCM200 in the parameter space of ALP coupling vs mass was obtained through simulations of pseudo-experiments under the background-only hypothesis and a standard Pearson ∆χ 2 test statistic. Fig. 11 shows CCM200 sensitivity to an ALP model with electromagnetic couplings in the parameter space of ALP g aγ coupling vs mass (see Eq. (1)). In the ALP mass range m a ≲ 10 keV, the constraints from NOMAD [86] and e + e − → γ +invisible states [87] are relevant. Here, the detection signal is dominated by inverse Primakoff scattering, and therefore CCM's reach is approximately mass-independent. At higher masses, where the signal is dominated by diphoton ALP decays, CCM200 will offer an improvement in reach over previous beam dump constraints, which is compelling since it will cover an unexplored region of ALP masses and couplings overlapping the KSVZ QCD axion parameter space [85,88,89], including the so-called "cosmological triangle" [90,91] formed by the boundary between astrophysical constraints [92][93][94][95][96][97][98][99], shown in gray, and the bounds from beam dumps [100,101]. Fig. 12 (top panel) shows CCM's reach to an ALP model with dominant electronic couplings in the parameter space of ALP g ae coupling vs mass (see Eq. (1)). Here, too, CCM200 will provide the strongest terrestrial constraints in the ALP mass The expected and actual 90% CLs from CCM120 for the ALP-electron coupling gae at tree-level (top) and with an effective photon coupling at loop level (bottom). Also included are projections for CCM200, using the same background specifications as in Fig. 11, for a 3-year run. QCD axion model parameter spaces for the DFSZ(I) and DFSZ(II) benchmark scenarios span the regions indicated by the arrows [85]. The region excluded by missing energy searches at NA64 is shown in gray, and the bound derived this work from the CCM120 engineering run is set at marginally lower couplings than the NA64 region. Even in the conservative assumption that loop-level a → γγ decays are not suppressed (bottom panel), CCM200 is projected to reach beyond the more stringent constraints set from E137 in this scenario.
range m a ≤ 2 m e , and will be able to probe parts of the QCD axion parameter space of DFSZ-I models 1 [103][104][105]. For masses within the range of 2 m e ≲ m a ≲ 10 MeV, CCM will be able to extend the reach to smaller g ae couplings beyond the parameter space excluded by E137 [106,107], including a QCD axion corner belonging to DFSZ-II models. In this region, the only competing constraints come from Supernova SN1987a bounds [108]. Other constraints shown in Fig. 12 were taken from [67] (Orsay), [68] (E141), [69] (E774), [109][110][111] (NA64), and [112] (stellar cooling). In Fig. 12 (bottom panel) we show the case in which an effective ALP-photon coupling is generated through an electron loop, permitting a → γγ decays. This modifies the parameter space such that E137 constrains lower mass parameter space below m a < 2m e from sensitivity to the γγ decay channel [113]. In this more conservative scenario, where loop-level a → γγ decays are not suppressed, CCM200 is projected to reach beyond the more stringent constraints set from E137 in this scenario.

V. CONCLUSION
This work illustrates the breadth of physics signatures that the CCM experiment will be able to explore. Being a proton beam experiment, CCM is uniquely suited to probe new physics signals initiated from hadronic processes. However, the high intensity of photons and electromagnetic cascades generated at the source will also allow CCM to probe new light particles that are coupled electromagnetically or electronically. A specially well-motivated class of models in this category are axion-like particles, which encompass the QCD axion that solves the strong CP problem, but also more generic light pseudoscalars from a dark sector.
In particular, given its unprecedented sensitivity to a lower energy kinematic range, CCM will be able to probe a variety of softer dark sector signatures ranging from coherent elastic dark matter/neutrino nuclear scattering, to O(0.1 − 100) MeV electromagnetic signals. This enables CCM to test a range of lighter and more weakly-coupled ALPs that was outside the reach of previous generations of beam dump/fixed target experiments. For the specific ALP models considered in this study, CCM can provide the leading terrestrial constraints in the sub-MeV ALP mass range. In this range, competing constraints stem mostly from bounds on stellar cooling, which are more model dependent and subject to large astrophysical uncertainties. Intriguingly, CCM can also probe surviving corners of the QCD axion parameter space in the m a ∼ O(0.1 − 1) MeV range that have proven difficult to exclude with past experiments and astrophysical observations.
In this paper, we analyzed the 2019 engineering run data for both g aγ and g ae couplings. We found that the CCM data has constrained regions of parameter space in the ≲ MeV mass range for the g ae coupling which so far has only been probed by astrophysical observations and in some regions of the parameter space, the CCM constraint is found to be marginally better than both astrophysical and terrestrial constraints. We also showed the projected sensitivity for a three year run and utilizing a background reduction of greater than two orders-ofmagnitude relative to the baseline measured during CCM's 2019 engineering run. With CCM200 ongoing detector performance improvements we expect to achieve these background levels starting in 2022 and with a 3-year physics run from 2023 to 2025 (collecting 2.25 × 10 22 POT) will achieve the sensitivities stated in this work. In the longer term, future planned upgrades that shorten the beam pulse width by an order of magnitude could provide even further background rejection and increased signal sensitivity. In the event of a discovery or observation of an excess, CCM's ability to detect visible energy (instead of missing energy/momentum) would also offer an advantage in narrowing down the possible BSM explanations for a putative signal.

ACKNOWLEDGEMENTS
We acknowledge support from the Department of Energy Office of Science, the National Science Foundation, Los Alamos National Laboratory LDRD program, the National Laboratories Office at Texas A&M University, and PAPIIT-UNAM grant No. IT100420. Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing. We also wish to acknowledge support from the LANSCE Lujan Center and LANL's Accelerator Operations and Technology (AOT) division. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001.  ALPs coupling to electrons via the pseudoscalar Yukawa operator, have many production and detection channels available to them that resemble similar scattering processes of the SM photon. The production and detection channels that are relevant for stopped-pion beam target experiments with this coupling are itemized below: ⋄ Compton scattering: The Compton-like scattering process for ALP production from photons scattering on electrons at rest, γ + e − → a + e − , is shown in Fig. 1c. It has a differential scattering cross-section, in light-cone coordinates, is the energy fraction distributed to the ALP in light-cone coordinates. This formula was presented in ref. [114], although we have corrected a sign mistake for the second term in square brackets in Eq. A2, also pointed out in ref. [31]. The differential cross section with respect to the outgoing ALP energy can be obtained by taking dσ/dE a = (1/E γ )dσ/dx, and the relevant bounds of integration over x can be taken from ref. [31] and references therein.

⋄ Inverse Compton Scattering:
The inverse-Compton scattering process (a + e − → γ + e − ) shown for ALP detection, shown in Fig. 1f, produces visible energy in the final state electron recoil and outgoing photon. The total cross section is [66,115]; σ(E a ) = g 2 ae α 8m e p a 2m 2 e (E a + m e )y (y + m 2 e ) 2 + 4m e (m 4 a + 2m 2 e m 2 a − 4m 2 e E 2 a ) y(y + m 2 e ) + 4m 2 e p 2 a + m 4 a p a y ln m e + E a + p a m e + E a − p a (A3) where y = m 2 a + 2m e E a .
⋄ Associated Production: The associated production of an ALP with a γ via electron-positron annihilation is shown in Figure 1e, with momenta e + (p + )e − (p − ) → a(k) + γ(q). The matrix element is computed below; The differential cross section is then given by the general expression for 2 → 2 scattering, which gives, in the CM frame (denoted with starred variables) Weighted samples of ALP 4-vectors can be drawn from this distribution, for example, by randomly drawing angles on the 2-sphere for the ALP in the CM frame, then transforming the resulting 4-vectors to the lab frame via a Lorentz boost. The MC weights are then given by dσ which is frame invariant in the limit of large sample size.
⋄ Axion Bremsstrahlung Production: ALP bremsstrahlung (e ± + N → e ± + N + a) shown in Fig. 1i via the interaction of electrons or positrons with the strong nuclear electric field of atoms in material was studied by Tsai [116]. The differential cross section as a function of outgoing ALP energy E a in the Weizsacker-Williams approximation is given as

⋄ Resonant Production:
The resonantly produced ALP from positrons annihilating on electrons at rest (e + e − → a, Fig. 1g) has a final state energy of This by itself is an interesting result, as it implies that the induced a flux from electron-positron annihilation will peak strongly close to ∼ m 2 a (MeV). The cross section in the electron rest frame with the narrow width approximation is where s = m 2 e + 2E + m e .
⋄ ALP external pair production: The process a(k)N (p 1 ) → N (p 2 )e + (ℓ + )e − (ℓ − ) shown in Fig. 1j, analogous to photon-lepton pair production in the SM, was computed in refs. [117,118] using the formalism and atomic form factors presented in ref. [119]. The cross section may be expressed in terms of invariants s = (k + p 1 ) 2 , s ℓℓ = (ℓ + + ℓ − ) 2 , and t = (p 1 − p 2 ) 2 and integrated in the dilepton CM frame; Here we define the leptonic and hadronic tensors; where F A (t) is the atomic form factor described in ref. [119].
⋄ Decays to electron-positron pairs: Decays to e + e − pairs are permitted for ALP masses m a > 2m e , i.e., m a ≳ 1 MeV. This decay is pictured in Fig. 1h, and the decay width is given by Eq. A15; (A15)

ALP-photon coupling Production and Detection Mechanisms
ALPs coupling to photons via the dimension-5 operator give rise to production and detection channels from Primakoff scattering as well as decays to two photons. The amplitudes for these processes are described below.
⋄ Primakoff and Inverse Primakoff Scattering: The matrix element for Primakoff scattering with a free, heavy fermion a(k)N (p) → γ(k ′ )N (p ′ ) with momentum transfer q = k − k ′ is depicted in Fig. 1a. The matrix element for scattering off a free fermion plane wave state is Squaring M free and evaluating the trace in terms of the usual mandelstam variables t = −q 2 and s = (k + p) 2 yields To apply this to a real atomic target, we can factorize the free matrix element and the atomic form factor separately to compute the interaction with the nuclear and electron cloud charge density in a straightforward way, as For a hydrogenic potential, we use Tsai's parameterization [116] of the atomic form factor F 2 (t), Z 2 a 4 q 4 (1 + a 2 q 2 ) 2 (A20) with a = 184.15e −1/2 Z −1/3 /m e . For inverse Primakoff scattering, the cross section is twice as big (since we no longer average over initial state photon polarizations) but is otherwise the same, with the appropriate re-definitions of Mandelstam invariant s = 2E a M + M 2 + m 2 a . If we instead use a simplified dipole form factor with screening length r 0 and express the differential cross section in the forward limit, the total cross section can be compactly expressed as [66,120,121] σ(k a ) = Z 2 αg 2 aγ 2 2r 2 0 k 2 a + 1 4r 2 0 k 2 a ln 1 + 4r 2 0 k 2 a − 1 , which can be also used for γ → a Primakoff scattering in the massless ALP limit after dividing by a factor of 2 to account for the polarization averaging.
⋄ Decays to photon pairs: The decay process a → γγ pictured in Fig. 1d is kinematically available at all ALP masses, although this decay is most relevant for masses in the MeV-GeV mass range at the energy and length scales of beam dump and stopped pion facilities. The width is

MCMC Method for Event Rate Distributions at CCM
We use an internal python library containing the matrix elements and cross sections listed above for input into a MCMC event generator. A schematic breakdown of the simulation pipeline is as follows: 1. Pass in fluxes of e ± , γ transport inside the tungsten target from GEANT4 simulation as arrays of energies (MeV) and weights (counts · POT −1 ).

2.
For each e + , e − , γ generate a weighted spectrum of ALPs using the appropriate production channel cross sections (for example, using the procedure outlined under Eq. A7).
3. Assume the ALP flux produced is isotropic, propagating the flux to the detector with a suppression factor of 1/(4πℓ 2 ) for ℓ = 23 m (CCM120).
4. Generate MCMC sampled distributions from the ALP fluxes propagated to the detector, appropriately reweighted to get the expected counts per exposure and visible energies in the ROI.
5. Convolve the distribution with the smearing matrix derived from the optical model and selection efficiencies described in Appendix B.
First, we review general formulae for predicting isotropic ALP fluxes produced inside the tungsten target. Photons produced and transported through the tungsten target are assumed to be entirely absorbed, and with minimal energy loss before getting stopped. In this case, the ALP production rate can be modeled by taking a ratio of the SM photon absorption cross section σ γ and the ALP production cross section σ a ; For ALPs produced from resonant production, bremsstrahlung, and associated production, the energy loss of the electrons and positrons in material must also be folded in. In the case of resonant production, this modifies the ALP production rate as where the track length probability function is the energy loss smearing function for the electron/positron radiation length t and target radiation thickness T [116]. The prefactor accounts for the radiation length X 0 in g/cm 2 , the atomic weight of tungsten A in g/mol, and Avogadro's number N A = 6.023 · 10 23 . For resonant production, the delta function in σ(E ′ ) removes the dE ′ integral and fixes the final state energy to the resonant energy. In the case of associated or bremsstrahlung production, we replace σ(E ′ ) with dσ/dE ′ and numerically carry out the integration over dtdE ′ dE ± . The propagation of these fluxes isotropically to a distance ℓ away from the production site will multiply the number event rate dN a /dt by a 1/(4πℓ 2 ) factor. In the limit that the distance ℓ from the source to the detector is much larger than the detector size, we can approximate the detector as a thin shell covering an area Ω d with thickness ∆ℓ. In this limit we treat ℓ as a constant over the transverse directions (θ, ϕ) over a sphere of length ℓ surrounding the source. The event rates for scattering and decays in terms of the ALP number rates produced in the W target, the detector particle number density n, and visible energy E vis can then be expressed simply as Scattering: The visible energy E vis is defined for scatterings (e.g. a + e − → γ + e − ) by the differential scattering cross section dσ/dE vis , and we have assumed for a → γγ and a → e + e − decays that the incoming ALP energy is totally converted inside the detector volume (e.g. E a = E vis ), before taking into account the smearing effects described in Appendix B. Here we have used the ALP lifetime in the laboratory, and the exponential law for decays, which together give us the probability that the ALP will survive until a distance ℓ from its production site, We then integrate the energy event rates given in Eq. A26 and Eq. A27 over the energy bins ROI and time exposure, convolving the true energy E vis with the smearing matrix described in Appendix B, to get the predicted signal events at CCM.
FIG. 13. The integral charge from laser pulses in each of the PMT rows, used for calculating the row ratios used for determining the OM. The difference in internal position and distance from the source between the rows was used to calibrate distance dependent OM properties.
light. However, due to the significant difference in photon energy and thus detector response between the laser UV -213 nm -and the argon scintillation peak at 128 nm, we also used the radioactive sources to probe the detector behavior at that scintillation peak. As the optical properties were primarily distance dependent, a method of obtaining this dependence was developed using the differing PMT positions of the height differentiated rows. Defining row 1 as that nearest the top of the detector, we created a set of 4 row ratios for each calibration run. These ratios were defined as the total light seen by each of the other 4 rows over that seen by row 1 over the course of a run. An example of the raw values used for this calculation are shown in Figure 13. Between the two sources and five laser positions and two frequencies and with the addition of a fifth internal ratio -light seen by uncoated tubes versus that of coated -this gave us a total of 60 row ratios to use for the calibrations of the OM. These row ratios were used to compare between the data and simulation outputs through a χ 2 analysis [84]. The process used for determining the parameter values involved generating simulations with varying OM parameters. Minimizing the χ 2 over the row ratios of a simulation would be used to determine the best OM. A range of simulations with χ 2 within 1σ of the minimum χ 2 would be taken to determine the OM uncertainty from the analysis. The χ 2 value of each passing simulation was used in combination with the value of the parameters for that simulation to produce a weighted best fit for each parameter. An overall best fit, including covariance and correlations between parameters, was also produced from these passing simulations. The generated optical model was used to determine the systematic error on potential signals. To determine these systematic errors, we generated signal events -high energy gammas for ALP -across a range of energies and OM uncertainties. As shown in Figure 14, the simulation OM produced a spread of results for the reconstructed energy even for a single true energy particle. This spread is a result of a number of factors, most notably significant differences in energy scale depending on the position inside the detector. However, this spread also includes the spread from OM uncertainty, as the OM parameters were varied in successive simulations according to the uncertainty from the best fit. We used the simulated events across the expected range of true energies -from 1 to 100 MeV -to build a smearing matrix (Figure 15) that would correlate between visible energy in the detector and reconstructed energy for the signal types. The smearing matrix was then applied to potential visible energy distributions as shown in Figure 4 to determine how a signal excess would be reconstructed. These reconstructed signals could be compared directly to the measured excess from the subtraction between ROI data and predicted background. Example excesses are shown overlaid onto the subtraction in Figure 10. Theses signal excesses include the OM uncertainty in the smearing matrix and thus the systematic error on potential signals.
FIG. 14. The reconstructed energy distribution from a large number of correlated OM simulations of a 1 MeV ALP event in CCM. The large width of the distribution represents the spread of outcomes from OM uncertainty, energy resolution, position resolution, and other sources of error. All these factors were encapsulated in the overall OM and OM covariance matrix, allowing the OM to be used for predicting signal efficiency and systematic uncertainty.
These systematic errors ranged from 30% at high energies to 50% at low energies. The primary contributors to this error are position variation, reconstruction resolution, and OM uncertainty. Position variation provides an uncertainty of 20-40% depending on energy. Reconstruction resolution is ±10 cm for position and ±15% for energy, adding another 20% uncertainty. Finally, OM uncertainty increases these uncertainties by around 10%. These uncertainties are added in quadrature to provide the final total systematic error on the signal.
FIG. 15. The generated smearing matrix for the ALP search across the expected potential ALP energies. The smearing matrix includes the efficiency for each energy, with the total for each column being equal to the total efficiency for that energy ALP.