Search for an axion-like particle in radiative J/ψ decays

We search for an axion-like particle (ALP) a through the process ψ (3686) → π + π − J/ψ , J/ψ → γa , a → γγ in a data sample of (2 . 71 ± 0 . 01) × 10 9 ψ (3686) events collected by the BESIII detector. No signiﬁcant ALP signal is observed over the expected background, and the upper limits on the branching fraction of the decay J/ψ → γa and the ALP-photon coupling constant g aγγ are set at 95% conﬁdence level in the mass range of 0 . 165 ≤ m a ≤ 2 . 84 GeV /c 2 . The limits on B ( J/ψ → γa ) range from 8 . 3 × 10 − 8 to 1 . 8 × 10 − 6 over the search region, and the constraints on the ALP-photon coupling are the most stringent to date for 0 . 165 ≤ m a ≤ 1 . 468 GeV /c 2 .


Introduction
Axion-like particles (ALPs) are pseudo-Goldstone bosons arising from some spontaneously broken global symmetry, addressing the strong CP [1][2][3][4] or hierarchy problems [5].ALPs could appear in theories beyond the Standard Model (SM), such as string theory [6] and extended Higgs models [7].ALPs could also provide a portal connecting SM particles to the dark sectors [8], and in certain situations, they are proposed as cold dark matter candidates [9][10][11].In the most common scenarios, the ALP a predominantly couples to photons with a photon coupling constant g aγγ .As a generalization of QCD axions, ALPs have arbitrary masses and couplings which are bounded by experiments.In the sub-MeV/c 2 mass range, the bounds on g aγγ are provided by laser experiments, solar photon instruments, as well as cosmological and astrophysical observations [12,13].In the MeV/c 2 to GeV/c 2 mass region, constraints on the photon coupling constant mainly come from beamdump experiments [14][15][16] and high-energy collider experiments.Long-lived ALPs have been searched for in J/ψ and Υ radiative decays [17,18] and through e + e − → γ + invisible processes [19].The limits for short-lived ALPs decaying to two photons are obtained from e + e − → γγ [20], e + e − → γγγ [21,22] via ALP-strahlung production, and the light-by-light scattering process γγ → γγ [23,24], as well as pp collisions [25][26][27].There is a striking gap in the limits on the ALP-photon coupling in the ALP mass range from 100 MeV/c 2 to roughly 10 GeV/c 2 , leaving the large area 10 −3 g aγγ 10 −5 GeV −1 still uncov-ered.Electron-positron colliders can play a unique role in further exploring this region [28].
We search for ALPs decaying into two photons at BESIII in J/ψ radiative decays via J/ψ → γa, a → γγ in the mass range of 0.165 ≤ m a ≤ 2.84 GeV/c 2 .In this Letter we ignore the ALP coupling to c-quarks [29] and assume the branching fraction of a decaying to photons is 100%, with a decay width Γ a = g 2 aγγ m 3 a /64π.The experimental signature of the ALP decaying into two photons depends on the correlation between its mass and coupling.Except for the mass region below 0.1 GeV/c 2 , the ALP has a negligible lifetime and decay width, and the opening angle between the decay photons is large enough to be resolved in the detector.The mass intervals of 0.10 < m a < 0.165 GeV/c 2 , 0.46 < m a < 0.60 GeV/c 2 and 0.90 < m a < 1.01 GeV/c 2 are excluded from the search due to the peaking backgrounds from the decays of π 0 , η and η mesons, respectively.In order to avoid pollution from non-resonant ALP production e + e − → γa and to probe resonant ALP production J/ψ → γa in a model-independent way [29], the decay ψ(3686) → π + π − J/ψ is exploited to isolate the J/ψ sample.

