Search for a $CP$-odd light Higgs boson in $J/\psi \to \gamma A^0$

Using $J/\psi$ radiative decays from 9.0 billion $J/\psi$ events collected by the BESIII detector, we search for di-muon decays of a $CP$-odd light Higgs boson ($A^0$), predicted by many new physics models beyond the Standard Model, including the Next-to-Minimal Supersymmetric Standard Model. No evidence for the $CP$-odd light Higgs production is found, and we set $90\%$ confidence level upper limits on the product branching fraction $\mathcal{B}(J/\psi \to \gamma A^0)\times \mathcal{B}(A^0 \to \mu^+\mu^-)$ in the range of $(1.2-778.0)\times 10^{-9}$ for $0.212 \le m_{A^0} \le 3.0$ GeV/$c^2$. The new measurement is a 6-7 times improvement over our previous measurement, and is also slightly better than the BaBar measurement in the low-mass region for $\tan \beta =1$.

Using J/ψ radiative decays from 9.0 billion J/ψ events collected by the BESIII detector, we search for di-muon decays of a CP -odd light Higgs boson (A 0 ), predicted by many new physics models beyond the Standard Model, including the Next-to-Minimal Supersymmetric Standard Model. No evidence for the CP -odd light Higgs production is found, and we set 90% confidence level upper limits on the product branching fraction B(J/ψ → γA 0 ) × B(A 0 → µ + µ − ) in the range of (1.2 − 778.0) × 10 −9 for 0.212 ≤ m A 0 ≤ 3.0 GeV/c 2 . The new measurement is a 6-7 times improvement over our previous measurement, and is also slightly better than the BaBar measurement in the low-mass region for tan β = 1.
The origin of mass is one of the most important questions in physics. The masses of the fundamental particles are generated through spontaneous breaking of electroweak symmetry by the Higgs mechanism [1]. The Higgs mechanism implies the existence of at least one new scalar particle, the Higgs boson, which was the last missing Standard Model (SM) particle. It was discovered by the Large Hadron Collider experiments at CERN [2] in July 2012 and has a profound effect on our fundamental understanding of matter.
Many models beyond the SM, such as the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [3][4][5], extend the Higgs sector to include additional Higgs fields. The NMSSM adds an additional singlet chiral superfield to the Minimal Supersymmetric Standard Model (MSSM) [6] to alleviate the so-called "little hierarchy problem" [7]. It contains three CP -even, two CP -odd, and two charged Higgs bosons [3,4]. The mass of the lightest Higgs boson, A 0 , may be smaller than twice the mass of the charmed quark, thus making it accessible via J/ψ → γA 0 [8].
The branching fraction of V → γA 0 (V = Υ, J/ψ) is expressed as [8][9][10] where α is the fine structure constant, G F is the Fermi coupling constant, l ≡ e or µ, m q is the quark mass, C QCD includes the leptonic width of B(V → l + l − ) [11,12] as well as m A 0 dependent QCD and relativistic corrections to B(V → γA 0 ) [10], and g q is the effective Yukawa coupling to the Higgs field to the up-or down-type quarkpair. In the NMSSM, g c = cos θ A / tan β (q = c) for the charm quark and g b = cos θ A tan β (q = b) for the bottom quark, where tan β is the ratio of up-and downtype Higgs doublets, and cos θ A is the fraction of the nonsinglet component of the A 0 [13,14]. The value of cos θ A is zero for a pure A 0 singlet state [15]. The branching fraction of J/ψ → γA 0 is predicted to be in the range of 10 −9 − 10 −7 depending upon the A 0 mass and the NMSSM parameters [4]. The branching fraction of A 0 → µ + µ − is predicted to be much larger for tan β ≥ 1 [13]. An experimental study of such a lowmass Higgs boson is desirable to test the SM [16] and to look for new physics beyond the SM [3,4,17]. The BaBar [18], CLEO [19], and CMS [20] experiment have searched for di-muon decays of A 0 , and placed a strong exclusion upper limit on g b . On the other hand, the BESIII measurements, sensitive on g c , is complementary to those by considering g b . The recent BESIII measurement [21], based on 225 million J/ψ events, is slightly lower than the BaBar measurement [18] in the low-mass region for tan β ≤ 0.6. The combined measurements of the BESIII and BaBar have revealed that the A 0 is mostly singlet in nature because of obtained upper limit on cos tan β, is very close to zero especially in the low-mass region [21]. However, this BESIII limit [21] is still an order of magnitude above the theoretical predictions [4]. BE-SIII has recently accumulated about 39 times more J/ψ events in comparison to the previous measurement [21], and these can be utilized to discover the A 0 or exclude parameter space of the NMSSM [22]. This paper describes the search for di-muon decays of a CP -odd light Higgs boson in radiative decays of J/ψ using 9 billion J/ψ events collected by the BESIII detector in 2009, 2018, and 2019 [22]. Because muon particle identification (PID) was not available for the J/ψ data collected in 2012, we exclude this data sample for the A 0 search.

I. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [23] records symmetric e + e − collisions provided by the BEPCII storage ring [24], which operates with a peak luminosity of 1 × 10 33 cm −2 s −1 in the center-of-mass (CM) energy range from 2.0 to 4.95 GeV. BESIII has collected large data samples in this energy region [25]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a heliumbased 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 MDC measures the momentum of charged particles with a resolution of 0.5% at 1 GeV/c. The EMC measures the photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region. The time resolution of the TOF in the barrel region is 68 ps. The time resolution of the TOF in the endcap region was 110 ps before 2015 and was improved to be 60 ps after upgrading with the multi-gap resistive plate chambers. Muons with momentum above 0.5 GeV/c are identified by the iron flux return of the magnet instrumented with about 1272 m 2 of resistive plate muon counters (MUC) arranged in nine (eight) layers in the barrel (endcaps).
Simulated Monte Carlo (MC) events based on Geant4 [26] are used to optimize the event selection criteria, to study the potential backgrounds, and to determine the detector acceptance. A MC sample of 9.0 billion inclusive J/ψ events is used for the background studies with the generic TopoAna tool [27]. The known J/ψ decay modes are generated by the EvtGen generator [28] with branching fractions taken from the Particle Data Group (PDG) [29], and the remaining unknown decay modes by LUNDCHARM [30]. The final state radiation corrections are included in the MC simulation using PHOTOS [31]. The production of the J/ψ resonance through e + e − annihilation including the beam-energy spread and the initial-state-radiation (ISR) is simulated by the KKMC [32]. A 2.93 fb −1 ψ(3770) data sample [33,34] is used to study the background from the quantum electrodynamics (QED) process of e + e − → γµ + µ − . To compute the detection efficiency, we generate 0.12 million simulated signal MC events at 23 different Higgs mass points ranging from 0.212 to 3.0 GeV/c 2 with a phase-space model for the A 0 → µ + µ − decay and a P -wave model for the J/ψ → γA 0 decay [28].

