Search for Axion-Like Particles produced in $e^+e^-$ collisions at Belle II

We present a search for the direct production of a light pseudoscalar $a$ decaying into two photons with the Belle II detector at the SuperKEKB collider. We search for the process ${e^+e^-\to\gamma a, a \to\gamma\gamma}$ in the mass range ${0.2} \,<m_a<{9.7}\,{\text{GeV/$c$}^2}$ using data corresponding to an integrated luminosity of $(445\pm 3)\,\text{pb}^{-1}$. Light pseudoscalars interacting predominantly with standard model gauge bosons (so-called axion-like particles or ALPs) are frequently postulated in extensions of the standard model. We find no evidence for ALPs and set 95% confidence level upper limits on the coupling strength $g_{a\gamma\gamma}$ of ALPs to photons at the level of $10^{-3}\,{\text{GeV}^{-1}}$. The limits are the most restrictive to date for $0.2\,<\,m_a\,<\,5\,{\text{GeV/$c$}^2}$.

We present a search for the direct production of a light pseudoscalar a decaying into two photons with the Belle II detector at the SuperKEKB collider. We search for the process e + e − → γa, a → γγ in the mass range 0.2 < ma < 9.7 GeV/c 2 using data corresponding to an integrated luminosity of (445 ± 3) pb −1 . Light pseudoscalars interacting predominantly with standard model gauge bosons (so-called axion-like particles or ALPs) are frequently postulated in extensions of the standard model. We find no evidence for ALPs and set 95% confidence level upper limits on the coupling strength gaγγ of ALPs to photons at the level of 10 −3 GeV −1 . The limits are the most restrictive to date for 0.2 < ma < 5 GeV/c 2 .
Axions and axion-like particles (ALPs) are predicted by many extensions of the standard model (SM) [1]. They occur, for example, in most solutions of the strong CP problem [2]. ALPs share the quantum numbers of axions, but differ in that their masses and couplings are independent. ALPs with sub-MeV/c 2 masses are interesting in the context of astrophysics and cosmology and are cold dark matter (DM) candidates, whereas ALPs with O(1 GeV/c 2 ) masses generally relate to several topics in particle physics [3][4][5]. Most notably, heavy ALPs can connect the SM particles to yet undiscovered DM particles [6]. ALPs that predominantly couple to γγ, γZ 0 , and Z 0 Z 0 are experimentally much less constrained than those that couple to gluons or fermions. The latter interactions typically lead to flavor-changing processes that can be probed in rare decays [7]. In this letter we will consider the case that the ALP a predominantly couples to photons, with coupling strength g aγγ , and has negligible coupling strength g aγZ to a photon and a Z 0 boson, so that B(a → γγ) ≈ 100%; we follow the notation for couplings introduced in Ref. [6]. In the MeV/c 2 to GeV/c 2 mass range, the current best limits for ALPs with photon couplings are derived from a variety of experiments. These limits come from e + e − → γ + invisible and beam-dump experiments for light ALPs [6,8,9], from e + e − → γγ [10] and coherent Primakoff production off a nuclear target [11] for intermediate-mass ALPs, and from peripheral heavy-ion collisions [12] for heavy ALPs.
We search for e + e − → γa, a → γγ in the ALP mass range 0.2 < m a < 9.7 GeV/c 2 in the three-photon final state. The signature in the center-of-mass (c.m.) system is a monoenergetic photon recoiling against the a → γγ decay. The energy of the recoil photon is where √ s is the c.m. collision energy. We search for an ALP signal as a narrow peak in the squared recoilmass distribution M 2 recoil = s − 2 √ sE c.m. recoilγ , or as a narrow peak in the squared-invariant-mass distribution M 2 γγ , computed using the two-photon system, depending on which provides the better sensitivity. We note that in the future a larger Belle II dataset will be available to calibrate the photon covariance matrix, which in turn will allow the use of kinematic fitting of the three photons to the known beam four-momentum, thus improving the sensitivity. In our search range, the width of the ALP is negligible with respect to the experimental resolution, and the ALP lifetime is negligible, thus it decays promptly. The dominant SM background process is e + e − → γγγ. The analysis selection, fit strategy, and limit-setting procedures are optimized and verified based on Monte Carlo simulation, i.e. without looking at data events, to avoid experimenter's bias.
We use a data set corresponding to an integrated lu-minosity of (496 ± 3) pb −1 [13] collected with the Belle II detector at the asymmetric-energy e + e − collider Su-perKEKB [14], which is located at the KEK laboratory in Tsukuba, Japan. Data were collected at the c.m. energy of the Υ(4S) resonance ( √ s = 10.58 GeV) from April to July 2018. The energies of the electron and positron beams are 7 GeV and 4 GeV, respectively, resulting in a boost of βγ = 0.28 of the c.m. frame relative to the laboratory frame. We use a randomly-chosen subset of the data, approximately 10%, to validate the selection, and we then discard it from the final data sample. The remaining data set is used for the search and corresponds to an integrated luminosity of (445 ± 3) pb −1 . The Belle II detector consists of several subdetectors arranged around the beam pipe in a cylindrical structure [15,16]. Only the components that are relevant to this analysis are described below. Photons are measured and identified in the electromagnetic calorimeter (ECL) consisting of CsI(Tl) crystals. The ECL provides both an energy and a timing measurement. A superconducting solenoid situated outside of the calorimeter provides a 1.5 T magnetic field. Charged-particle tracking is done using a silicon vertex detector (VXD) and a central drift chamber (CDC). Only one azimuthal octant of the VXD was present during the 2018 operations. The z-axis of the laboratory frame coincides with that of the solenoid and its positive direction is approximately that of the incoming electron beam. The polar angle θ is measured with respect to this direction. Events are selected only by the hardware trigger, and no further software trigger selection is applied. Trigger energy thresholds are very low and no vetoes for abundant QED scattering processes are applied.
We use Babayaga@nlo [17][18][19][20] to generate SM background processes e + e − → e + e − (γ), e + e − → γγ(γ). We use Phokhara9 [21] to generate SM background processes e + e − → P γ(γ), where P is a SM pseudoscalar meson (π 0 , η, η ). This includes production via the radiative decay of the intermediate vector resonances ρ, ω, and φ. The largest pseudoscalar background contribution for this analysis comes from e + e − → ωγ, ω → π 0 γ with a boosted π 0 decaying into overlapping photons. We use the same generators to calculate the cross sections of the respective processes. We use MadGraph5 [22] to simulate signal events, including the effects of initial-state radiation (ISR) in event kinematics [23], for different hypotheses for m a in step sizes approximately equal to the signal resolution in our search range.
We use Geant4 [24] to simulate the interactions of particles in the detector, taking into account the nominal detector geometry and simulated beam-backgrounds adjusted to match the measured beam conditions. We use the Belle II software framework [25] to reconstruct and analyze events.
All selection criteria are chosen to maximize the Punzi figure of merit for 5 σ discovery [26]. Quantities are defined in the laboratory frame unless otherwise specified. Photon candidates are reconstructed from ECL clusters with no associated charged tracks. We select events with at least three photon candidates with energy E γ above 0.65 GeV (for m a > 4 GeV/c 2 ) or 1.0 GeV (for m a ≤ 4 GeV/c 2 ). This ALP-mass-dependent threshold is used to avoid shaping effects on the background distribution in the mass fit range. The following selection variables are not dependent on the ALP mass. All three photon candidates must be reconstructed with polar angles 37.3 < θ γ < 123.7 • . This polar-angle region provides the best calorimeter energy resolution, avoids regions close to detector gaps, and offers the lowest beam background levels. If more than three photons pass the selection criteria, we select the three most energetic ones and the additional photons are ignored in the calculation of any variables. This occurs in fewer than 0.2% of all events. We reduce contamination from beam backgrounds by requiring that each photon detection time t i is compatible with the average weighted photon timē where ∆t i is the energy-dependent timing range that includes 99% of all signal photons, and is between 3 ns (high E γ ) and 15 ns (low E γ ). The requirement is |(t i −t)/∆t i | < 10, which is insensitive to global time offsets. The invariant mass M γγγ of the three-photon system must satisfy 0.88 √ s ≤ M γγγ ≤ 1.03 √ s to eliminate kinematically unbalanced events coming from cosmic rays, beam-gas backgrounds, or two-photon production. We reject events that have tracks originating from the interaction region to suppress background from e + e − → e + e − γ. We require a θ γ separation between any two photons of ∆θ γ > 0.014 rad, or an azimuthal angle separation of ∆φ γ > 0.400 rad to reduce background from photon conversions outside of the tracking detectors. Following a data-sideband analysis using M γγγ < 0.88 √ s, we additionally apply a loose selection, based on a multivariate shower-shape classifier that uses multiple Zernike moments [27], on the most isolated of the three photons. This criterion reduces the number of clusters produced by neutral hadrons and by particles that do not originate from the interaction point. The selection procedure results in three ALP candidates per event from all possible combinations of the three selected photons.
The resulting M 2 recoil and M 2 γγ distributions are shown in Fig. 1 together with the stacked contributions from the luminosity-normalized simulated samples of SM backgrounds. The expected background distributions are dominated by e + e − → γγγ with a small contribution from e + e − → e + e − γ due to tracking inefficiencies. We find contributions from cosmic rays, assessed in data-taking periods without colliding beams, neither significant nor peaking in photon energy or invariant mass. The data shape agrees well with simulation except for a small and localized excess seen in the low-mass region M 2 γγ < 1 GeV 2 /c 4 . The excess is broad (see the inset in Fig. 1 (b)) and not consistent with an ALP signal, for which we expect a much smaller width in this region (see the inset in Fig. 2). As described later, the signal extraction does not directly depend on the background predictions because we fit the background only using data, thus any discrepancy between data and simulation has little impact on the result. Triggers based on 1 GeV threshold energy sums in the calorimeter barrel are found to have ε trg = 1.0 for the ALP selection, based upon studies of radiative Bhabha events. The ALP selection efficiency is determined using large simulated signal samples, and varies smoothly between 20% (low m a ) and 34 % (high m a ). The number of candidates in data is 3.6 ± 0.9% (4.2 ± 1.1%) higher than in the simulation for the E γ > 0.65 GeV (E γ > 1.0 GeV) selection. No correction is applied and we assign the sum of the full difference and its uncertainty as a systematic uncertainty for the selection efficiency. We assess the difference in the photon-energy reconstruction between data and simulation by using radiative muon-pair events in which we compare the predicted recoil energy calculated from the muon-pair momenta with the energy of the photon candidate. We correct for the observed linear energy bias that ranges from 0 (low energy) to 0.5% (high energy). We vary the energy selection by ±1% and the angular-separation selection by the approximate position resolution of ± 5 mrad, and take the respective full difference in the signal selection efficiency with respect to the nominal selection as a systematic uncertainty. We add these three uncertainties in quadrature assuming no correlations amongst them. The total relative uncertainty due to the selection efficiency is approximately 5.5% for ALP masses above 0.5 GeV/c 2 , and increases to approximately 8% for the lightest ALP masses considered. As additional systematic checks we vary the photon-timing selection by ±1 and the shower-shape classifier selection by ± 5% to account for possible between data and simulation samples, the invariant mass M γγγ selection by ± 0.002 GeV/c 2 to account for uncertainties in the beam energy, and the polar-angle-acceptance selection by propagating the effect of a ±2 mm shift of the interaction point relative to the calorimeter to account for maximal possible misalignment of the ECL. For all of these checks, we find that they have a negligible effect on the signal selection efficiency, so we do not associate any systematic uncertainty with them.
We extract the signal yield as a function of m a by performing a series of independent binned maximum-likelihood fits. We use 100 bins for each fit range.
The fits are performed in the range 0.2 < m a < 6.85 GeV/c 2 for the M 2 γγ spectrum, and in the range 6.85 < m a < 9.7 GeV/c 2 for the M 2 recoil spectrum. The resolution of M 2 γγ worsens with increasing m a , while that of M 2 recoil improves with increasing m a (see Fig. 2). The transition between M 2 γγ and M 2 recoil fits is determined as the point of equal sensitivity obtained using background simulations.
The signal probability density function (PDF) has two components: a peaking contribution from correctly reconstructed signal photons and a combinatorialbackground contribution from the other two combinations of photons. We model the peaking contribution using a Crystal Ball (CB) function [28]. The massdependent CB parameters used in the fits to data are fixed to those obtained by fitting simulated events. For the simulated M 2 recoil distribution, the CB mean is found to be unbiased. For the simulated M 2 γγ distribution, we observe a linear bias of the CB mean of about 0.5% resulting from the combination of two photons with asymmetric reconstructed-energy distributions. This bias is determined to have negligible impact on the signal yield and mass determination; therefore, no attempt to correct for it is made. Combinatorial-background contributions from the wrong combinations of photons in signal events are taken into account by adding a massdependent, one-dimensional, smoothed kernel density estimation (KDE) [29] PDF obtained from signal simulation. The fits are performed in steps of m a that correspond to half the CB width (σ CB ) for the respective squared mass. This results in a total of 378 fits to the M 2 γγ distribution and 124 fits to the M 2 recoil distribution. CB signal parameters are interpolated between the known simulated masses, and the KDE shape is taken from the simulation sample generated with the closest value of m a to that assumed in the fit.
The photon-energy resolution σ(E γ )/E γ in simulation is about 3% for E γ = 0.65 GeV and improves to about 2% for E γ > 1 GeV. Using the same muon-pair sample as used for the photon-energy bias study, we find that the photon energy resolution in simulation is better than that in data by at most 30% at low energies. Therefore, we apply an energy-dependent additional resolution smearing to our simulated signal samples before determining the CB resolution parameter σ CB ; we assume conservatively that the full observed difference between data and simulation is due to the photon-energy-resolution difference. We assign half of the resulting mass-resolution difference as a systematic uncertainty. The effect of a ±2 mm shift of the interaction point relative to the calorimeter is found to have a negligible impact on the the mass resolution and is not included as a systematic uncertainty.
We describe the backgrounds by polynomials of the minimum complexity consistent with the data features. Polynomials of 2nd to 5th order are used: 2nd for 0.2 < m a ≤ 0.5 GeV/c 2 , 4th for 0.5 < m a ≤ 6.85 GeV/c 2 , and 5th for 6.85 < m a ≤ 9.7 GeV/c 2 . The background poly-nomial parameters are not fixed by simulation but are free parameters of each data fit. Each fit is performed in a mass range that corresponds to −20 σ CB to +30 σ CB for M 2 γγ , and −25 σ CB to +25 σ CB for M 2 recoil . In addition, the fit ranges are constrained between M 2 γγ > 0 GeV 2 /c 4 and M 2 recoil < 100.5 GeV 2 /c 4 . The choice of the order of background polynomial and fit range is optimized based on the following conditions: giving a reduced χ 2 close to one, providing locally smooth fit results, and being consistent with minimal variations between adjacent fit ranges. Peaking backgrounds from e + e − → P γ are very small compared to the expected statistical uncertainty on the signal yield and found to be modeled adequately by the polynomial background PDF.
The systematic uncertainties due to the signal efficiency and the signal mass resolution are included as Gaussian nuisance parameters with a width equal to the systematic uncertainty. The systematic uncertainty due to the background shape, which is the dominant source of systematic uncertainty, is estimated by repeating all fits with alternative fit ranges changed by ±5σ CB and with the polynomial orders modified by ±1. For each mass value m a , we report the smallest of all signal significance values determined from each background model. The local significance including systematic uncertainties is given by S = 2 ln(L/L bkg ), where L is the maximum likelihood for the fit, and L bkg is the likelihood for a fit to the background-only hypothesis. The local significances, multiplied by the sign of the signal yield, are shown in Fig. 3. The largest local significance, including systematic uncertainties, is found near m a = 0.477 GeV/c 2 with a value of S = 2.8 σ. By dividing the signal yield by the signal efficiency and the integrated luminosity, we obtain the ALP cross section σ a . We compute the 95% confidence level (CL) upper limits on σ a as a function of m a using a one-sided frequentist profile-likelihood method [30]. For each m a fit result, we report the least stringent of all 95% confidence level (CL) upper limits determined from the variations of background model and fit range. We convert the crosssection limit to the coupling limit using where α QED is the electromagnetic coupling [6]. This calculation does not take into account any energy dependence of α QED and g aγγ itself [31]. An additional 0.2% collision-energy uncertainty when converting σ a to g aγγ results in a negligible additional systematic uncertainty. Our median limit expected in the absence of a signal and the observed upper limits on σ a are shown in Fig. 4. The observed upper limits on the photon couplings g aγγ of ALPs, as well as existing constraints from previous experiments, are shown in Fig. 5. Additional plots and numerical results can be found in the Supplemental Material [32]. Our results provide the best limits for 0.2 < m a < 5 GeV/c 2 . This region of ALP parameter space is completely unconstrained by cosmological considerations [33]. The remaining mass region below 0.2 GeV/c 2 is challenging to probe at colliders due to the poor spatial resolution of photons from highly boosted ALP decays, and irreducible peaking backgrounds from π 0 production. In conclusion, we search for e + e − → γa, a → γγ in the ALP mass range 0.2 < m a < 9.7 GeV/c 2 using Belle II data corresponding to an integrated luminosity of 445 pb −1 . We do not observe any significant excess of events consistent with the signal process and set 95% CL upper limits on the photon coupling g aγγ at the level of Upper limit (95% CL) on the ALP-photon coupling from this analysis and previous constraints from electron beam-dump experiments and e + e − → γ + invisible [6,9], proton beam-dump experiments [8], e + e − → γγ [10], a photonbeam experiment [11], and heavy-ion collisions [12].
10 −3 . These limits are almost one order of magnitude more restrictive than existing limits from LEP [10]. In the future, with increased luminosity, Belle II is expected to improve the sensitivity to g aγγ by more than one order of magnitude [6].
We thank the SuperKEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group for on-site computing support. This work was supported by the following funding sources: Sci- Comparisons between data and simulated sample distributions of variables used in the selection and fit are given in Fig. 1-8. Fig. 9 shows the fit for m a = 0.477 GeV/c 2 that resulted in the largest significance observed. We provide an additional textfile with numerical results the observed 95% CL upper limit on the ALP cross section σ a (pb), and the observed 95% CL upper limit on the photon coupling g aγγ (GeV −1 ) as a function of ALP mass m a (GeV).  , signal PDF (black) and total PDF (red). The signal peak corresponds to a local significance including systematic uncertainties of S = 2.8 σ.