Search for a dark photon and an invisible dark Higgs boson in $\mu^+\mu^-$ and missing energy final states with the Belle II experiment

The dark photon $A^\prime$ and the dark Higgs boson $h^\prime$ are hypothetical particles predicted in many dark sector models. We search for the simultaneous production of $A^\prime$ and $h^\prime$ in the dark Higgsstrahlung process $e^{+}e^{-}\rightarrow A^\prime \, h^\prime$ with $A^\prime \rightarrow \mu^+\mu^-$ and $ h^\prime$ invisible in electron-positron collisions at a center-of-mass energy of 10.58 GeV collected by the Belle II experiment in 2019. With an integrated luminosity of 8.34 fb$^{-1}$, we observe no evidence for signal. We obtain exclusion limits at 90% Bayesian credibility in the range of 1.7--5.0 fb on the cross section and in the range of $1.7 \times10^{-8}$--$200 \times10^{-8}$ on the effective coupling $\varepsilon^2 \times \alpha_D$ for the $A^\prime$ mass in the range of 4.0 GeV/$c^2$ $<M_{A^\prime}<9.7$ GeV/$c^2$ and for the $h^\prime$ mass $M_{h^\prime}<M_{A^\prime}$, where $\varepsilon$ is the mixing strength between the standard model and the dark photon and $\alpha_D$ is the coupling of the dark photon to the dark Higgs boson. Our limits are the first in this mass range.

