Search for a Dark Leptophilic Scalar in e + e − Collisions

J. P. Lees, V. Poireau, V. Tisserand, E. Grauges, A. Palano, G. Eigen, D. N. Brown, Yu. G. Kolomensky, M. Fritsch, H. Koch, T. Schroeder, R. Cheaib, C. Hearty, T. S. Mattison, J. A. McKenna, R. Y. So, V. E. Blinov, A. R. Buzykaev, V. P. Druzhinin, V. B. Golubev, E. A. Kozyrev, E. A. Kravchenko, A. P. Onuchin, S. I. Serednyakov, Yu. I. Skovpen, E. P. Solodov, K. Yu. Todyshev, A. J. Lankford, B. Dey, J. W. Gary, O. Long, A. M. Eisner, W. S. Lockman, W. Panduro Vazquez, D. S. Chao, C. H. Cheng, B. Echenard , K. T. Flood, D. G. Hitlin, J. Kim, Y. Li, D. X. Lin, T. S. Miyashita, P. Ongmongkolkul, J. Oyang, F. C. Porter, M. Röhrken, Z. Huard, B. T. Meadows, B. G. Pushpawela, M. D. Sokoloff, L. Sun, J. G. Smith, S. R. Wagner, D. Bernard, M. Verderi, D. Bettoni, C. Bozzi, R. Calabrese, G. Cibinetto, E. Fioravanti, I. Garzia, E. Luppi, V. Santoro, A. Calcaterra, R. de Sangro, G. Finocchiaro, S. Martellotti, P. Patteri, I. M. Peruzzi, M. Piccolo, M. Rotondo, A. Zallo, S. Passaggio, C. Patrignani, B. J. Shuve, H. M. Lacker, B. Bhuyan, U. Mallik, C. Chen, J. Cochran, S. Prell, A. V. Gritsan, N. Arnaud, M. Davier, F. Le Diberder, A. M. Lutz, G. Wormser, D. J. Lange, D. M. Wright, J. P. Coleman, E. Gabathuler, D. E. Hutchcroft, D. J. Payne, C. Touramanis, A. J. Bevan, F. Di Lodovico, R. Sacco, G. Cowan, Sw. Banerjee, D. N. Brown, C. L. Davis, A. G. Denig, W. Gradl, K. Griessinger, A. Hafner, K. R. Schubert, R. J. Barlow,32,∥ G. D. Lafferty, R. Cenci, A. Jawahery, D. A. Roberts, R. Cowan, S. H. Robertson, R. M. Seddon, N. Neri, F. Palombo, L. Cremaldi, R. Godang, D. J. Summers, P. Taras, G. De Nardo, C. Sciacca, G. Raven, C. P. Jessop, J. M. LoSecco, K. Honscheid, R. Kass, A. Gaz, M. Margoni, M. Posocco, G. Simi, F. Simonetto, R. Stroili, S. Akar, E. Ben-Haim, M. Bomben, G. R. Bonneaud, G. Calderini, J. Chauveau, G. Marchiori, J. Ocariz, M. Biasini, E. Manoni, A. Rossi, G. Batignani, S. Bettarini, M. Carpinelli, G. Casarosa, M. Chrzaszcz, F. Forti, M. A. Giorgi, A. Lusiani, B. Oberhof, E. Paoloni, M. Rama, G. Rizzo, J. J. Walsh, L. Zani, A. J. S. Smith, F. Anulli, R. Faccini, F. Ferrarotto, F. Ferroni, A. Pilloni, G. Piredda, C. Bünger, S. Dittrich, O. Grünberg, M. Heß, T. Leddig, C. Voß, R. Waldi, T. Adye, F. F. Wilson, S. Emery, G. Vasseur, D. Aston, C. Cartaro, M. R. Convery, J. Dorfan, W. Dunwoodie, M. Ebert, R. C. Field, B. G. Fulsom, M. T. Graham, C. Hast, W. R. Innes, P. Kim, D.W. G. S. Leith, S. Luitz, D. B. MacFarlane, D. R. Muller, H. Neal, B. N. Ratcliff, A. Roodman, M. K. Sullivan, J. Va’vra, W. J. Wisniewski, M. V. Purohit, J. R. Wilson, A. Randle-Conde, S. J. Sekula, H. Ahmed, M. Bellis, P. R. Burchat, E. M. T. Puccio, M. S. Alam, J. A. Ernst, R. Gorodeisky, N. Guttman, D. R. Peimer, A. Soffer, S. M. Spanier, J. L. Ritchie, R. F. Schwitters, J. M. Izen, X. C. Lou, F. Bianchi, F. De Mori, A. Filippi, D. Gamba, L. Lanceri, L. Vitale, F. Martinez-Vidal, A. Oyanguren, J. Albert, A. Beaulieu, F. U. Bernlochner, G. J. King, R. Kowalewski, T. Lueck, I. M. Nugent, J. M. Roney, R. J. Sobie, N. Tasneem, T. J. Gershon, P. F. Harrison, T. E. Latham, R. Prepost, and S. L. Wu

