Search for an Axion-Like Particle in B Meson Decays

Axion-like particles (ALPs) are predicted in many extensions of the Standard Model, and their masses can naturally be well below the electroweak scale. In the presence of couplings to electroweak bosons, these particles could be emitted in flavor-changing $B$ meson decays. We report herein a search for an ALP, $a$, in the reaction $B^\pm\rightarrow K^\pm a$, $a\rightarrow\gamma\gamma$ using data collected by the BABAR experiment at SLAC. No significant signal is observed, and 90% confidence level upper limits on the ALP coupling to electroweak bosons are derived as a function of ALP mass, improving current constraints by several orders of magnitude in the range $0.175\,\,\mathrm{GeV}

The physics of spontaneous symmetry breaking drives much of the phenomenology of the Standard Model (SM).For instance, the Higgs mechanism gives mass to the fermions and weak gauge bosons of the SM, while the spontaneous breaking of approximate chiral global symmetries gives rise to pseudo-Goldstone bosons, such as the pions.Many extensions of the SM feature anomalous global symmetries whose spontaneous breaking leads to new pseudo-Goldstone bosons known as axion-like particles (ALPs) [1][2][3][4].Such particles are ubiquitous in beyond-the-SM theories, such as supersymmetry [5][6][7], as well as in string theory [8][9][10][11].Potentially, ALPs could resolve several outstanding issues related to the naturalness of SM parameters, such as the strong CP problem [1][2][3][4] or the hierarchy problem [12], and they may also serve as mediators to dark sectors [13][14][15][16].Consequently, ALPs have motivated a large number of searches in experimental particle physics and cosmology [17][18][19][20].
In the simplest models, ALPs predominantly couple to pairs of SM gauge bosons.While the photon and gluon couplings are already significantly constrained by collider and beam-dump experiments for ALP masses in the MeV-GeV range [21][22][23][24][25][26][27][28], the coupling to W ± bosons is less explored.This coupling leads to ALP production in flavor-changing neutral-current decays, which can serve as powerful discovery modes.For example, flavorchanging B meson and kaon decays already provide the most stringent bounds on invisibly decaying ALPs over a range of masses [29].The search presented here is the first for visibly decaying ALPs produced in B meson decays.Its sensitivity complements existing studies of K → πγγ [30][31][32], which have been conservatively reinterpreted to obtain limits on ALP couplings [29].
In the following, we consider a minimal ALP (a) model with coupling g aW to the SU(2) W gauge-boson field strengths, W b µν , and Lagrangian where W bµν is the dual field-strength tensor.This coupling leads to the production of ALPs at one loop in the process B ± → K ± a, where the ALP is emitted from an internal W ± boson [29].Electroweak symmetry breaking and the resulting gauge-boson mixing generates an ALP coupling to a pair of photons, and the branching fraction for a → γγ in this model is nearly 100% for m a < m W .
The same ALP production and decay modes also occur in models with axion couplings to gluons [33].
We report herein the first search for an ALP in the reaction B ± → K ± a, a → γγ in the range 0.175 GeV < m a < m B + − m K + ≈ 4.78 GeV, excluding the mass intervals 0.45-0.63GeV and 0.91-1.01GeV because of large peaking backgrounds from η and η mesons, respectively.Note that existing searches already constrain m a < 0.1 GeV in the range of couplings to which our search is sensitive [22][23][24][25][26], while the mass range 0.1 GeV < m a < 0.175 GeV is excluded from our analysis due to large peaking π 0 contributions.The B ± → K ± a, a → γγ product branching fraction is measured assuming all signal observed is produced in B ± → K ± a with a decaying promptly.However, the ALP has a decay width Γ a = g 2 aW m 3 a sin 4 θ W /64π, where θ W is the weak mixing angle, and the present search has sensitivity to couplings predicting long-lived ALPs for m a < 2.5 GeV.We additionally determine the branching fraction for cτ a values of 1, 10, and 100 mm in this mass range.
The search is based on a sample of 4.72×10 8 B B meson pairs corresponding to 424 fb −1 of integrated luminosity collected at the Υ (4S) resonance by the BABAR detector at the PEP-II e + e − storage ring at the SLAC National Accelerator Laboratory [34].The BABAR detector is described in detail elsewhere [35,36].A small sample, corresponding to 8% of the total data set, is used to optimize the search strategy and is subsequently discarded.
We reconstruct signal B ± candidates by combining a pair of photons with a track identified as a kaon by particle identification algorithms [35].All other reconstructed tracks and neutral clusters in the event are collectively referred to as the rest of the event (ROE).To suppress backgrounds, we require an energy-substituted is the B ± energy in the CM frame, and E i and p i are the energy and momentum of the initial state in the lab frame.A kinematic fit is performed on the selected B ± candidates, requiring the photon and kaon candidates to originate from the measured beam interaction region, and constraining their total energy and invariant mass to the nominal B ± meson mass and measured CM beam energy.
Two boosted decision trees (BDTs) [43] are used to further separate signal from each of the main backgrounds: one BDT is trained using continuum MC background events and the other using B + B − MC background events.For the signal sample, we combine events from all simulated ALP masses with prompt decays to obtain a uniform distribution in diphoton invariant mass (m γγ ).Each BDT includes the following 13 observables: invariant mass of the ROE; cosine of the angle between two sphericity axes, one computed with the B ± constituents and the other with the ROE; second Legendre moment of the ROE, calculated relative to the B ± thrust axis; m ES and ∆E; particle identification information for the K ± ; helicity angle of the K ± , which is the angle between the K ± and the Υ (4S) as measured in the B ± frame; helicity angle and energy of the most energetic photon forming the a; three invariant masses m(γ i γ P j ), where γ i is an ALP-daughter photon, γ P j is a photon in the ROE, and γ i and γ P j are chosen so that m(γ i γ P j ) is closest to the nominal mass of each of P = π 0 , η, η ; and, multiplicity of neutral candidates in the event.
The BDT score distributions for data, signal MC, and background MC are provided in Ref. [44].For our final signal region selection, we apply the criteria on the two BDT scores shown in Ref. [44], allowing multiple candidates per event.The BDT selection criteria are independent of the ALP mass hypothesis.The signal efficiency estimated from MC varies between 2% for m a = 4.78 GeV to 33% for m a = 0.3 GeV.The resulting m γγ distribution is shown in Fig. 1.
The background is dominated by continuum events and by peaking contributions from B ± → K ± h 0 and B ± → π ± h 0 decays where h 0 = π 0 , η, η .Non-resonant B ± → K ± γγ decays and B ± → K * γ, K * → Kγ decays are negligible, as they have total branching fractions < ∼ 10 −7 [45,46] and do not give a peak in m γγ .The B ± → K ± η c , η c → γγ decay is not included in our background MC; we observe an excess at the η c mass, with a local significance of 2.6σ as determined by the signal extraction procedure defined below.The measured product branching fraction is consistent with the world average value of B(B ± → K ± η c )B(η c → γγ) [46].Because of the relatively small η c background compared to the π 0 , η, and η , we do not exclude signal mass hypotheses in the vicinity of the η c mass.
We extract the signal yield of promptly decaying ALPs by performing a series of unbinned maximum likelihood fits of a hypothetical signal peak over a smooth background to the data shown in Fig. 1.We perform fits for 461 signal mass hypotheses with a scan step size equal to the signal resolution, σ γγ .The latter is determined by fitting the signal sample at each simulated ALP mass with a double-sided Crystal Ball function [47] and interpolating the results to the remaining mass hypotheses.The resolution ranges from 8 MeV near m a = 0.175 GeV to 14 MeV near m a = 2 GeV, and decreasing back to 2 MeV near m a = 4.78 GeV as a result of the constraint imposed on the mass of the B ± meson candidate in the kinematic fit.The MC predictions are validated using a sample of B ± → K ± π 0 and B ± → K ± η decays.The simulated π 0 and η mass resolutions agree with the data to within 3%.
Each unbinned likelihood fit is performed over an m γγ interval with a width in the range (24-60)σ γγ .The mass-dependent interval width is chosen to be sufficiently broad as to fix the continuum background shape.We have verified that our results are independent of minor variations of the fit interval widths.The probability density function (pdf) includes contributions from signal, continuum background components, and, where needed, peaking components describing the π 0 , η, η , and η c .
The signal pdf is described by a non-parametric kernel density function modeled from signal MC and extrapolated between simulated mass points [48].The continuum background is modeled for m a < 4 GeV by the sum of a template derived from background MC and a first-order polynomial, with the normalization determined from the fit.At higher masses, only the first-order polynomial is needed to model the background.The data-to-MC ratio is approximately constant over each fit interval, and the residual differences are accommodated by the linear poly- nomial.The shapes of the π 0 , η, and η resonances are also modeled from background MC, while the η c is modeled using the signal MC mass distribution with a width broadened to match the η c natural linewidth.For the π 0 , η, and η background components, the normalization is determined from the fit to data, while the normalization of the η c component is fixed to the product of the worldaverage value of B(B ± → K ± η c )B(η c → γγ) and the signal efficiency evaluated at this mass.This allows us to measure an ALP signal rate for m a ≈ m ηc while simultaneously accounting for events from B ± → K ± η c , η c → γγ decays.We have verified that our signal extraction procedure is robust against changes in the background model by varying the order of the polynomial component of the continuum background.
To assess systematic uncertainties in the MC-derived continuum and peaking background components, we fit the relative normalizations of different background components (continuum q q, B + B − , B 0 B0 ) to data rather than fixing each component's normalization to match the luminosity of the total data set, and we repeat our signal extraction procedure with the re-weighted MC-derived templates.We also propagate the uncertainties in the resolution of the peaking components and in the uncertainties in the world-average value of the η c linewidth.For the η c model, we assess a systematic uncertainty originating from uncertainties in B(B ± → K ± η c )B(η c → γγ) by varying the η c normalization within the uncertainties in the world-average value.The systematic uncertainty in the signal yield resulting from variations in the continuum (peaking) background shape due to re-fitting the component normalizations is estimated to be, on average, 1% (2%) of the corresponding statistical uncertainty.
We further assess systematic uncertainties associated with our signal model.The systematic uncertainty in the signal yield derived from our extrapolation procedure between simulated mass points is estimated to be, on average, 4% of the corresponding statistical uncertainty.We assess a signal resolution systematic uncertainty by repeating our fits with a signal shape whose width is varied within the mass resolution uncertainty, leading to a signal resolution systematic uncertainty that is, on average, 3% of the statistical uncertainty.A 6% relative systematic uncertainty in the signal efficiency is derived from the data-to-MC ratio for events in the vicinity of the η resonance.The fitted signal yields and statistical significances are shown in Fig. 2. The largest local significance of 3.3σ is observed near m a = 3.53 GeV with a global significance of 1.1σ after including trial factors [49], consistent with the null hypothesis.Background-only fits to the m γγ spectrum are shown over the whole mass range in Ref. [44].
To further validate the signal extraction procedure, we measure the B ± → K ± h 0 , h 0 → γγ (h 0 = π 0 , η, η , η c ) product branching fractions by treating the peaks as signal, extracting the number of events in the peak using the fitting procedure described above, and subtracting non-peaking background whose magnitude is determined from MC.The results are found to be compatible with the current world averages [46] within uncertainties.
In the absence of significant signal, Bayesian upper limits at 90% confidence level (CL) on B(B ± → K ± a)B(a → γγ) are derived with a uniform positive prior in the product branching fraction.We have verified that the limits are robust with respect to the choice of prior.The systematic uncertainty is included in the limit calculation by convolving the likelihood function with a Gaussian having a width equal to the systematic uncertainty.Uncertainties in the luminosity (0.6%) [34] and from the limited statistical precision of simulated samples (1%) are included as well.The resulting limits on the branching fraction product assuming promptly decaying ALPs are displayed in Fig. 3.
Our search targets promptly decaying ALPs.However, ALPs can be long lived at small masses and coupling, and we assess how our signal efficiency and resolution are affected for ALP proper decay lengths of cτ a = 1, 10, and 100 mm.These decay lengths range from nearly prompt decays for which the efficiency and resolution are comparable to the zero-lifetime signal, through to the longest values to which our analysis is sensitive.We measure the B ± → K ± a branching fraction for each decay length.We restrict this study to the mass range for which we obtain sensitivity to couplings that give rise to long-lived ALPs, namely m a < 2.5 GeV.Long-lived ALPs induce a non-negligible bias in the measurement of m γγ , and the resolution is significantly impacted, ranging from 15 MeV near m a = 0.175 GeV to 28 MeV near 2 GeV for cτ a = 100 mm.For cτ a = 100 mm, we only consider mass hypotheses m a ≥ 0.2 GeV, because there is a significant overlap between the signal mass distribution and the π 0 background for lower ALP masses.
The signal is extracted in the same manner as for the promptly decaying ALP, and the fitted signal yields and local statistical significances are shown in Ref. [44].The largest local significance is found to be at m a = 1.10 GeV and cτ a = 10 mm, with a global significance of less than one standard deviation.Systematic uncertainties are assessed in the same manner as for the prompt analysis.The systematic uncertainty in the signal yield resulting from variations in the continuum (peaking) background shape due to re-fitting the component normalizations is larger for long-lived ALPs because of the long tail induced by the bias in the measurement of the signal m γγ The 90% CL upper limits on the coupling gaW as a function of the ALP mass (red), together with existing constraints [29] (blue, green, brown, and grey).
distribution, and is estimated to be, on average, 16% (24%) of the corresponding statistical uncertainty for cτ a = 100 mm.The other systematic uncertainties are comparable in magnitude to the values for prompt ALPs, and the total systematic uncertainty is subdominant to the statistical uncertainty for all signal mass hypotheses.
The 90% CL upper limits on B(B ± → K ± a)B(a → γγ) are plotted in Fig. 4. The limits degrade at cτ a = 100 mm because of the broadening of the signal shape and lower efficiency.The cτ a dependence of the limit is less pronounced at higher masses because the ALP is less boosted, leading to a shorter decay length in the detector.We use an interpolating function to obtain product branching fraction limits for intermediate lifetimes between those shown in Fig. 4.
The 90% CL limits on the ALP coupling g aW are presented in Fig. 5.For each ALP mass hypothesis, we determine the value of g aW such that the calculated branching fraction is equal to the 90%-CL-excluded branching fraction for the lifetime predicted using the same value of g aW .This is the excluded value of g aW shown in Fig. 5.The 90% CL bounds on g aW extend below 10 −5 GeV −1 for many ALP masses, improving current constraints by more than two orders of magnitude.The strongest limit on the coupling at m a = 0.2 GeV corresponds to a lifetime of cτ a = 100 mm, decreasing to cτ a = 1 mm at m a = 2.5 GeV.Along with our limit, we show in Fig. 5 existing constraints derived in Ref. [29] from LEP, beam dump, and K → πγγ searches.We have also reinterpreted a search for K ± → π ± X with invisible X [50], which applies to our model if the ALP is sufficiently long lived that it decays outside of the detector.
In summary, we report the first search for axion-like particles in the process B ± → K ± a, a → γγ.The results strongly constrain ALP couplings to electroweak gauge bosons, improving upon current bounds by several or-ders of magnitude, except in the vicinity of the π 0 , η, and η resonances.Our results demonstrate the sensitivity of flavor-changing neutral current probes of ALP production, which complement existing searches for the ALP coupling to photons below the B meson mass.
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible.The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR.The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them.This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), the Commissariat à l'Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the Research Council of Norway, the Ministry of Education and Science of the Russian Federation, Ministerio de Economía y Competitividad (Spain), the Science and Technology Facilities Council (United Kingdom), and the Binational Science Foundation (U.S.-Israel).Individuals have received support from the Marie-Curie IEF program (European Union) and the A. P. Sloan Foundation (USA).