Detector, data sets and Monte Carlo simulation
The BESIII detector [30] records symmetric e + e − collisions provided by the BEPCII storage ring [31], which operates in the center-of-mass energy range from 2.00 to 4.95 GeV, with a peak luminosity of 1 × 10 33 cm −2 s −1 achieved at √ s = 3.77 GeV.BESIII has collected large data samples in this energy region [32].
The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field.The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel.The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for electrons from Bhabha scattering.The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region.The time resolution in the TOF barrel region is 68 ps, while that in the end-cap region is 110 ps.The end-cap TOF system was upgraded in 2015 using multigap resistive plate chamber technology, providing a time resolution of 60 ps [33].
We perform the search using the ψ(3686) data samples collected by the BESIII detector during three datataking periods: 2009, 2012, and 2021.The numbers of ψ(3686) events are determined by counting inclusive hadronic events whose branching fraction is known rather precisely [34,35].There are (107.0± 0.8) × 10 6 ψ(3686) events collected in 2009, (341.1 ± 2.9) × 10 6 events collected in 2012 and (2.26 ± 0.01) × 10 9 events collected in 2021, for a total of (2.71 ± 0.01) × 10 9 ψ(3686) events.In this analysis, the J/ψ sample originates from the decay ψ(3686) → π + π − J/ψ.The analysis strategy is to first tag J/ψ events by selecting two oppositely charged pions, and then to search for the signal candidate events with three additional photons originating from the tagged J/ψ sample.A semi-blind procedure is performed to avoid possible bias, where approximately 10% of the full data set is used to optimize the event selection and validate the fit approach.The final result is then obtained with the full data set only after the analysis strategy is fixed.
Simulated data samples produced with a GEANT4based [36] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate backgrounds.The simulation models the beam-energy spread and initialstate radiation in e + e − annihilations with the generator KKMC [37].An inclusive MC sample with approximately the same number of ψ(3686) decays as in data is used to check for potential backgrounds.The inclusive MC sample includes the production of the ψ(3686) resonance, the initial-state radiation production of the J/ψ, and the continuum processes incorporated in KKMC [37].All particle decays are mod-eled with EVTGEN [38] using branching fractions either taken from the Particle Data Group (PDG) [39], when available, or otherwise estimated with LUND-CHARM [40].Final-state radiation from charged finalstate particles is incorporated using the PHOTOS package [41].To simulate the signal, the decay ψ(3686) → π + π − J/ψ is modeled according to the partial wave analysis results of ψ(3686) → π + π − J/ψ [42]; a P -wave model [38] is used for the subsequent decay J/ψ → γa and a phase-space model is used for the a → γγ decay.The signal events are simulated individually for different mass hypotheses of ALP m a in the range of 0.165 ≤ m a ≤ 2.84 GeV/c 2 with a step size of 5 MeV/c 2 , which is smaller than the signal resolution.Study of the inclusive MC sample with a generic event type analysis tool, TopoAna [43], indicates that the dominant backgrounds are from ψ(3686) → π + π − J/ψ with subsequent decays of J/ψ → γπ 0 , γη and γη .These backgrounds are each generated exclusively with the angular distribution of 1 + cos 2 θ γ [38], where θ γ is the angle of the radiative photon relative to the positron beam direction in the J/ψ rest frame.Some potential backgrounds of the form ψ(3686) → π + π − J/ψ with J/ψ decaying into purely neutral particles in the final states are generated exclusively with different generators: J/ψ → γη c with the angular distribution of 1 + cos 2 θ γ [38], J/ψ → γπ 0 π 0 according to the partial wave analysis results of J/ψ → γπ 0 π 0 [44], as well as J/ψ → γγγ with a phase-space distribution.Each background MC sample is weighted to match the integrated luminosity of the data set.