II. DATA ANALYSIS
We select events with two oppositely charged tracks and at least one photon candidate. A photon candidate, reconstructed with clusters of energy deposited in the EMC, is selected with a minimum energy of 25 MeV in the barrel region (| cos θ| < 0.8) or 50 MeV in the endcap region (0.86 < | cos θ| < 0.92). The energy deposited in the nearby TOF is included to improve the energy resolution and reconstruction efficiency. The angle between a photon and the nearest extrapolated track in the EMC is required to be larger than 10 degrees to remove bremsstrahlung photons. The EMC timing is required to be within 700 ns relative to the event start time to suppress electronic noise and energy deposits unrelated to the signal events.
Charged tracks are reconstructed from the ionization signals measured by the MDC and are required to be in the MDC detection acceptance region of | cos θ| < 0.93, where θ is the angle of the charged track with the z axis, which is the axis of the MDC. Further, their points of closest approach to the z-axis must be within ±10 cm from the interaction point along the z direction and within ±1 cm in the plane perpendicular to z. To suppress contamination by electrons and pions, both charged tracks are required to satisfy the following selection criteria: 1) E µ cal /p < 0.9 c, 2) 0.1 < E µ cal < 0.3 GeV, and 3) the absolute value of the time difference between the TOF and expected muon time (∆t TOF ) must be less than 0.26 ns. Here, E µ cal is the energy deposited in the EMC by the µ + /µ − particle, and p is the momentum of the charged muon track. To further improve the purity of muons, one of the charged tracks is required to have its penetration depth in the MUC be greater than (−40.0 + 70 × p/(GeV/c)) cm for 0.5 ≤ p ≤ 1.1 GeV/c and 40 cm for p > 1.1 GeV/c.
The two muon tracks are required to originate from a common vertex by performing a vertex fit to form an A 0 candidate. A four-constraint (4C) kinematic fit is performed with the two charged tracks and one of the photon candidates in order to improve the mass resolution of the A 0 candidate. If there is more than one γµ + µ − candidate, the candidate with the minimum value of the χ 2 from the 4C kinematic fit (χ 2 4C ) is selected, and the χ 2 4C is required to be less than 40 to reject backgrounds from J/ψ → π + π − π 0 and e + e − → γπ + π − π 0 . We reject fake photons by requiring the di-muon invariant mass (m µ + µ − ) obtained from the 4C kinematic fit to be less than 3.04 GeV/c 2 . To suppress backgrounds from e + e − → γµ + µ − and J/ψ → µ + µ − (γ), the absolute value of the cosine of the muon helicity angle (cos θ hel µ ), defined as the angle between the direction of one of the muons and the direction of the J/ψ in the A 0 rest frame, is required to be less than 0.92.
After the above selection criteria, we determine the signal yield as a function of m A 0 in the interval of 0.212 ≤ m A 0 ≤ 3.0 GeV/c 2 by performing a series of one-dimensional unbinned extended maximum likelihood (ML) fits to the reduced mass, m red = m 2 µ + µ − − 4m 2 µ distribution of surviving events. Fig. 1 shows the m red distribution of surviving events together with the background predictions from various simulated MC samples and 2.93 fb −1 of ψ(3770) data [33,34]. We use m red rather than m µ + µ − because it is easier to model the nonpeaking background across the entire m A 0 region, in particular, the kinematic threshold region m µ + µ − ≈ 2m µ (m red ≈ 0). The non-peaking background is dominated by e + e − → γµ + µ − and J/ψ → µ + µ − (γ), and the peaking background by J/ψ → ρ/ωπ and J/ψ → γf (f = f 2 (1270), f 0 (1500), f 0 (1710)) decays, where both ρ/ω and f decay to π + π − . The m red distribution of data is generally well described by the background predictions, except in the low-mass region, where KKMC [32] fails to reproduce the ISR events for the e + e − → γJ/ψ, J/ψ → µ + µ − process. This disagreement has little impact on the search because the signal extraction procedure does not depend on the background predictions. The fit function includes the contributions of signal, continuum background and peaking background components from ρ/ω, f 2 (1270), f 0 (1500), and f 0 (1710) mesons. Table I summarizes the m red fit intervals for the various m A 0 points used to handle both non-peaking and peaking backgrounds smoothly.
Simulated MC samples are used to develop the probability density functions (PDFs) of signal and backgrounds. The A 0 is assumed to be a scalar or pseudoscalar particle with a very narrow decay width in compar-  ison to the experimental resolution [35]. We describe the m red distribution of the signal PDF by the sum of two Crystal Ball (CB) functions [36] with a common peak value and opposite side tails. The m red resolution varies from 2 MeV/c 2 to 12 MeV/c 2 while the signal efficiency varies between 27% and 53% depending upon the muon momentum values at different A 0 mass points. We interpolate the signal efficiency and signal PDF parameters linearly between the mass points of the generated signal MC events. The non-peaking background PDF is described by a function tanh( 5 l=1 p l m l red ) in the threshold mass region of 0.212 ≤ m A 0 ≤ 0.40 GeV/c 2 , where p l are the polynomial coefficients. This function provides a threshold like behavior in the low-mass region of the m red distribution and passes through the origin when m red = 0. In the other mass regions, we use second, third, fourth, fifth or sixth-order Chebyshev polynomial function to describe the non-peaking background PDFs detailed in Table I. We determine the initial parameters of these background PDFs using a cocktail MC sample of all possible non-peaking backgrounds to achieve better agreement between data and the fit models.
To take into account the well-known structure of the ρω interference, we describe the peaking background PDF of the m red distribution with the Gounaris and Sakurai (GS) function in the range of 0.4 ≤ m red ≤ 1.06 GeV/c 2 [37]. The fit formula, detailed in Ref. [38], is the same as that used previously by the BaBar [38] and BESIII [34] experiments in the measurement of the e + e − → π + π − cross-section in the ρ/ω mass region. The amplitudes for the higher ρ states, ρ(1450), ρ(1700), and ρ(2150), as well as the massses and widths of those states are taken from Ref. [38]. We fix the ω width according to the PDG [29] value and float the other parameters during the fit. We describe the peaking background PDFs corresponding to f 2 (1270), f 0 (1500) and f 0 (1710) resonances by the sum of the two CB functions [36] using the parameters determined from MC samples of J/ψ → γf , f → π + π − decays, where f = f 2 (1270), f 0 (1500), and f 0 (1710) mesons.
The search for the A 0 narrow resonance is performed in steps of approximately half the m red resolution, i.e., 1 MeV/c 2 in the mass range of 0.22 ≤ m A 0 ≤ 1.5 GeV/c 2 and 2.0 MeV/c 2 in the other m A 0 regions, with a total of 2,035 m A 0 points. The PDF parameters of the signal and peaking backgrounds of J/ψ → γf are fixed while the non-peaking background PDF, and the numbers of the signal, peaking, and non-peaking background events are floated. Plots of the fit to the m red distribution for two selected mass points are shown in Fig. 2.
The fit is repeated with alternative signal, peaking, and non-peaking background PDFs to determine the systematic uncertainties for the numbers of signal events associated with the corresponding PDFs at each m A 0 point. The uncertainty associated with the signal PDF is studied by replacing the sum of the two CB functions with a 'Cruijff' function [39]. The uncertainty associated with the ρ − ω peak is evaluated by varying the ρ and ω contributions in the formula of Eq.(26) of Ref. [38]. The uncertainty due to the peaking background of J/ψ → γf is studied by replacing the sum of the two CB functions with the simulated MC samples of the corresponding decay processes convolved with a Gaussian function whose parameters are floated during the fit. The uncertainty due to the non-peaking background PDF is studied by replacing the tanh( 5 l=1 p l m l red ) and n th order Chebyshev polynomial function with tanh( 6 l=1 p l m l red ) and (n + 1) th order Chebyshev polynomial functions, respectively, in the fit. The one with the largest signal yield among these fit scenarios is considered to produce the final result.
The product branching fraction of J/ψ → γA 0 and A 0 → µ + µ − as a function of m A 0 is calculated as where N sig is the number of signal events, is the sig- The corresponding local significance values at these mass points are observed to be 3.3σ and 3.5σ, respectively. Black dots with error bars represent the data, the red long-dashed curve the non-peaking background, the pink dotted curve the peaking background, the green dashed curve the signal PDF, and the solid blue curve the total fit results. In the bottom figure, the wellknown ρ − ω interference is taken care of by describing the peaking background PDF of the m red distribution by a GS function [37,38], as described in the text. nal selection efficiency, and N J/ψ = (8.998 ± 0.039) × 10 9 is the number of J/ψ events. Fig. 3 shows the plots of the product branching fractions B(J/ψ → γA 0 ) × B(A 0 → µ + µ − ) and the statistical significance, defined as S = sign(N sig ) 2ln(L max /L 0 ), where L max (L 0 ) is the maximum likelihood value for a fit with number of signal events being floated (fixed at zero). The largest upward local significance value is determined to be 3.5σ at m A 0 = 0.696 GeV/c 2 . Based on a large ensemble of pseudo experiments [18], the probability of observing a fluctuation of S ≥ 3.5σ is estimated to be 12%. The corresponding global significance value is determined to be at the level of 1σ. Thus, we conclude that no evidence of Higgs production is found within the searched m A 0 regions.