EPAPS Material
The following includes supplementary material for the Electronic Physics Auxiliary Publication Service.

FIG. 1 :
FIG.1:The diphoton mass distribution of ALP candidates, together with Monte Carlo predictions of leading background processes normalized to the data luminosity.

5 FIG. 2 :
FIG. 2:The distribution of signal events (Ns) and local signal significance (Ss) from fits as a function of ma for prompt ALP decays.The vertical gray bands indicate the regions excluded from the search in the vicinity of the π 0 , η, and η masses.

FIG. 3 :
FIG.3: 90% CL upper limits on the B ± → K ± a branching fraction assuming promptly decaying ALPs.The vertical gray bands indicate the regions excluded from the search in the vicinity of the π 0 , η, and η masses.

FIG. 4 :
FIG.4:The 90% CL upper limits on the B ± → K ± a branching fraction for ma < 2.5 GeV and cτa between 0 and 100 mm.The vertical gray bands indicate the regions excluded from the search in the vicinity of the π 0 , η, and η masses.
FIG.5:The 90% CL upper limits on the coupling gaW as a function of the ALP mass (red), together with existing constraints[29] (blue, green, brown, and grey).

25 •
Adopted same BDT cuts for all signal masses• Cut of > 0.13 on continuum-trained BDT output, > 0.15 on -trained BDT output e + e !B + B < l a t e x i t s h a 1 _ b a s e 6 4 = " FIG. 6: Boosted decision tree classifier scores for ALP candidates in data, signal MC with ma = 1 GeV normalized to a branching fraction of 10 −4 , and background MC normalized to the data luminosity.The classifier scores are for BDTs trained on (left) continuum MC; (right) B + B − MC.The selections used in the analysis are indicated for each distribution by the arrow.

4 FIG. 7 :
FIG. 7: Fits of our background model to data for a sample of fit intervals used in the ALP search.The green line indicates the best-fit background model, and the points represent data.The fit interval width is determined by the signal resolution as outlined in the main-body text.

5 FIG. 8 :
FIG. 8:The distribution of signal events (Ns) and local signal significance (Ss) from fits as a function of ma for ALP proper decay lengths of 1 mm, 10 mm, and 100 mm.The vertical gray bands indicate the regions excluded from the search in the vicinity of the π 0 , η, and η masses.
GeV, where √ s denotes the center-of-mass (CM) energy, p B and E B are the B ± momentum and energy in the lab frame, E *