Many scenarios of physics beyond the standard model predict the existence of new gauge singlets, which might be substantially lighter than the weak scale. The experimental constraints on additional scalars with masses in the MeV to GeV range could be significantly weakened if they interact predominantly with leptons rather than quarks. At an e þ e − collider, such a leptophilic scalar (ϕ L ) would be produced predominantly through radiation from a τ lepton. We report herein a search for e þ e − → τ þ τ − ϕ L , ϕ L → l þ l − (l ¼ e, μ) using data collected by the BABAR experiment at SLAC. No significant signal is observed, and we set limits on the ϕ L coupling to leptons in the range 0.04 < m ϕ L < 7.0 GeV. These bounds significantly improve upon the current constraints, excluding almost entirely the parameter space favored by the observed discrepancy in the muon anomalous magnetic moment below 4 GeV at 90% confidence level. DOI: 10.1103/PhysRevLett.125.181801 Many theories beyond the standard model (SM) predict the existence of additional scalars, and discovering or constraining their existence might shed light on the physics of electroweak symmetry breaking and the Higgs sector (e.g., see Ref. [1]). Some of these particles may be substantially lighter than the weak scale, notably in the next-to-minimal supersymmetric standard model [2], but also in more generic singlet-extended sectors [3,4]. In the MeV-GeV range, new scalars could mediate interactions between the SM and dark matter, as well as account for the discrepancy in the observed value of the muon anomalous magnetic dipole moment [5][6][7].
The possible coupling of a new scalar ϕ L to SM particles is constrained by SM gauge invariance. In the simplest case, the mixing between the scalar and the SM Higgs boson gives rise to couplings proportional to SM fermion masses. Because the new scalar couples predominantly to heavy-flavor quarks, this minimal scenario is strongly constrained by searches for rare flavor-changing neutral current decays of mesons, such as B → Kϕ and K → πϕ [8]. However, these bounds are evaded if the coupling of the scalar to quarks is suppressed and the scalar interacts preferentially with heavy-flavor leptons [3,4,9,10]. We refer to such a particle as a leptophilic scalar, ϕ L . Its interaction Lagrangian with leptons can be described by where ξ denotes the flavor-independent coupling strength to leptons and v ¼ 246 GeV is the SM Higgs vacuum expectation value [11]. Model independent constraints relying exclusively on the coupling to leptons are derived from a BABAR search for a muonic dark force [12] and beam dump experiments [13,14]. A large fraction of the parameter space, including the region favored by the measurement of the muon anomalous magnetic moment, is still unexplored [3,9,10]. Examples of model dependent bounds from B and h decays for a specific UV completion of the theory can be found in Ref. [3]. The large sample of τ þ τ − pairs collected by BABAR offers a clean environment to study model independent ϕ L production via final-state radiation in e þ e − → τ þ τ − ϕ L . The mass proportionality of the coupling, in particular the feeble interaction with electrons, dictates the experimental signature. For 2m e < m ϕ L < 2m μ , the scalar decays predominantly into electrons, leading to displaced vertices for sufficiently small values of the coupling. Prompt decays into a pair of muons (taus) dominate when 2m μ ≤ m ϕ L < 2m τ (2m τ < m ϕ L ).
We report herein the first search for a leptophilic scalar in the reaction e þ e − → τ þ τ − ϕ L , ϕ L → l þ l − (l ¼ e, μ) for 0.04 < m ϕ L < 7.0 GeV. The cross section for m ϕ L < 2m μ is measured separately for ϕ L lifetimes corresponding to cτ ϕL values of 0, 1, 10, and 100 mm. Above the dimuon threshold, we determine the cross section for prompt ϕ L → μ þ μ − decays. In all cases, the ϕ L width is much smaller than the detector resolution, and the signal can be identified as a narrow peak in the dilepton invariant mass. The search is based on 514 fb −1 of data collected at the ϒð2SÞ, ϒð3SÞ, ϒð4SÞ resonances and their vicinities [15] by the BABAR experiment at the SLAC PEP-II e þ e − collider. The BABAR detector is described in detail elsewhere [16,17]. A sample corresponding to about 5% of the data, called the optimization sample, is used to optimize the search strategy and is subsequently discarded. The remaining data are examined only once the analysis procedure has been finalized.
We select events containing exactly four charged tracks with zero net charge, focusing on τ lepton decays to single tracks and any number of neutral particles. The ϕ L → l þ l − candidates are formed by combining two oppositesign tracks identified as an electron or muon pair by particle identification (PID) algorithms [12,16]. We do not attempt to select a single ϕ L candidate per event, but simply consider all possible combinations. Radiative Bhabha and dimuon events in which the photon converts to an e þ e − pair are suppressed by rejecting events with a total visible mass greater than 9 GeV. We further veto e þ e − → e þ e − e þ e − events by requiring the cosine of the angle between the momentum of the ϕ L candidate and that of the nearest track to be less than 0.98, the missing momentum against all tracks and neutral particles to be greater than 300 MeV, and that there be three or less tracks identified as electrons. We perform a kinematic fit to the selected ϕ L candidates, constraining the two tracks to originate from the same point in space. The dimuon production vertex is required to be compatible with the beam interaction region, while we only constrain the momentum vector of the e þ e − pair to point back to the beam interaction region since the dielectron vertex can be substantially displaced. We select dielectron (dimuon) combinations with a value of the χ 2 per degree of freedom of the fit, χ 2 =n:d:f., less than 3 (12).
A multivariate selection based on boosted decision trees (BDT) further improves the signal purity [28]. The BDTs include variables capturing the typical τ and ϕ L decay characteristics: a well-reconstructed l þ l − vertex, either prompt or displaced; missing energy and momentum due to neutrino emission; relatively large track momenta; low neutral particle multiplicity; and two or more tracks identified as electrons or muons. A few variables are also targeted at specific backgrounds, such as ψð2SÞ → π þ π − J=ψ; J=ψ → μ þ μ − production in initial-state radiation (ISR) events. The ϕ L mass is specifically excluded to limit potential bias in the classifier. A full description of these variables can be found in the Supplemental Material [29]. We train a separate BDT for each of the different final states and cτ ϕL values with signal events modeled using a flat m ϕ L distribution and background events modeled using the optimization sample data.
The final selection of ϕ L candidates for each lifetime selection and decay channel is made by applying a massdependent criterion on the corresponding BDT score that maximizes signal sensitivity. The distributions of the resulting dielectron and dimuon masses for prompt decays are shown in Fig. 1, and spectra for other lifetimes for ϕ L → e þ e − decays are shown in Fig. 2, together with the dominant background components among the set of simulated MC samples. The differences between the data and summed-MC distributions are mainly due to processes that  02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22  are not simulated, dominated by ISR production of highmultiplicity QED and hadronic events as well as twophoton processes. Peaking contributions from J=ψ → μ þ μ − and ψð2SÞ → μ þ μ − decays are also seen, and the corresponding regions are excluded from the signal search. In addition, the dielectron spectrum for cτ ϕL ¼ 1 mm features a broad enhancement from π 0 → γγ decays in which one or both photons convert to e þ e − pairs. Since this feature is much broader than the signal, we do not exclude this mass region but instead treat it as an additional background component. No statistically significant π 0 component is observed for other values of cτ ϕL .
We extract the signal yield for the different lifetimes and final states separately by scanning the corresponding mass spectrum in steps of the signal mass resolution, σ. The latter is estimated by performing fits of a double-sided Crystal Ball function [30] to each signal MC sample and interpolating the results to the full mass range. The resolution ranges from 1 MeV near m ϕ L ¼ 40 MeV for cτ ϕL ¼ 100 mm to 50 MeV near m ϕ L ¼ 7.0 GeV for prompt decays. The signal MC predictions are validated with samples of K 0 S → π þ π − and ψð2SÞ → π þ π − J=ψ; J=ψ → μ þ μ − decays; agreement with the data is observed. For each mass hypothesis, we perform an unbinned likelihood fit over an interval varying between 20 − 50σ (fixed to 60σ) for the dielectron (dimuon) final state. To facilitate the background description, the reduced dimuon mass, m R ¼ ðm 2 μμ − 4m 2 μ Þ 1=2 , is used for 2m μ <m ϕ L < 260 MeV. In that region, fits are performed over a fixed interval m R < 0.2 GeV.
The likelihood function includes contributions from signal, continuum background, and, where needed, peaking components describing the π 0 , J=ψ → μ þ μ − , and ψð2SÞ → μ þ μ − resonances. The signal probability density function (pdf) is described by a nonparametric kernel density function modeled directly from the signal MC mass distribution. An algorithm based on the cumulative density function [31] is used to interpolate the pdf between simulated mass points. The uncertainty associated with this procedure is on average 4% (3%) of the corresponding statistical uncertainty for the dielectron (dimuon) analysis.
The dielectron continuum background is modeled by a second-order polynomial for the cτ ϕL ¼ 100 mm sample and by a second-order polynomial plus an exponential function for the other lifetimes. The peaking π 0 shape for the cτ ϕL ¼ 1 mm sample is determined from sideband data obtained by applying all selection criteria, but requiring the χ 2 =n:d:f: of the kinematic fit to be greater than 3. The peaking π 0 yield and all the continuum background parameters are determined in the fit. To assess systematic uncertainties, we repeat the fits with a third-order polynomial for the continuum background, vary the width of the π 0 shape within its uncertainty, or include a π 0 component for all lifetime samples. The resulting systematic uncertainties are typically at the level of the statistical uncertainty, but dominate the total uncertainty for several mass hypotheses in the vicinity of the π 0 peak.
The reduced mass distribution of the dimuon continuum background is modeled by a third-order polynomial constrained to intersect the origin, and the dimuon continuum is described by a second-order polynomial at higher masses. The shape of the J=ψ and ψð2SÞ resonances are fixed to the predictions of the corresponding MC samples, but their yields cannot be accurately estimated from MC  02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22  simulations and are therefore left to float freely in the fit. A range of AE50 MeV around the nominal J=ψ and ψð2SÞ masses is therefore excluded from the search. The systematic uncertainty associated with the choice of the background model, assessed by repeating the fits with alternative descriptions, is typically at the level of a few events, but can be as large as half the statistical uncertainty for a few points in the high mass region, where statistical precision is limited.
The fitted signal yields and statistical significances are presented in the Supplemental Material [29], together with a few examples of fits. The bias in the fitted values is determined from pseudoexperiments to be negligible compared to the statistical uncertainties. Since the systematic uncertainty associated with the choice of background model can be large in the dielectron channel, we define the signal significance as the smallest of the significance values determined from each background model. Including trial factors, the largest significance is 1.4σ observed near m ϕ L ¼ 2.14 GeV, consistent with the null hypothesis.
The signal efficiency varies between 0.2% for m ϕ L ¼ 40 MeV and cτ ϕL ¼ 100 mm, to 26% around m ϕ L ¼ 5 GeV for prompt decays. The effect of ISR, not included in the samples generated by MADGRAPH, is assessed by simulating events with PYTHIA8 using the matrix elements calculated by MADGRAPH, and reweighting this sample to match the p T distribution of the ϕ L predicted by MADGRAPH. The resulting change in efficiency is found to be about 4% over the full mass range covered by the dielectron channel, and varies from 7% near the dimuon threshold to less than 1% at m ϕ L ∼ 7 GeV. Half the value of each of these differences is propagated as a systematic uncertainty in the signal yield. A correction factor of 0.98 (0.93) on the signal efficiency is included for the dielectron (dimuon) final state to account for differences between data and simulation in track and neutral reconstruction efficiencies, charged particle identification, and trigger efficiencies. The correction for the dielectron channel is derived from a sample of K 0 S → π þ π − produced in τ decays, while that for the dimuon channel is assessed from the BDT score distribution for events in which the missing transverse momentum is greater than 2 GeV, a region where the contribution of unsimulated background components can be neglected. An uncertainty of 3.8% (4.0%) in the dielectron (dimuon) efficiency correction is propagated as a systematic uncertainty.
The e þ e − → τ þ τ − ϕ L cross section at the ϒð4SÞ energy is derived for each lifetime and final state by taking into account the variation of the cross section and signal efficiencies with the beam energy and the ϕ L → l þ l − branching fraction: where N sig denotes the number of signal events, and σ th;nS , ϵ nS and L nS (n ¼ 2, 3, 4) are the theoretical e þ e − → τ þ τ − ϕ L cross section, signal efficiency, and data luminosity at the ϒðnSÞ center-of-mass energy, respectively. In the absence of a significant signal, Bayesian upper limits at 90% confidence level (CL) on the cross sections are derived by assuming a uniform prior in the cross section. Systematic effects are taken into account by convolving the likelihood with a Gaussian having a width equal to the systematic uncertainty. The uncertainties due to the luminosity (0.6%) [15] and the limited statistical precision of the signal MC sample (1%-4%) are incorporated. The resulting limits are shown in Fig. 3. The sharp increase just above the ditau threshold is a reflection of the ϕ L → μ þ μ − branching fraction decreasing quickly in favor of the τ þ τ − final state. The limit on the production cross section of a scalar S without any assumptions on other decay modes is presented in the Supplemental Material [29]. The limits on the scalar coupling ξ, presented in Fig. 4, are derived with an iterative procedure that accounts for a potentially long ϕ L lifetime. An estimate of ξ is first chosen, and the corresponding lifetime and cross section are calculated [3]. These values are compared to the cross section limit interpolated at that lifetime, and the estimate of the coupling is updated. The procedure is iterated until convergence is obtained. Bounds at the level of 0.5-1 are set on the dielectron final state, corresponding to cτ ϕL values of the order of 1 cm, and limits down to 0.2 are derived for dimuon decays. These results are approximately an order of magnitude smaller than the couplings favored by the muon anomalous magnetic moment below the ditau threshold [3] and rule out a substantial fraction of previously unexplored parameter space at 90% C.L.
In summary, we report the first model-independent search for the direct production of a new dark leptophilic scalar. The limits significantly improve upon the previous constraints over a large range of masses, almost entirely ruling out the remaining region of parameter space below the dimuon threshold. More significantly, this search excludes the possibility of the dark leptophilic scalar accounting for the observed discrepancy in the muon magnetic moment for almost all ϕ L masses below 4 GeV. Since these results rely only on ϕ L production in association with tau leptons and its subsequent leptonic decay, they can also be reinterpreted to provide powerful constraints on other leptonically decaying new bosons interacting with tau leptons. The Belle II experiment should be able to further probe these possibilities, and cover the remaining parameter space above the beam dump constraints.