III. SYSTEMATIC UNCERTAINTIES
According to Eq. 2, the systematic uncertainties for the branching fraction measurement include those from the number of signal events, the reconstruction efficiency, and the number of J/ψ events. The uncertainties associated with the number of signal events originating from the PDF parameters of signal and backgrounds are considered by performing alternative fits at each m A 0 point.
Pseudo experiments are utilized to test the reliability of the fit procedures and compute the fit bias, which may appear due to imperfect signal and background modeling. The same fit procedure is performed in each pseudo experiment. The resultant average difference between the input and output signal yields is determined to be 0.3 events. We consider it as an additive systematic uncertainty (σ add ), which may affect the significance of any observation but does not scale with the reconstructed signal yield.
The uncertainties associated with the reconstruction efficiency and the number of J/ψ events don't affect the significance of any observation. Thus, we consider them as multiplicative systematic uncertainties (σ mult ) and scale with the number of reconstructed signal events. The uncertainty associated with the reconstruction efficiency includes those from tracking, PID, and the photon The uncertainty due to MDC tracking is determined to be 1% per track using the high statistics control samples of J/ψ → ρπ and J/ψ → ppπ + π − . A total of 2.0% systematic uncertainty is assigned for the two charged tracks in this analysis. The systematic uncertainty associated with the photon reconstruction efficiency is determined using a control sample of e + e − → γµ + µ − in which the ISR photon is predicted using the four momenta of the two charged tracks. This sample also includes the dominant contribution from J/ψ → γπ + π − decay, including all the possible intermediate resonances. The relative difference in efficiency between data and MC is found to be 0.2%, which is considered as a systematic uncertainty.
A control sample of J/ψ → µ + µ − (γ) is used to evaluate the systematic uncertainty due to the muon PID, cos θ hel µ , and χ 2 4C requirements. In this sample, one track is tagged with a tight muon PID. The final uncertainty associated with the muon PID also takes into account the fraction of events with one or two tracks identified as muons obtained from the simulated signal MC sample. The corresponding uncertainties, computed as the relative change in efficiency between data and MC, are determined to be (2.9 − 4.1)%, 0.8% and 1.8%, respectively. The systematic uncertainty due to the number of J/ψ events is 0.3% using J/ψ inclusive hadronic events. Table II summarizes the fit bias and multiplicative sources of the systematic uncertainties, where we obtain the total σ mult by adding the individual ones in quadrature. The total σ mult varies between 4.1% to 5.0% depending on the Higgs mass point. The final systematic uncertainty is calculated as σ 2 add + (σ mult × N sig ) 2 . The inner and outer bands correspond to 68% and 95% of the expected limit values and include the statistical uncertainties only.

