Search for an axion-like particle in radiative $J/\psi$ decays

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\gamma\gamma}$. As a generalization of QCD axions, ALPs have arbitrary masses and couplings which are bounded by experiments. In the sub-$\mbox{MeV}/c^{2}$ mass range, the bounds on $g_{a\gamma\gamma}$ are provided by laser experiments, solar photon instruments, as well as cosmological and astrophysical observations [12, 13]. In the $\mbox{MeV}/c^{2}$ to $\mbox{GeV}/c^{2}$ mass region, constraints on the photon coupling constant mainly come from beam-dump experiments [14, 15, 16] and high-energy collider experiments. Long-lived ALPs have been searched for in $J/\psi$ and $\Upsilon$ radiative decays [17, 18] and through $e^{+}e^{-}\rightarrow\gamma+\text{invisible}$ processes [19]. The limits for short-lived ALPs decaying to two photons are obtained from $e^{+}e^{-}\rightarrow\gamma\gamma$ [20], $e^{+}e^{-}\rightarrow\gamma\gamma\gamma$ [21, 22] via ALP-strahlung production, and the light-by-light scattering process $\gamma\gamma\rightarrow\gamma\gamma$ [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\,\mbox{MeV}/c^{2}$ to roughly $10\,\mbox{GeV}/c^{2}$, leaving the large area $10^{-3}\lesssim g_{a\gamma\gamma}\lesssim 10^{-5}\,\mbox{GeV}^{-1}$ still uncovered. 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/\psi$ radiative decays via $J/\psi\rightarrow\gamma a$, $a\rightarrow\gamma\gamma$ in the mass range of $0.165\leq m_{a}\leq 2.84\,\mbox{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 $\Gamma_{a}=g_{a\gamma\gamma}^{2}m_{a}^{3}/64\pi$. 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\,\mbox{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\,\mbox{GeV}/c^{2}$, $0.46<m_{a}<0.60\,\mbox{GeV}/c^{2}$ and $0.90<m_{a}<1.01\,\mbox{GeV}/c^{2}$ are excluded from the search due to the peaking backgrounds from the decays of $\pi^{0}$, $\eta$ and $\eta^{\prime}$ mesons, respectively. In order to avoid pollution from non-resonant ALP production $e^{+}e^{-}\rightarrow\gamma a$ and to probe resonant ALP production $J/\psi\rightarrow\gamma a$ in a model-independent way [29], the decay $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ is exploited to isolate the $J/\psi$ sample.

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\times 10^{33}$ cm^-2s^-1 achieved at $\sqrt{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\leavevmode\nobreak\ {\rm GeV}/c$ is $0.5\%$, and the ${\rm d}E/{\rm d}x$ 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 $\psi(3686)$ data samples collected by the BESIII detector during three data-taking periods: 2009, 2012, and 2021. The numbers of $\psi(3686)$ events are determined by counting inclusive hadronic events whose branching fraction is known rather precisely [34, 35]. There are $(107.0\pm 0.8)\times 10^{6}$ $\psi(3686)$ events collected in 2009, $(341.1\pm 2.9)\times 10^{6}$ events collected in 2012 and $(2.26\pm 0.01)\times 10^{9}$ events collected in 2021, for a total of $(2.71\pm 0.01)\times 10^{9}$ $\psi(3686)$ events. In this analysis, the $J/\psi$ sample originates from the decay $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$. The analysis strategy is to first tag $J/\psi$ 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/\psi$ 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 geant4-based [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 initial-state radiation in $e^{+}e^{-}$ annihilations with the generator kkmc [37]. An inclusive MC sample with approximately the same number of $\psi(3686)$ decays as in data is used to check for potential backgrounds. The inclusive MC sample includes the production of the $\psi(3686)$ resonance, the initial-state radiation production of the $J/\psi$, and the continuum processes incorporated in kkmc [37]. All particle decays are modeled with evtgen [38] using branching fractions either taken from the Particle Data Group (PDG) [39], when available, or otherwise estimated with lundcharm [40]. Final-state radiation from charged final-state particles is incorporated using the photos package [41]. To simulate the signal, the decay $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ is modeled according to the partial wave analysis results of $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ [42]; a $P$-wave model [38] is used for the subsequent decay $J/\psi\rightarrow\gamma a$ and a phase-space model is used for the $a\rightarrow\gamma\gamma$ decay. The signal events are simulated individually for different mass hypotheses of ALP $m_{a}$ in the range of $0.165\leq m_{a}\leq 2.84\,\mbox{GeV}/c^{2}$ with a step size of $5\,\mbox{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 $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ with subsequent decays of $J/\psi\rightarrow\gamma\pi^{0}$, $\gamma\eta$ and $\gamma\eta^{\prime}$. These backgrounds are each generated exclusively with the angular distribution of $1+\cos^{2}\theta_{\gamma}$ [38], where $\theta_{\gamma}$ is the angle of the radiative photon relative to the positron beam direction in the $J/\psi$ rest frame. Some potential backgrounds of the form $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ with $J/\psi$ decaying into purely neutral particles in the final states are generated exclusively with different generators: $J/\psi\rightarrow\gamma\eta_{c}$ with the angular distribution of $1+\cos^{2}\theta_{\gamma}$ [38], $J/\psi\rightarrow\gamma\pi^{0}\pi^{0}$ according to the partial wave analysis results of $J/\psi\rightarrow\gamma\pi^{0}\pi^{0}$ [44], as well as $J/\psi\rightarrow\gamma\gamma\gamma$ with a phase-space distribution. Each background MC sample is weighted to match the integrated luminosity of the data set.

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 ($\theta$) range of $|\rm{cos\theta}|<0.93$, where $\theta$ 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\,\mbox{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(\pi^{+}\pi^{-})_{\rm recoil}$ must be in the range of [3.080, 3.114] $\,\mbox{GeV}/c^{2}$.

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\theta|<0.80$) and more than 50 MeV in the end-cap region ($0.86<|\!\cos\theta|<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 $\pi^{+}$ and $\pi^{-}$ 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 $\psi(3686)$ is applied for all possible $\pi^{+}\pi^{-}\gamma\gamma\gamma$ combinations. The combination with the smallest fit $\chi_{\rm 4C}^{2}$ is retained, with the requirement of $\chi_{\rm 4C}^{2}<40$. The distribution of $\chi_{\rm 4C}^{2}$ is shown in Fig. 1. To further reject background from $J/\psi\rightarrow\gamma\pi^{0}\pi^{0}$ decays, we form all $\pi^{+}\pi^{-}5\gamma$ combinations if 5 or more photons are present, and we require the minimum $\chi^{2}_{\rm 4C}(5\gamma)$ of the corresponding kinematic fits to satisfy $\chi^{2}_{\rm 4C}<\chi^{2}_{\rm 4C}(5\gamma)$. For the remaining photons, the sum of their energy deposition in the EMC $E_{\rm other}$ is required to be less than $0.1\,\mbox{GeV}$. The requirements of $\chi_{\rm 4C}^{2}$ and $E_{\rm other}$ are optimized according to the Punzi significance [45] defined as $\epsilon/(1.5+\sqrt{B})$, where $\epsilon$ denotes the signal efficiency obtained from signal MC samples and $B$ is the number of background events obtained from background MC samples.

The two-photon invariant mass $M_{\gamma\gamma}$ 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/\psi\rightarrow\gamma\pi^{0}$, $\gamma\eta$ and $\gamma\eta^{\prime}$ decays. There is a small contribution from residual $J/\psi\rightarrow\gamma\pi^{0}\pi^{0}$ background with two soft photons peaking around the $\pi^{0}$ mass region, and a contamination from $J/\psi\rightarrow\gamma f_{2}(1270)$ decays, followed by $f_{2}(1270)\rightarrow\pi^{0}\pi^{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_{\gamma\gamma}$ distribution to determine the signal yields with different ALP mass hypotheses in the mass range of $0.165\leq m_{a}\leq 2.84\,\mbox{GeV}/c^{2}$. The search step is $3\,\mbox{MeV}/c^{2}$ for $0.165\leq m_{a}<1.20\,\mbox{GeV}/c^{2}$ and $4\,\mbox{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 $\pi^{0}$, $\eta$ and $\eta^{\prime}$ mesons. Different $M_{\gamma\gamma}$ 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.

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\,\mbox{MeV}/c^{2}$ near $m_{a}=0.165\,\mbox{GeV}/c^{2}$ to $11\,\mbox{MeV}/c^{2}$ near $m_{a}=2.2\,\mbox{GeV}/c^{2}$, and it decreases back to $7\,\mbox{MeV}/c^{2}$ near $m_{a}=2.84\,\mbox{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 $\pi^{0}$, $\eta$ and $\eta^{\prime}$ 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\leq m_{a}<1.20\,\mbox{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 signal-plus-background fit to the background-only $M_{\gamma\gamma}$ 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_{\gamma\gamma}$ distributions for selected $m_{a}$ points and the fit results are shown in Fig. 3.

The branching fractions of $J/\psi\rightarrow\gamma(\pi^{0},\eta,\eta^{\prime})\rightarrow\gamma\gamma\gamma$ 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.

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 systematic uncertainty in the analysis and is incorporated in the overall likelihood assuming a Gaussian distribution: 8.09 events for $0.165\leq m_{a}<0.35\,\mbox{GeV}/c^{2}$, 1.41 events for $0.35\leq m_{a}<0.75\,\mbox{GeV}/c^{2}$, 17.95 events for $0.75\leq m_{a}<1.20\,\mbox{GeV}/c^{2}$ and 39.93 events for $1.20\leq m_{a}\leq 2.84\,\mbox{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/\psi\rightarrow p\bar{p}\pi^{+}\pi^{-}$ 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/\psi\rightarrow\rho^{0}\pi^{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(\pi^{+}\pi^{-})_{\rm recoil}$ is determined to be 0.3% according to the study of $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$, $J/\psi\rightarrow e^{+}e^{-}\mu^{+}\mu^{-}$ decays. The systematic uncertainties associated with the 4C kinematic fit, $E_{\rm other}$ and $N_{\rm trk}$ (number of charged tracks) requirements are studied with the control sample of $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$, $J/\psi\rightarrow\gamma\eta$ decays and are determined to be 1.8%, 1.4% and 0.1%, respectively. The uncertainty on the $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ branching fraction is taken from the PDG [39], and the systematic uncertainty due to the number of $\psi(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.

For each $m_{a}$ hypothesis, the local significance is determined by $\mathcal{S}=\mathrm{sign}(N_{\mathrm{sig}})\cdot\sqrt{2\ln(\mathcal{L}_{\mathrm{max}}/\mathcal{L}_{0})}$, where $\mathcal{L}_{\mathrm{max}}$ and $\mathcal{L}_{0}$ are the likelihood values with and without the signal hypothesis included in the fit, respectively. The fitted signal yields $N_{\rm sig}$ and the corresponding significances are shown in Fig. 4. The largest local significance of $2.6\sigma$ is observed near $m_{a}=2.208\,\mbox{GeV}/c^{2}$, consistent with the null hypothesis. The fit results for $m_{a}=2.208\,\mbox{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 $\mathcal{B}(J/\psi\rightarrow\gamma 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 $\mathcal{B}(J/\psi\rightarrow e^{+}e^{-})=(5.971\pm 0.032)\%$ is the world average value from the PDG [39] and $\alpha_{\mathrm{em}}$ is the electromagnetic coupling. An additional 0.5% uncertainty arising from the knowledge of $\mathcal{B}(J/\psi\rightarrow e^{+}e^{-})$ is included when converting $\mathcal{B}(J/\psi\rightarrow\gamma a)$ to $g_{a\gamma\gamma}$.

The expected and observed upper limits at 95% CL on $\mathcal{B}(J/\psi\rightarrow\gamma a)$ are shown in Fig. 5. The observed limits range from $8.3\times 10^{-8}$ to $1.8\times 10^{-6}$ in the ALP mass region of $0.165\leq m_{a}\leq 2.84\,\mbox{GeV}/c^{2}$. The exclusion limits in the ALP-photon coupling $g_{a\gamma\gamma}$ 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\gamma\gamma}>3\times 10^{-4}\,\mbox{GeV}^{-1}$ for an ALP mass $m_{a}$ around $0.25\,\mbox{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\leq m_{a}\leq 1.468\,\mbox{GeV}/c^{2}$.

Based on a data sample of $(2.71\pm 0.01)\times 10^{9}$ $\psi(3686)$ events collected by the BESIII detector, we search for ALPs decaying into two photons produced in $J/\psi$ radiative decays using the $\psi(3686)\rightarrow\pi^{+}\pi^{-}J/\psi$ process. No significant ALP signal is observed and we set 95% CL upper limits on the branching fraction of the decay $J/\psi\rightarrow\gamma a$ and the ALP-photon coupling $g_{a\gamma\gamma}$. The observed limits on $\mathcal{B}(J/\psi\rightarrow\gamma a)$ range from $8.3\times 10^{-8}$ to $1.8\times 10^{-6}$ in the ALP mass region of $0.165\leq m_{a}\leq 2.84\,\mbox{GeV}/c^{2}$, and the exclusion limits on the ALP-photon coupling are the most stringent to date for $0.165\leq m_{a}\leq 1.468\,\mbox{GeV}/c^{2}$.