The dark photon A and the dark Higgs boson h are hypothetical particles predicted in many dark sector models. We search for the simultaneous production of A and h in the dark Higgsstrahlung process e + e − → A h with A → µ + µ − and h invisible in electron-positron collisions at a center-ofmass energy of 10.58 GeV in data collected by the Belle II experiment in 2019. With an integrated luminosity of 8.34 fb −1 , we observe no evidence for signal. We obtain exclusion limits at 90% Bayesian credibility in the range of 1.7-5.0 fb on the cross section and in the range of 1.7 × 10 −8 -200 × 10 −8 on the effective coupling ε 2 × αD for the A mass in the range of 4.0 GeV/c 2 < M A < 9.7 GeV/c 2 and for the h mass M h < M A , where ε is the mixing strength between the standard model and the dark photon and αD is the coupling of the dark photon to the dark Higgs boson. Our limits are the first in this mass range. This paper is dedicated to the memory of Cate MacQueen.
Several astrophysical observations in the last decades suggest the existence of a large quantity of dark matter in the universe coupled with ordinary matter, at least through gravitational interactions. In recent years, theoretical models (commonly called dark sector models) with light particles mediating new interactions have gained considerable attention as solutions to the long-standing problem of reproducing the observed relic density of dark matter. A well-motivated model predicts the existence of an additional massive vector gauge boson, a dark photon A , coupled to the standard model (SM) only through its kinetic mixing ε with the hypercharge field [1][2][3][4][5][6][7][8][9][10]. In this model, dark matter particles can annihilate to SM particles, and vice versa, through the exchange of dark photons: this process contributes to the relic density of dark matter so as to match the observed value. The mass of the A can arise from spontaneous symmetry breaking, which introduces a new scalar particle: a dark Higgs boson h [11]. Several searches for the A have set upper limits on ε [12][13][14][15][16][17][18][19][20][21][22][23][24]; they depend on the mass of the dark photon and are typically of the order ε 2 ≤ 5 × 10 −7 for masses below 10.6 GeV/c 2 .
In this Letter, we search for the so-called dark Higgsstrahlung process e + e − → A h using data collected by the Belle II experiment at the SuperKEKB collider. We consider the minimal secluded model of Ref. [11], in which either the dark Higgs boson does not mix with the SM Higgs boson or its mixing can be neglected and any additional particles (in particular dark matter candidates) are heavier than both A and h . The cross section for dark Higgsstrahlung is proportional to ε 2 ×α D , where α D is the coupling constant of the secluded model [11]. Two scenarios exist, differentiated by the hierarchy of the masses M A and M h . If M h > M A , then the h decays dominantly to an A A ( * ) pair (where A ( * ) can be virtual), which is a final state searched for by the BaBar [25] and Belle [26] experiments. If M h < M A , then the h is long lived and invisible because it does not interact with the detector material. We focus on the latter scenario, which was previously investigated by KLOE [27] but in a much smaller mass range than is accessible to Belle II. Our search is limited to the decay of the A into a muon pair.
The dark Higgsstrahlung process studied here produces a pair of oppositely charged muons with a mass M µµ distributed around M A and a recoil against them with a mass M recoil that, in the absence of radiated photons, is distributed around M h . In the following, recoil quantities are computed against the dimuon system. The accessible search region is the mass plane M recoil -M µµ : it has a triangular shape, limited by the conditions M recoil < M µµ (corresponding to M h < M A ) and M recoil + M µµ ≤ √ s/c 2 , resulting from energy conservation. Our analysis uses events with exactly two tracks, identified as muons, and negligible extra energy. The backgrounds are SM processes that produce final states with two tracks identified as muons and missing energy. These are dominantly e + e − → µ + µ − (γ) with typically one or more photons undetected due to inefficiency or limited acceptance; e + e − → τ + τ − (γ) with both τ leptons decaying to muons and neutrinos; and e + e − → e + e − µ + µ − with electrons outside acceptance. We search for the signal as a narrow enhancement in the two-dimensional M 2 recoil -M 2 µµ distribution. We scan the search region for local excesses above the expected background with a counting technique that uses a set of twodimensional M 2 recoil -M 2 µµ overlapping windows. Event selection is optimized using simulated events prior to examining data.
The Belle II detector [28] operates at the SuperKEKB electron-positron e + e − collider [29], located at KEK in Tsukuba, Japan. The beam energies are 7 GeV for e − and 4 GeV for e + , resulting in a boost of βγ = 0.28 of the center-of-mass (c.m.) frame relative to the laboratory frame. Data used in this analysis were collected in 2019 at a c.m. energy √ s corresponding to the mass of the Υ(4S) resonance, for an integrated luminosity of 8.34 ± 0.08 fb −1 (see Ref. [30] for the description of the luminosity measurement technique).
The Belle II detector consists of several subdetectors arranged around the beam pipe in a cylindrical structure. Subdetectors relevant for this analysis are briefly described here in the order of innermost out; a description of the full detector is given in Refs. [28,31]. The innermost subdetector is the vertex detector, which consists of two inner layers of silicon pixels and four outer layers of silicon strips. The second pixel layer was only partially installed for the data sample we analyze, covering only one sixth of the azimuthal angle. The main tracking subdetector is a large helium-based small cell drift chamber (CDC). The relative charged-particle transverse momentum resolution is typically 0.1% p T ⊕ 0.3%, with p T expressed in GeV/c. An electromagnetic calorimeter (ECL) consists of a barrel and two end caps made of CsI(Tl) crystals. A superconducting solenoid, situated outside of the calorimeter, provides a 1.5 T magnetic field. A K 0 L and muon subdetector (KLM) is made of iron plates, which serve as a magnetic flux-return yoke, alternated with resistive-plate chambers and plastic scintillators in the barrel as well as with plastic scintillators only in the end caps. The longitudinal direction, the transverse plane and the polar angle θ are defined with respect to the detector's solenoidal axis in the direction of the electron beam. In the following, quantities are defined in the laboratory frame unless specified otherwise.
The identification of muons uses criteria that rely mostly on track penetration in the KLM for momenta larger than 0.7 GeV/c and on information from the CDC and ECL otherwise. Electrons are identified mostly by comparing measured momenta in the CDC with energies of the associated ECL clusters. Photons are identified as ECL clusters with energies greater than 100 MeV that are not associated with tracks. Details of the particle reconstruction and identification algorithms are in Refs. [31][32][33].
The detector geometry and interactions of final-state particles with detector materials are simulated using Geant4 [44] and the Belle II Analysis Software Framework (basf2) [45]. Both real data and simulated events are reconstructed and analyzed using basf2.
The search uses an online event selection (trigger) that requires at least one pair of tracks in a restricted polar angle acceptance, θ ∈ [37, 120] • , with an azimuthal opening angle ∆φ larger than 90 • and rejects events consistent with Bhabha scattering through a dedicated veto based on the pattern of energy depositions in the ECL. The efficiency of this trigger for ∆φ > 90 • is 89%, measured in µ + µ − (γ) events. The trigger requirement on the opening angle leads to large signal inefficiencies for To suppress misreconstructed and beam-induced background tracks, we require that the transverse and longitudinal projections of their distances of closest approach to the interaction point be smaller than 0.5 and 2.0 cm, respectively. We require that events have exactly two oppositely charged particles, identified as muons, with polar angles θ ∈ [37, 120] • and with ∆φ larger than 90 • to match the trigger requirements. In addition, the energy of each ECL cluster associated with a muon track must be below 1.5 GeV to suppress background events from muon pairs, with final-state radiation and unresolved photons close to the muons. We require that the recoil momentum point into the ECL barrel, θ ∈ [32, 125] • , to exclude regions where photons from radiative backgrounds can escape undetected. This selection also increases the signal-to-background ratio by suppressing µ + µ − (γ) and e + e − µ + µ − processes that, unlike the signal, have recoil momenta dominantly in the forward direction. To reduce radiative muon-pair backgrounds, we require that the total energy of all photons be less than 0.4 GeV and no photon be within 15 • of the recoil momentum. To suppress background events from µ + µ − (γ) and e + e − µ + µ − processes, we require that the transverse c.m. frame momentum of the muon pair be greater than 0.1 GeV/c.
We count events in 9003 partially overlapping regions (search windows) of the two-dimensional space of squared dimuon and recoil masses, which span the M recoil -M µµ search plane. We search for signal by comparing the observed yields with expectations from known backgrounds. The squared dimuon and recoil masses are negatively correlated in signal events and, to a lesser extent, in background events, with the correlation varying across the plane. Initial-state radiation partially spoils this correlation because it affects the M recoil distribution only. The search window boundaries are chosen as ellipses in the M 2 recoil -M 2 µµ plane that take the local correlation into account. Each window is centered at one of the values of (M 2 A , M 2 h ) used to produce a simulated dataset. We fit a sum of two two-dimensional Gaussian distributions that share a common mean and correlation to simulated data without ISR. The elliptical search window boundaries correspond to the two-dimensional two-standarddeviation contours resulting from the fit. Search windows partially overlap to maximize signal efficiency, with an overlap in area typically around 75%. The average fraction of signal events retained in a search window is 71%, with variations due to the different effects induced by ISR depending on M A and M h . Accounting for correlation in defining the window increases the signal-tobackground ratio by a factor of three to five.
A final selection is based on the helicity angle η, defined as the angle in the dimuon rest frame between the momentum direction of the c.m. system and the momentum direction of the µ − . The η distribution for the signal is that of a massive vector particle decaying into two fermions. For an unpolarized A , the distribution of C η = | cos η | is uniform. In background events, C η peaks at one because the muons come either from independent decays, as in τ + τ − (γ), or from physics processes of a different nature, as in µ + µ − (γ) and e + e − µ + µ − . We exclude high values of C η , typically larger than 0.9, by maximizing the figure of merit of Ref. [46] in each search window.
The resulting signal efficiency depends on the masses of the A and h and ranges between 10 and 25% for M A > 4 GeV/c 2 . The efficiency drops considerably for M A < 4 GeV/c 2 and becomes too low to allow analysis for M A < 1.65 GeV/c 2 . We restrict our analysis to the range M A > 1.65 GeV/c 2 . Backgrounds are typically reduced by factors of 10-1000. The surviving background events come 78.5% from µ + µ − (γ), 18.5% from τ + τ − (γ), and 3% from e + e − µ + µ − : they dominantly populate the regions of the search plane close to the kinematical limit M recoil + M µµ ≈ √ s/c 2 , where the C η selection is less effective and µ + µ − (γ) contributes more, leaving most of the mass plane almost background free.
The distribution of the observed event counts in the search windows is shown in Fig. S6. Because of overlap, events can be counted in multiple search windows: this correlation creates groups of isolated, sparsely populated windows. A total of 28 985 events pass all the selection criteria, which is in agreement with the expectation from simulation, 28 486 ± 331. The sum of the event counts inside all mass windows is 78 740. On average, each event is contained in 2.7 search windows.
We search for an excess of signal above the expected background independently in the 9003 search windows. Event counts N in a search window are interpreted ac-cording to the relation N = sig × L × σ + B, where σ is the cross section for the dark Higgsstrahlung process e + e − → A h with A → µ + µ − and h invisible, L = L dt is the integrated luminosity, and sig and B are the signal efficiency and the expected background inside the window. Both sig and B are determined from simulation and are subject to systematic uncertainties. Several sources of systematic uncertainties affecting the signal efficiency and the background estimate are taken into account. They are studied by comparing data and simulation on two control samples that emulate the two main background processes µ + µ − (γ) and τ + τ − (γ) and on the search µ + µ − sample. The µµγ control sample contains events that pass all selection criteria except for the veto on the presence of photons, which is replaced by a requirement that a photon with energy greater than 1 GeV be reconstructed in the barrel of the ECL. This sample is dominantly composed of µ + µ − (γ) events. The eµ control sample contains events passing all the selection criteria but with an identified electron replacing an identified muon. This sample is almost entirely composed of τ + τ − (γ) events. We split the mass plane into six nonoverlapping macroregions, with each dominantly populated by a single background source. Events are counted in the data and simulation for each macroregion and their discrepancies used to evaluate systematic uncertainties. When using the search µ + µ − sample, a region ten times larger than the search window under study is excluded from the counting to avoid that the presence of a signal may bias the result.
Uncertainties affecting the background due to the trigger, luminosity, tracking efficiency, muon identification, cross sections, and the selection criteria are collectively evaluated through the macroregion studies before applying the C η selection. Over most of the mass plane, discrepancies between the data and simulation, for both the control and the search µ + µ − samples, are of the order of 2%, which are treated as relative systematic uncertainties. In a small region, where M A > 9 GeV/c 2 , the search µ + µ − sample data yields are 9.1% lower than the simulation yields. In this region, there are also discrepancies between the data and simulation in the shapes of mass distributions that lead to an additional relative systematic uncertainty of 9.3%: we sum it quadratically with the 9.1% normalization uncertainty. We assume the 2% and 9.1% uncertainties also hold for the signal efficiency below and above M A = 9 GeV/c 2 , respectively.
Uncertainties affecting the background due to the C η selection are evaluated by comparing the data and simulation. The numbers of observed and expected events for both the control and the search µ + µ − samples agree within 1%, which we use as a relative systematic uncertainty due to this source. We evaluate the contribution due to this effect on the signal efficiency to be negligible.
We also include systematic uncertainties due to discrepancies in the dimuon and recoil mass resolutions in data and simulation. A modified µµγ control sample is used to check the mass distributions in the region of the J/ψ resonance: the requirement on the opening angle ∆φ is released and an ECL-only trigger is used, which is allowed by the presence of the photon. This trigger requires that the total energy deposition in the barrel and in the forward end cap exceed 1 GeV. The search µ + µ − sample is also used to check the mass distributions in the region where the dimuon mass is close to √ s/c 2 . We find differences between the data and simulation in the dimuon and recoil mass resolutions of no more than 10%. Their effects on the relative signal efficiency range between 1 and 5%, depending on the masses, with an average of 2.4%.
We evaluate a systematic uncertainty on the signal efficiency due to M A and M h not coinciding with a window center by recalculating the efficiency with the two mass values randomly varied to points near the window center. The signal efficiency varies 2% on average, and no more than 5%, which we assign as a relative systematic uncertainty.
Finally, a relative systematic uncertainty of 4% on the theoretical prediction of the A decay branching fraction to muons is used when interpreting results in terms of the coupling product ε 2 × α D . This uncertainty comes dominantly from uncertainties on the measured ratio of cross sections for the production of hadrons or muons in e + e − collisions, which enter in the A width theoretical calculation [11].
The average total relative systematic uncertainties are 2.2 and 5.4% on the background and signal efficiencies, respectively. They rise up to 12.7 and 11.3% in the region M A > 9 GeV/c 2 .
We search for excesses in data in each window separately with both a Bayesian technique based on Bayes factors [47] and a frequentist technique based on signif-icances from one-sided Gaussian integral transformation of p values. Background expectations and signal efficiencies are assumed from the simulation. Systematic uncertainties are taken into account as correlated Gaussian smearings of background expectations and signal efficiencies, with widths equal to the estimated uncertainties. We choose thresholds of 80 for the Bayes factor and of 3.5σ for the significance before inspecting the data: they are larger than normally used, because we expect a relevant look-elsewhere effect [48,49] due to the high number of search windows. We find only one case of a local significance above the threshold, 3.7σ, which also corresponds to the highest Bayes factor of 45.6. It is in the search window centered at M µµ = 5.44 GeV/c 2 and M recoil = 3.18 GeV/c 2 . Taking into account the look-elsewhere effect, this excess has a global significance below 1σ, showing no evidence for signal.
We compute Upper Limits (UL) at a 90% Bayesian credibility level (CL) on the cross section for the dark Higgsstrahlung process e + e − → A h with A → µ + µ − and h invisible as a function of M A and M h using the Bayesian Analysis Toolkit software package [50]. We assume uniform priors for all positive values of the cross section, Poissonian likelihoods for the number of observed and simulated events, and Gaussian smearing to model systematic uncertainties, accounting for their correlations. The result is shown in Fig. 2.  Our results are dominated by their statistical uncertainties. In most of the search plane, systematic uncertainties degrade the upper limits by less than 1%. Only in the small region where M A > 9 GeV/c 2 are system-atic uncertainties significant, worsening the upper limits by 25%. We test for prior dependence of the results by using logarithmic priors for the cross section and find differences smaller than 3%. Additional plots and detailed numerical results are provided in the Supplemental Material [51]. In summary, we search for the dark Higgsstrahlung process e + e − → A h with A → µ + µ − and h invisible in a data sample of electron-positron collisions at 10.58 GeV collected by Belle II at SuperKEKB in 2019, corresponding to an integrated luminosity of 8.34 fb −1 . We find no significant excess above the expected background and set upper limits on the cross section and coupling ε 2 × α D for M A between 1.65 and 10.51 GeV/c 2 and M h < M A . Our limits are the first in this mass range. The excluded region is much larger than that previously covered by other experiments [27]. Our 90% CL upper limits range between 1.7 and 5.0 fb for the cross section and between 1.7×10 −8 and 200×10 −8 for the coupling for 4.0 GeV/c 2 < M A < 9.7 GeV/c 2 and M h < M A . For specific values of α D and assuming the existence of a light invisible dark Higgs, our results can be interpreted as upper limits on ε 2 and compared with limits obtained by other experiments. With α D = 1, our constraints would improve on previous searches [22] across almost the full mass range. For α D = 0.1, this conclusion would still hold in a substantial part of the mass range. These results can be interpreted in a wider class of models compared to that of Ref. [11], for example those with a long-lived invisible h that mixes with the SM Higgs boson [52,53]. Figure S1: Example of a two-dimensional squared mass distribution for a given signal hypothesis (M A =6.5 GeV/c 2 , M h =2.7 GeV/c 2 ). The upward tail in the squared recoil mass distribution is due to ISR. Also shown is the elliptical search window contour.