IV. RESULT
Since no evidence of Higgs production is found, we set 90% confidence level (C.L.) upper limits on the product branching fractions B(J/ψ → γA 0 ) × B(A 0 → µ + µ − ) as a function of m A 0 using a Bayesian method [29] after incorporating the systematic uncertainty by smearing the likelihood curve with a Gaussian function having a width equal to the systematic uncertainty. The limits vary in the range of (1.2 − 778.0) × 10 −9 for the Higgs mass region of 0.212 ≤ m A 0 ≤ 3.0 GeV/c 2 depending on the m A 0 point, as shown in Fig. 4. The new measurement has a 6-7 times improvement over the previous BESIII measurement [21].
To compare our results with the BaBar measurement [18], we also compute 90% C.L. upper limits on the effective Yukawa coupling of the Higgs fields to the bottom-quark pair g b (= g c tan 2 β)× B(A 0 → µ + µ − ) as a function of m A 0 for different values of tan β using Eq.(1) as shown in Fig. 5. Our new measurement is slightly better than the BaBar measurement [18] in the low-mass region for tan β = 1.0.

V. SUMMARY
We search for di-muon decays of A 0 in J/ψ → γA 0 using 9.0 billion J/ψ events collected by the BESIII detector. No evidence of Higgs production is found, and we set 90% C.L. upper limits on product branching fractions B(J/ψ → γA 0 ) × B(A 0 → µ + µ − ) in the range of