Event selection
Signal candidates are required to have two oppositely charged tracks and at least three photon candidates.Charged tracks detected in the MDC are required to be within a polar angle (θ) range of |cosθ| < 0.93, where θ is defined with respect to the z-axis, which is the symmetry axis of the MDC.For each charged track, the momentum must be less than 0.45 GeV/c and the distance of closest approach to the interaction point (IP) along the z-axis |V z | must be less than 10 cm, and in the transverse plane |V xy | less than 1 cm.Both charged tracks are assumed to be pion candidates, and the recoil mass in the center-of-mass system M (π + π − ) recoil must be in the range of [3.080, 3.114 Photon candidates are identified using showers in the EMC.The deposited energy of each shower must be more than 25 MeV in the barrel region (|cos θ| < 0.80) and more than 50 MeV in the end-cap region (0.86 < |cos θ| < 0.92).To exclude showers that originate from charged tracks, the angle subtended by the EMC shower and the position of the closest charged track at the EMC must be greater than 10 degrees as measured from the IP.To suppress electronic noise and showers unrelated to the event, the difference between the EMC time and the event start time is required to be within [0, 700] ns.Events satisfying the above requirements are retained for further analysis.
The π + and π − tracks are constrained to a common vertex to determine the event interaction point, and a four-constraint (4C) kinematic fit to the initial four-momentum of the ψ(3686) is applied for all possible π + π − γγγ combinations.The combination with the smallest fit χ 2 4C is retained, with the requirement of χ 2 4C < 40.The distribution of χ 2 4C is shown in Fig. 1.To further reject background from J/ψ → γπ 0 π 0 decays, we form all π + π − 5γ combinations if 5 or more photons are present, and we require the minimum χ 2 4C (5γ) of the corresponding kinematic fits to satisfy χ 2 4C < χ 2 4C (5γ).For the remaining photons, the sum of their energy deposition in the EMC E other is required to be less than 0.1 GeV.The requirements of χ 2 4C and E other are optimized according to the Punzi significance [45] defined as /(1.5+ √ B), where denotes the signal efficiency obtained from signal MC samples and B is the number of background events obtained from background MC samples.4C for data, signal and MC-simulated backgrounds.The requirement of χ 2 4C < 40 is applied to suppress the background from J/ψ → γπ 0 π 0 decays.

ALP signal search
The two-photon invariant mass M γγ distribution of the events selected by the above criteria is shown in Fig. 2. The MC simulation shows a rather good agreement with the data, corroborating the analysis procedure.There are three entries per event from all possible combinations of the three selected photons.The background is dominated by contributions from J/ψ → γπ 0 , γη and γη decays.There is a small contribution from residual J/ψ → γπ 0 π 0 background with two soft photons peaking around the π 0 mass region, and a contamination from J/ψ → γf 2 (1270) decays, followed by f 2 (1270) → π 0 π 0 .However, since the f 2 (1270) state has a large width and the contribution is relatively small, this contribution can be treated as non-peaking background and is well described by a polynomial function when extracting the signal yields.
A series of one-dimensional unbinned extended maximum-likelihood fits are performed to the M γγ distribution to determine the signal yields with different ALP mass hypotheses in the mass range of 0.165 ≤ m a ≤ 2.84 GeV/c 2 .The search step is 3 MeV/c 2 for 0.165 ≤ m a < 1.20 GeV/c 2 and 4 MeV/c 2 for other mass regions.A total of 674 mass hypotheses are probed.The likelihood function is a combination of signal, non-peaking background, and peaking components of the π 0 , η and η mesons.Different M γγ fit intervals are used for various m a hypotheses in order to handle the non-peaking and peaking backgrounds properly, which are listed in Table 1.

FIG. 2:
The diphoton invariant-mass distributions for data and the MC-simulated backgrounds, which are normalized to the integrated luminosity of the data.The signal probability density function (PDF) is constructed from a peaking component from the ALP decay, and a combinatorial contribution from the other two combinations of photons.The peaking component is parameterized by the sum of two Crystal Ball (CB) [46] functions with opposite-side tails with the same mean, resolution, and relative weight.These values and the detection efficiency are determined by fitting the signal MC samples.The resolution ranges from 6 MeV/c 2 near m a = 0.165 GeV/c 2 to 11 MeV/c 2 near m a = 2.2 GeV/c 2 , and it decreases back to 7 MeV/c 2 near m a = 2.84 GeV/c 2 .The detection efficiency ranges from 30% to 35%.The combinatorial component of the signal is described by a smoothed kernel density estimation PDF [47] obtained from the signal MC sample with the mass hypothesis m a closest to the fitted assumption.The ratio between the combinatorial part and the peaking component of the signal is determined for each fit interval.For any search point, the parameters of the CB functions, the ratio of the combinatorial part and the efficiency are interpolated between the mass points by a fit with a polynomial function.The peaking background components at the π 0 , η and η masses are described by the sum of two CB functions using parameters determined from data.A fifth-order Chebyshev polynomial function is used to describe the non-peaking background for 0.75 ≤ m a < 1.20 GeV/c 2 and a third-order Chebyshev polynomial function is used for the remaining mass regions.The choice of the order of polynomial function is optimized through a spurious signal test by performing a signalplus-background fit to the background-only M γγ distribution [48].The obtained signal yield (spurious signal) is required to be less than 30% of its expected statistical uncertainty.When performing the fit, the shapes of the signal and the peaking background PDFs are fixed, while the non-peaking background PDF shape and the yields of signal, peaking and non-peaking background events are free parameters.The M γγ distributions for selected m a points and the fit results are shown in Fig. 3.
The branching fractions of J/ψ → γ(π 0 , η, η ) → γγγ are measured to validate the signal extraction procedure.The peaks are treated as signals, and the fitting procedure described above is applied to extract the number of peak events.After accounting for the contribution of peaking backgrounds, the results are found to be compatible with the published BESIII measurement [49] within uncertainties.

Systematic uncertainties
The systematic uncertainties are divided into two parts, which are additive and multiplicative in nature.The additive systematic uncertainty arises from the uncertainties in the signal parameterization and background modeling.The multiplicative systematic uncertainty includes contributions coming from MDC tracking, photon reconstruction, and selection criteria, which are related to the signal efficiency.The uncertainty due to signal and peaking background shapes is considered by performing the same fitting procedure with alternative signal and peaking background shapes.For the alternative signal shape, the parameters of the CB functions (except the mean value) and the ratio of the combinatorial part are obtained from the signal MC sample with the mass hypothesis m a closest to the search point.The alternative peaking background shape is modeled with the sum of a CB function and a Gaussian function.For each mass value m a , the fit is performed three times in total with different methods, and the minimum significance value and the maximum upper limit are recorded.The non-peaking background is described by polynomials.By performing the spurious signal test, the maximum of the absolute values of the spurious signal over the search ranges constitutes the related systematic uncertainty, which is the dominant source of sys-tematic uncertainty in the analysis and is incorporated in the overall likelihood assuming a Gaussian distribution: 8.09 events for 0.165 ≤ m a < 0.35 GeV/c 2 , 1.41 events for 0.35 ≤ m a < 0.75 GeV/c 2 , 17.95 events for 0.75 ≤ m a < 1.20 GeV/c 2 and 39.93 events for 1.20 ≤ m a ≤ 2.84 GeV/c 2 .
The multiplicative systematic uncertainties are listed in Table 2.The tracking efficiency of charged pions is investigated using control samples of J/ψ → ppπ + π − decays [50].The difference in tracking efficiencies between data and MC simulation is found to be 1.0% per track, which is taken as the uncertainty for the tracking efficiency.The photon detection efficiency is studied with a clean sample of J/ψ → ρ 0 π 0 decays [51].The difference in detection efficiencies between data and the MC simulation is 1.0% per photon.The systematic uncertainty due to the requirement on M (π + π − ) recoil is determined to be 0.3% according to the study of ψ(3686) → π + π − J/ψ, J/ψ → e + e − µ + µ − decays.The systematic uncertainties associated with the 4C kinematic fit, E other and N trk (number of charged tracks) requirements are studied with the control sample of ψ(3686) → π + π − J/ψ, J/ψ → γη decays and are determined to be 1.8%, 1.4% and 0.1%, respectively.The uncertainty on the ψ(3686) → π + π − J/ψ branching fraction is taken from the PDG [39], and the systematic uncertainty due to the number of ψ(3686) events is determined to be 0.6% by studying inclusive hadronic decays [34,35].The multiplicative systematic uncertainty is included in the overall likelihood as a Gaussian nuisance parameter with a width equal to the uncertainty.

Result
For each m a hypothesis, the local significance is determined by S = sign(N sig ) • 2 ln(L max /L 0 ), where L max and L 0 are the likelihood values with and without the signal hypothesis included in the fit, respectively.The fitted signal yields N sig and the corresponding significances are shown in Fig. 4. The largest local significance of 2.6σ is observed near m a = 2.208 GeV/c 2 , consistent with the null hypothesis.The fit results for m a = 2.208 GeV/c 2 are shown in Fig. 3(b).
Since no significant ALP signal is observed, we compute 95% confidence level (CL) upper limits on B(J/ψ → γa) as a function of m a using a one-sided frequentist profile-likelihood method [52].Then the branching fraction limit is converted to the ALP-photon coupling limit using [29] where B(J/ψ → e + e − ) = (5.971± 0.032)% is the world average value from the PDG [39] and α em is the electromagnetic coupling.An additional 0.5% uncertainty arising from the knowledge of B(J/ψ → e + e − ) is included when converting B(J/ψ → γa) to g aγγ .The expected and observed upper limits at 95% CL on B(J/ψ → γa) are shown in Fig. 5.The observed limits range from 8.3 × 10 −8 to 1.8 × 10 −6 in the ALP mass region of 0.165 ≤ m a ≤ 2.84 GeV/c 2 .The exclusion limits in the ALP-photon coupling g aγγ versus ALP mass m a plane obtained from this analysis are shown in Fig. 6, together with the constraints of other experiments.Our limits exclude the region in the ALP-photon coupling range g aγγ > 3 × 10 −4 GeV −1 for an ALP mass m a around 0.25 GeV/c 2 , with an improvement by a factor of 2-3 over the previous Belle II measurement [22].In addition, the constraints on the ALP-photon coupling are the most stringent to date for 0.165 ≤ m a ≤ 1.468 GeV/c 2 .6: Exclusion limits at 95% CL in the ALP-photon coupling gaγγ versus ALP mass ma plane obtained from this analysis.The constraints from BABAR [19], OPAL [20], Belle II [22], CMS [23], ATLAS [24], PrimEx [14], and beam-dump experiments [15,16] are also shown.All measurements assume a 100% ALP decay branching fraction into photons.

Summary
Based on a data sample of (2.71 ± 0.01) × 10 9 ψ(3686) events collected by the BESIII detector, we search for ALPs decaying into two photons produced in J/ψ radiative decays using the ψ(3686) → π + π − J/ψ process.No significant ALP signal is observed and we set 95% CL upper limits on the branching fraction of the decay J/ψ → γa and the ALP-photon coupling g aγγ .The observed limits on B(J/ψ → γa) range from 8.3 × 10 −8 to 1.8 × 10 −6 in the ALP mass region of 0.165 ≤ m a ≤ 2.84 GeV/c 2 , and the exclusion limits on the ALP-photon coupling are the most stringent to date for 0.165 ≤ m a ≤ 1.468 GeV/c 2 .

FIG. 3 :
FIG. 3: The Mγγ distribution for (a) ma = 0.327 GeV/c 2 and (b) ma = 2.208 GeV/c 2 with the fit results overlaid.The black dots with error bars are data, and the contribution of the non-peaking background is represented by a blue dashed curve.The green shaded region and the red solid curve represent the signal PDF and the total PDF, respectively.The inset of the figure (a) displays an enlargement of the Mγγ region between 0.25 and 0.4 GeV/c 2 .The largest local significance is 2.6σ at the ma = 2.208 GeV/c 2 hypothesis.

FIG. 4 :
FIG. 4:The distribution of the signal yields N sig (upper plot) and the local significance S (lower plot) obtained from the maximumlikelihood fits as a function of ma.The vertical gray bands indicate the excluded regions close to the π 0 , η, and η masses.

FIG. 5 :
FIG.5: Expected and observed upper limits at 95% CL on B(J/ψ → γa).The black curve is for the data, the black dashed curve represents the expected values and the green (yellow) band represents the ±1σ (±2σ) region.

TABLE 1 :
The Mγγ fit intervals for various ma points.

TABLE 2 :
The multiplicative systematic uncertainties.The total systematic uncertainty is obtained by adding all individual uncertainties in quadrature, assuming all sources to be independent.