Search for photoproduction of axion-like particles at GlueX

We present a search for axion-like particles, $a$, produced in photon-proton collisions at a center-of-mass energy of approximately 4 GeV, focusing on the scenario where the $a$-gluon coupling is dominant. The search uses $a\to\gamma\gamma$ and $a\to\pi^+\pi^-\pi^0$ decays, and a data sample corresponding to an integrated luminosity of 168 pb$^{-1}$ collected with the GlueX detector. The search for $a\to\gamma\gamma$ decays is performed in the mass range of $180<m_a<480$ MeV, while the search for $a\to\pi^+\pi^-\pi^0$ decays explores the $600<m_a<720$ MeV region. No evidence for a signal is found, and 90% confidence-level exclusion limits are placed on the $a$-gluon coupling strength. These constraints are the most stringent to date over much of the mass ranges considered.

decays explores the 600 < ma < 720 MeV region. No evidence for a signal is found, and 90% confidence-level exclusion limits are placed on the a-gluon coupling strength. These constraints are the most stringent to date over much of the mass ranges considered.

I. INTRODUCTION
Axion-like particles (ALPs), a, are hypothetical pseudoscalars found in many proposed extensions to the Standard Model (SM) [1][2][3][4]. ALPs can naturally address the Strong CP [5][6][7][8] and Hierarchy problems [9], and could provide a portal that connects SM particles to dark matter [10][11][12][13]. The couplings of ALPs to the SM are highly suppressed at low energies by a large cutoff scale Λ; however, since ALPs are pseudo-Nambu-Goldstone bosons, their mass m a can be much smaller than the scale that controls their dynamics, i.e., m a Λ. Recently, ALPs with Λ QCD -scale masses whose dominant coupling to the SM is with the gluonic field have received considerable interest [14][15][16][17][18][19][20][21]. Such m a values can be obtained, e.g., by introducing a mirror strongly-coupled sector that generates a large ALP potential which aligns with the QCD-generated axion potential [22][23][24][25]. In this scenario, the ALP solves the Strong CP problem, and because m a Λ QCD , it is robust against UV contributions that would otherwise give rise to the well-known Quality problem [26][27][28][29].
The effective Lagrangian describing ALP-gluon interactions is where c g is the dimensionless agg vertex coupling constant andG µν ≡ 1 2 µναβ G αβ . Reference [19] presented a novel data-driven method for determining the hadronic interaction strengths of such ALPs, and showed that the dominant decay modes for m a 0.8 GeV (c = 1 throughout this article) are to the γγ, 3π, and π + π − γ final states. In a follow-up article [20], photoproduction of ALPs was explored, including the discovery potential of the GlueX experiment. In addition, Ref. [20] presented a search based on the published GlueX result [30] on γp → pγγ that set world-leading limits in some regions of the ALP parameter space.
In this article, we present a search for ALPs produced in photon-proton collisions at a center-of-mass energy √ s ≈ 4 GeV, focusing on the scenario where the a-gluon coupling is dominant. The search uses a → γγ and a → π + π − π 0 decays, and a data sample corresponding to an integrated luminosity of 168 pb −1 collected with the GlueX detector in the beam-energy range 8 < E beam < 9 GeV. The search for a → γγ decays is performed in the mass range of 180 < m a < 480 MeV, while the search for a → π + π − π 0 decays explores the 600 < m a < 720 MeV * Current address: University of York, York YO10 5DD, United Kingdom region. These mass ranges are chosen to avoid the narrow peaks from the π 0 , η, and ω mesons, and focus on where GlueX is expected to have world-leading sensitivity. Only the charged 3π ALP final state is considered, since the a → 3π 0 decay is more challenging to study using the GlueX detector. Similarly, the other dominant decay mode predicted by Ref. [19] at low masses, namely a → π + π − γ, is not considered due to the huge background from final states such as γp → pπ + π − π 0 with an undetected photon, or γp → pπ + π − with an additional photon from hadronic split-off interactions of the charged pions in the calorimeter. Following Ref. [19], we denote the m a -dependent mixing terms between the ALP and pseudoscalar mesons as aπ 0 and aη . The expected ALP yield in final state F in a small bin of [s, t], where t ≡ (p a − p beam ) 2 , is related to the observed π 0 → F and η → F yields in the bin, n π 0 (s, t) and n η (s, t), by [20] n a (s, where f π and f a ≡ −Λ/32π 2 c g are the pion and ALP decay constants, B(π 0 , η → F) are the known meson-decay branching fractions [31], denotes the m a -dependent product of the detector acceptance and efficiency, and B(a → F) is the ALP-decay branching fraction [19]. Equation (2) assumes that t-channel processes are dominant, which is known to be true at GlueX energies for −t 1 GeV 2 [30]. For the π + π − π 0 decay, the ALP-pion mixing term is negligible and is ignored.
The use of Eq. (2) facilitates a largely data-driven search where most experimental systematic uncertainties cancel. For example, knowledge of the luminosity and absolute efficiencies is not required; only the relative efficiency to reconstruct the ALP and pseudoscalar-meson decays to the same final state is needed. The fiducial regions used, defined in Table I, ensure that the detector response is sufficiently model independent and that t-channel production processes are dominant. This not only reduces the systematic uncertainties that contribute to the results presented here, but enables recasting our results for other models, for which sufficient information has been provided as Supplemental Material [32].

II. EXPERIMENT & SIMULATION
This search is performed using the GlueX spectrometer located in Hall D at Jefferson Lab. A tagged linearly polarized photon beam is created from the 11.6 GeV electron beam at the Continuous Electron Beam Accelerator I. Fiducial regions of the searches for both a → γγ and a → π + π − π 0 decays, defined in terms of the photon beam energy, E beam , the momentum of the outgoing proton, pp, Mandelstam t, the polar angle of the outgoing photons, θγ, the opening angle between the outgoing photon pair, α(γ1, γ2), and the energy of the outgoing photons, Eγ.
All searches Facility [33] by coherent bremsstrahlung off a diamond radiator. This search does not use the polarization information. The momenta of the scattered beam electrons are measured using a scintillating-fiber array, which enables measuring the photon-beam energy to about 0.1% precision in the energy region used in this search. The photon beam passes through a collimator that highly suppresses its incoherent component, before reaching the liquid hydrogen target. This collimation results in most of the high-energy flux being in the small energy range considered here, which motivates our choice of using only a single s bin. The GlueX spectrometer [34], which has a nearly hermetic angular coverage, consists of the liquid hydrogen target surrounded by a scintillator start counter [35], a straw-tube central drift chamber [36], and a lead and scintillating-fiber barrel calorimeter [37], all of which are inside the bore of a superconducting solenoid. In the forward region downstream of the central drift chamber, there are four sets of planar wire drift chambers [38] inside of the solenoid, with a time-of-flight scintillator wall and a forward lead-glass calorimeter [39] located further downstream and outside of the solenoid. The trigger selects events with sufficient energy deposited in the calorimeters, while the event selection for the reactions of interest is performed offline. The drift chambers provide momentum and energy loss measurements for charged particles. The calorimeters measure the energy and position of electromagnetic showers induced by both charged and neutral particles. Finally, particle identification and beam-photon selection are based on timeof-flight measurements obtained using information from the start counter, the calorimeters, and the time-of-flight wall.
Simulation is required to model the effects of the GlueX detector acceptance and its response to the reactions studied here. Specifically, the relative efficiency to reconstruct the ALP and pseudoscalar-meson decays to the same final state is required in Eq. (2). In the simulation, the reactions are generated using the Genr8 package [40], which assumes purely t-channel production.
The effects of the interactions of other beam photons, including accidental coincidences and random detector backgrounds, are included in the simulation. The interaction of the generated particles with the GlueX detector, and its response, are implemented using the Geant4 toolkit [41].

III. EVENT SELECTION
The search is based on events containing candidates for the exclusive reactions γp → pγγ or γp → pπ + π − π 0 , with the subsequent decay π 0 → γγ. Candidate reactions must satisfy all of the fiducial requirements in Table I. In order to avoid experimenter bias, all aspects of the selection are fixed without examining the evidence for an ALP signal.
Events must have the exact number of positively and negatively charged tracks required for each reaction, i.e. events with additional tracks are discarded. Events must also have at least one tagged beam photon candidate, and at least two neutral shower candidates. Additional tagged beam photons and showers are permitted in the initial selection; however, the total energy of the latter is required to be less than 100 MeV.
Charged particles, namely protons and charged pions, are identified using time-of-flight information, along with the energy loss of their tracks in the drift chambers. The absolute value of the squared missing mass for each reaction is required to be less than 0.05 GeV 2 . For a → π + π − π 0 candidates, the diphoton invariant mass must be consistent with that of the π 0 meson. A kinematic fit is used to select particle combinations that are consistent with conservation of energy and momentum, and subsequently, to improve the experimental resolution by enforcing these conservation laws (see Ref. [34] for details). Finally, the photon-proton collision must be consistent with having occurred within the liquid hydrogen target, as determined using the tracking information from the final-state charged particles.
The use of time-of-flight information to perform charged-particle identification enforces that beamphoton candidates have an arrival time, based on the precise electron-beam timing information, at the photonproton collision point that is consistent with when the final-state particles were produced, obtained using information recorded by the GlueX spectrometer. However, incorrect beam-photon candidates that satisfy both the timing and kinematic requirements do occur, typically due to additional photons in the same beam bunch with similar energies. This accidental background is statistically subtracted using a data sample of out-of-time beam photons. Figure 1 shows the γγ and π + π − π 0 invariant mass spectra obtained after applying the full selection. Fits are performed to these spectra in bins of t to obtain the observed π 0 and η yields needed in Eq. (2); Fig. 1 shows the fit results integrated over t. It is worth noting that the bump-hunt procedure we use to search for ALP signals, described in detail in the next section, introduces additional complexity to the background models beyond what is used here. The fits shown in Fig. 1 are only used to determine the meson yields and the mass resolution. For the γγ channel, the fit model consists of the following components. The π 0 and η components are modeled by sums of Gaussian and Crystal Ball [42] functions with power-law tails on both sides of the peaks. The ω → π 0 γ → γγ(γ) component, where one of the three photons is not detected, is modeled by a sum of two Crystal Ball functions with power-law tails on opposite sides of the peak. A combinatorial background component is modeled by a linear function. The π 0 and η yields are 4.4 ± 0.1 million and 0.62 ± 0.02 million, respectively.
Here, the systematic uncertainties, which are obtained by varying both the pseudoscalar and background models, are dominant. The variations include using different numbers of Gaussian and Crystal Ball functions for the pseudoscalar models and using different orders of the polynomial for the background model. In addition, the γγ mass resolution, defined as half the width of the region containing 68% of the signal probability, is determined to be about 6 (9) MeV at m π 0 (m η ). Simulation is used to interpolate between these values to obtain the resolution in the m γγ region considered in the ALP search with a precision of 2%.
The components used in the fit model for the π + π − π 0 channel are as follows. The η component is modeled by a double Gaussian. The known ω lineshape is taken from Ref. [43] and then convolved with a resolution function modeled by a sum of Crystal Ball functions with power-law tails on both sides of the ω peak. A combinatorial background component is modeled by a powerlaw function of the energy released in the π + π − π 0 system. The η yield is found to be 70 ± 1 thousand, where the stated uncertainty is dominated by systematic errors. The π + π − π 0 mass resolution is determined to be about 6 (11) MeV at m η (m ω ). Simulation is again used to interpolate between these values to obtain the resolution in the m π + π − π 0 region considered in the ALP search with 2% precision.

IV. SIGNAL SEARCHES
The signal-search strategy and method, which were first introduced in Ref. [44], are similar to those used in Refs. [45][46][47]. The mass spectra for ALP final states F = γγ and π + π − π 0 are scanned in steps of about half the mass resolution, σ(m F )/2, searching for ALP contributions. At each m a hypothesis, a binned extended maximum-likelihood fit is performed in a ±12.5σ(m γγ ) or ±7.5σ(m π + π − π 0 ) window around m a ; a narrower window is used in the π + π − π 0 final state due to the small distance between the η and ω peaks compared to σ(m π + π − π 0 ). The profile likelihood method is used to determine the local p-values and the ALP signal-yield confidence intervals. The trial factors are obtained using pseudoexperiments for each final state, i.e. by generating a large ensemble of background-only datasets, and rerunning the full signal-search procedure on each sample. The bounded likelihood approach [48] is used when determining the confidence intervals, which defines ∆ log L relative to zero signal, instead of the best-fit value, if the best-fit signal value is negative. This approach has two benefits: it enforces that only physical (nonnegative) upper limits are placed on the ALP yields, and it prevents these limits from being much better than the experimental sensitivity if a large deficit in the background yield is observed.
The ALP signal mass distributions are well modeled by a Gaussian function, whose resolution is determined precisely as described in the previous section. The mass-resolution uncertainty is included in the profile likelihood. A small correction is applied to remove the bias due to neglecting non-Gaussian components of the signal shape, which is determined to be about 1% from the large π 0 and η peaks observed in the data.
The background models include the meson components described in the previous section and shown in Fig. 1. In addition, Legendre polynomial terms up to = 4 (2) for F = γγ (π + π − π 0 ) are taken as inputs, then the datadriven model-selection process of Ref. [44] is performed. This approach is necessary due to the unknown origins of much of the backgrounds. The uncertainty of the model-selection process is included in the profile likelihood following Ref. [49]. Specifically, the aic-o method in Ref. [44] is used, which penalizes the log-likelihood of each background model according to its complexity (number of parameters). The confidence intervals are obtained from the penalized profile likelihoods treating the model index as a discrete nuisance parameter [49].
The maximum values are chosen for each final state to ensure adequate description is possible for any peaking background that may contribute. Where such complexity is unnecessary, the data-driven model-selection procedure effectively reduces the complexity to increase the sensitivity because the associated penalty terms result in overly complex models being ignored when constructing the confidence intervals (see Ref. [44] for more detailed discussion). Following Ref. [44], all fit regions are transformed onto the interval [−1, 1] with the scan m a value mapped to zero. The signal model is an even function after this transformation; therefore, the presence of odd Legendre modes, which are orthogonal to the signal component, has minimal impact on the variance of the ALP yield. Thus, all odd Legendre modes are included in every background model, while only a subset of the even modes is selected. This procedure results in a mass-dependent background-model uncertainty, which on average is about 5%. Figure 2 shows the signed local significances for all m a values scanned for both ALP decays considered. The largest local excess in the γγ spectrum is 2.3σ at 213 MeV. Similarly, the largest local excess in the π + π − π 0 spectrum is 1.6σ at 669 MeV. The global p-value is found to be 0.28 (0.50) for the γγ (π + π − π 0 ) channel, after accounting for the trials factor due to the number of signal hypotheses.

V. ACCEPTANCE & EFFICIENCY
The use of Eq. (2) for normalization means that only the relative efficiency to reconstruct the ALP and pseudoscalar-meson decays to the same final state is needed. Furthermore, since the normalization is done in narrow [s, t] bins, the production mechanisms do not need to be well understood. The acceptance is defined here as the probability that a reaction producing an ALP in a [s, t] bin will have all final-state particles in the fidu- Signed local significances at each scan mass from (black points) all fits and (red lines) the expected distribution for the (top) γγ and (bottom) π + π − π 0 channels; if the bestfit signal-yield estimator is negative, the signed significance is negative and vice versa. The error bars account for the correlation between nearest-neighbor fit results, which often produces outliers in pairs due to the σ(m)/2 step size. cial region defined in Table I. This acceptance is strongly dependent on m a and requires careful treatment as described below. For accepted reactions, the reconstruction efficiency has minimal dependence on m a or t and is taken from simulation. Indeed, our choice of fiducial region is designed to minimize the m a and t dependence, since only this dependence enters into Eq. (2). The uncertainty due to the relative reconstruction efficiency is, therefore, negligible compared to that of the acceptance.
The acceptance varies strongly with both m a and t. The t bins used for the normalization are only 0.05 GeV 2 wide, which reduces the impact of the t dependence, but it still must be accounted for. To do this, we numerically sample from the phase-space for each t bin and obtain the probability that the reaction satisfies the fiducial requirements. In each t bin, we consider three scenarios: (i) t fixed to the lower edge of the bin, (ii) t fixed to the upper edge of the bin, and (iii) t sampled uniformly over the bin range. Since the bins are narrow, scenario iii is used to obtain the nominal acceptance. Half the difference of i or ii from nominal, whichever difference is larger is used, is taken as the systematic uncertainty in the acceptance in each t bin. Bins whose acceptance uncertainty is larger than 10% are excluded from the fiducial region. The product of the acceptance and efficiency for each ALP decay as a function of m a and t is provided in Supplemental Material. The acceptance uncertainties in each bin, obtained as described in the previous paragraph, are propagated to the expected ALP yield using Eq. (2) for each m a value considered in the search. These uncertainties, which vary with m a , are typically about 5%.

VI. RESULTS
The upper limits on the ALP signal yields obtained in Sec. IV are normalized using Eq. (2), which takes as input the pseudoscalar-meson yields from Sec. III, relative efficiency from Sec. V, the pseudoscalar-meson decay branching fractions from the PDG [31], and the ALP decay branching fractions from Ref. [19]. The systematic uncertainties on the ALP signal yield, the pseudoscalarmeson yields and branching fractions, and on the relative efficiency are included in the profile likelihood when determining the upper limits on the ALP yield. These uncertainties, which were described previously, are summarized in Table II. ALPs are excluded at 90% confidence level (CL) when the upper limit on the observed ALP yield is less than the expected ALP yield in Eq. 2. Figure 3 shows the constraints placed on the ratio of ALP parameters c g /Λ for each final state. Taking c g to be O(1), our results correspond to O(TeV) constraints on Λ. Figure 4 compares our results to the best existing constraints on the ALP-gluon coupling. The kaon-decay and B-lifetime constraints involve penguin decays that proceed via loops that are sensitive to the unknown UV completion of the full theory. The constraints shown in Fig. 4 are taken from Refs. [19,50] which assume an O(TeV) UV scale, though both references note that these constraints have O(1) uncertainties induced by the unknown UV physics. The searches presented here place more robust limits-which are also the most stringent constraints over much of the mass ranges considered. These results demonstrate the power of using high-intensity photon beams to search for low-mass physics beyond the Standard Model.

VII. SUMMARY
In summary, we present a search for axion-like particles produced in photon-proton collisions at a center-ofmass energy of approximately 4 GeV. The search looked for a → γγ and a → π + π − π 0 decays in a data sample corresponding to an integrated luminosity of 168 pb −1 collected with the GlueX detector. The search for a → γγ decays was performed in the mass range of 180 < m a < 480 MeV, while the search for a → π + π − π 0 decays explored the 600 < m a < 720 MeV region. No evidence for an ALP signal was found, leading to 90% confidence-level exclusion limits on the ALP-gluon coupling strength. These constraints are the most stringent to date over much of the mass ranges considered.  [19] from LEP [51,52], φ and η decays [31], and the (cyan line) GlueX limits recast from Ref. [30]. In addition, limits obtained from kaon decays [19,50,[53][54][55] and the B-meson lifetime [19,50], which have O(1) uncertainties induced by the unknown UV physics, are shown as hashed regions.