# Decoding the $B\to K\nu\nu$ excess at Belle II: kinematics, operators, and
masses

Kåre Fridell<EMAIL_ADDRESS>KEK Theory Center, Tsukuba, Ibaraki
305–0801, Japan  Department of Physics, Florida State University, Tallahassee,
FL 32306-4350, USA  Mitrajyoti Ghosh<EMAIL_ADDRESS>Department of Physics,
Florida State University, Tallahassee, FL 32306-4350, USA  Takemichi Okui
<EMAIL_ADDRESS>Department of Physics, Florida State University, Tallahassee, FL
32306-4350, USA  KEK Theory Center, Tsukuba, Ibaraki 305–0801, Japan  Kohsaku
Tobioka<EMAIL_ADDRESS>Department of Physics, Florida State University,
Tallahassee, FL 32306-4350, USA  KEK Theory Center, Tsukuba, Ibaraki 305–0801,
Japan

###### Abstract

An excess in the branching fraction for $B^{+}\to K^{+}\nu\nu$ recently
measured at Belle II may be a hint of new physics. We perform thorough
likelihood analyses for different new physics scenarios such as $B\to KX$ with
a new invisible particle $X$, or $B\to K\chi\chi$ through a scalar, vector, or
tensor current with $\chi$ being a new invisible particle or a neutrino. We
find that vector-current 3-body decay with $m_{X}\simeq 0.6$ GeV—which may be
dark matter—is most favored, while 2-body decay with $m_{X}\simeq 2$ GeV is
also competitive. The best-fit branching fractions for the scalar and tensor
cases are a few times larger than for the 2-body and vector cases. Past BaBar
measurements provide further discrimination, although the best-fit parameters
stay similar.

††preprint: KEK-TH-2587

##  I Introduction

The decay $B^{+}\to K^{+}\nu\nu$ is a rare process in the Standard Model (SM).
It has long been a sought-after experimental target Lutz  _et al._ (2013);
Grygier  _et al._ (2017); del Amo Sanchez  _et al._ (2010); Lees  _et al._
(2013) as it is one of the cleanest channels for the search for New Physics
(NP) owing to well-controlled theoretical uncertainties in the SM prediction
for its branching fraction Parrott  _et al._ (2023) (also see Bečirević  _et
al._ (2023)):

${\cal B}_{\text{SM}}(B^{+}\to K^{+}\nu\nu)=(5.58\pm 0.37)\times 10^{-6}\,.$
(1)

The Belle II collaboration recently reported Adachi  _et al._ (2023) a
measurement of this branching fraction by employing the novel Inclusive
Tagging Analysis (ITA) method Prim  _et al._ (2020); Abudinén  _et al._
(2021); Adachi  _et al._ (2023):

${\cal B}(B^{+}\to K^{+}\nu\nu)=(2.7\pm 0.5\text{(stat)}\pm
0.5\text{(syst)})\times 10^{-5}\,,$ (2)

which is obtained by scaling and fitting the normalization of a SM-like event
distribution. Ref. Adachi  _et al._ (2023) also combines this with results
from the well-established Hadronic Tagging Analysis (HTA) method, and obtains
${\cal B}(B^{+}\to K^{+}\nu\nu)=(2.3\pm
0.5\text{(stat)}^{+0.5}_{-0.4}\text{(syst)})\times 10^{-5}$. This is higher
than the SM prediction by about $2.7\sigma$, prompting a burst of theoretical
activities to look for NP via this mode Abdughani and Reyimuaji (2023); Datta
_et al._ (2023a, b); Berezhnoy and Melikhov (2023); Felkl  _et al._ (2023);
Bause  _et al._ (2023); Dreiner  _et al._ (2023); Allwicher  _et al._ (2023);
Athron  _et al._ (2023); He  _et al._ (2023); Crivellin  _et al._ (2022);
Altmannshofer  _et al._ (2023); McKeen  _et al._ (2023).

Assuming the excess is due to NP, there are two simple interpretations of the
“$\nu\nu$”, which is undetected experimentally. One (the “2-body scenario”) is
that “$\nu\nu$” is a single new particle, $X$, which is invisible or decays
invisibly experimentally Altmannshofer  _et al._ (2023); McKeen  _et al._
(2023). Another (the “3-body scenario”) is that it is a pair of particles
$\chi_{1}\chi_{2}$, where each $\chi_{i}$ may be either a new particle or a SM
(anti-)neutrino (e.g. Bause  _et al._ (2023); He  _et al._ (2023); Abdughani
and Reyimuaji (2023); Datta  _et al._ (2023a)). The 3-body scenario can be
further categorized into three sub-scenarios depending on whether the $b\to s$
transition happens via a scalar, vector, or tensor operator. The question,
therefore, is which scenario seems to be most favored by data on statistical
grounds.

In this work, we perform likelihood analyses mainly using the binned ITA data
of Ref. Adachi  _et al._ (2023), and identify the preferred ranges of the mass
and branching fraction in each scenario. Then, comparing the best-fit points
between different scenarios, we find that the vector-current 3-body scenario
with $m_{\chi_{1,2}}\simeq 0.6\>\text{GeV}$ is most statistically favored, but
the 2-body decay with $m_{X}\simeq 2$ GeV is also comparable. We also find
that in the scalar- and tensor-current 3-body scenarios the best-fit branching
fractions are a few times larger than the aforementioned value by Belle II.
This is because the scalar and tensor event distributions have a very
different shape from a SM-like (vector) event distribution assumed in the
Belle II analysis. Finally, we also explore further discrimination between
scenarios by including the past BaBar data.

##  II Testing New Physics with Belle II

###  II.1 Data Analysis

Given the latest ITA, Belle II reports the results in terms of a variable
($q_{\rm rec}^{2}$) that is only loosely related to
$q^{2}\equiv(p_{B}-p_{K})^{2}$, we carefully analyze distortions in the
reported signal shape compared to the true $q^{2}$ distributions expected in
the various scenarios. We will also find it essential to take into account the
effects of a very non-uniform signal efficiency. We further combine the
results with past measurements to narrow down the possibilities.

Since the 4-momentum of the tagged $B$ meson is not reconstructed in ITA, the
missing mass-squared, $q^{2}=(p_{B}-p_{K})^{2}$, is approximated by $q_{\rm
rec}^{2}$:

$\displaystyle q^{2}_{\text{rec}}$ $\displaystyle\equiv
E_{B}^{2}+m_{K}^{2}-2E_{B}E_{K}$
$\displaystyle=q^{2}+(E_{B}^{2}-m_{B}^{2})-2\vec{\bf{p}}_{K}(q^{2})\\!\cdot\\!\vec{\bf{p}}_{B},$
(3)

where $E_{B,K}$ and $\vec{\bf{p}}_{B,K}$ are the energies and momenta of the
observed $B$ and $K$ in the $B\bar{B}$ center-of-mass frame. The resolution in
$E_{K}$ measurement is $O(1)$%, leading to negligible smearing in $q_{\rm
rec}^{2}$. Given the bin size of $1\text{ GeV}^{2}$ chosen in Ref. Adachi  _et
al._ (2023), the smearing in $q_{\rm rec}^{2}$ almost entirely stems from the
absence of the $\vec{\bf{p}}_{K}(q^{2})\cdot\vec{\bf{p}}_{B}$ term in the
inference of $q_{\rm rec}^{2}$. By numerical simulation, we find that events
with a given value of $q^{2}$ get distributed almost uniformly in $q_{\rm
rec}^{2}$ as follows:

$\displaystyle f^{{\rm rec}}_{q^{2}}(q_{\rm
rec}^{2})=\begin{cases}(\Delta_{+}-\Delta_{-})^{-1}&\mbox{if
}\Delta_{-}<q_{\rm rec}^{2}-q^{2}<\Delta_{+}\\\
0&\mbox{otherwise}\end{cases},$ (4)
$\displaystyle\Delta_{\pm}=(E_{B}^{2}-m_{B}^{2})\pm
2\left(1\pm\frac{|\vec{\bf{p}}_{B}|}{E_{B}}\right)|\vec{\bf{p}}^{*}_{K}(q^{2})||\vec{\bf{p}}_{B}|$
(5)

where $\vec{\bf{p}}^{*}_{K}(q^{2})$ is the momentum of the $K$ in the rest
frame of the $B$ meson, $|\vec{\bf{p}}_{B}|=0.33$ GeV, and $E_{B}=5.29$ GeV.
For example, $B\to KX$ events with $m_{X}=2$ GeV would be spread over
$2.5\,{\rm GeV^{2}}<q_{\rm rec}^{2}<5.5\,{\rm GeV^{2}}$ as shown in Fig. 1.

Figure 1:  Number of $B^{+}\to K^{+}\nu\nu$ events at Belle II Adachi  _et
al._ (2023) (black dots, with error bars shown) after background subtraction,
as a function $q^{2}_{\text{rec}}$, using inclusive tagging in the
$\eta(\text{BDT}_{2})>0.98$ signal region. In blue is shown the SM
distribution Parrott  _et al._ (2023); Bečirević  _et al._ (2023). The red
line shows the predicted distribution of events for a 3-body decay $B^{+}\to
K^{+}\chi\chi$ by a vector current, with $m_{\chi}=0.6$ GeV and
$\mathcal{B}(B^{+}\to K^{+}\chi\chi)=3.2\times 10^{-5}$ in addition to the SM.
The green line shows $B^{+}\to K^{+}\nu_{L}\chi$ via a scalar current for
$m_{\chi}=0$ and $\mathcal{B}(B^{+}\to K^{+}\chi\chi)=7.3\times 10^{-5}$. The
yellow line shows the distribution for the 2-body decay $B^{+}\to K^{+}X$ with
$\mathcal{B}(B^{+}\to K^{+}X)=0.7\times 10^{-5}$ and $m_{X}=2$ GeV.  Figure 2:
Inferred efficiency in the $\eta(\text{BDT}_{2})>0.98$ signal region with the
inclusive tagging at Belle II.

Figure 3: The upper panel shows 2D plots of $\chi^{2}-\chi^{2}_{\min}$ in the
${\cal B}$-${m_{\rm NP}}$ plane for different scenarios. For each plot in the
upper panel, the likelihoods profiling one variable (along the ${\cal B}$
direction and along the $m_{\text{NP}}$ direction) are provided in the lower
panel. The three columns correspond to different NP scenarios: on the left
panel is the 2-body scenario $B^{+}\to K^{+}X$, for which $m_{\rm NP}=m_{X}$,
the center panel is for the 3-body decay $B^{+}\to K^{+}\chi\chi,\,m_{\rm
NP}=m_{\chi}$ mediated by a vector current, while the right column corresponds
to the 3-body decay $B^{+}\to K^{+}\bar{\chi}_{R}\nu_{L},m_{\rm NP}=m_{\chi}$,
where $\nu_{L}$ is a SM neutrino.

The signal efficiency, $\epsilon(q^{2})$, is another essential input to
process the various scenarios. Although two efficiencies are reported in the
latest Belle II analysis, they are not applicable to the samples with the
highest “BDT2” cut, BDT${}_{2}>0.98$, where the excess from the SM background
is visible. Therefore we deduce the efficiencies to reproduce the shape of the
reported $B\to K\nu\nu$ distribution, including the $q_{\rm rec}^{2}$ smearing
effect. The resulting efficiency, shown in Fig. 2 after interpolation in
$q^{2}$, is similar in shape to the one with BDT${}_{2}>0.95$ reported in the
earlier Belle II analysis Abudinén  _et al._ (2021).

For the statistical combination, we construct the binned likelihood,
$L=\prod_{i}f_{\rm P}(N_{i}^{{\rm obs}};N_{i}^{\rm ex})$, using Poisson
statistics for the $i$th bin $q^{2}_{i}\leq q_{\rm rec}^{2}<q^{2}_{i}+1\,{\rm
GeV^{2}}$. The expectation $N_{i}^{\rm ex}=N_{i}^{\rm BG}+N_{i}^{K\nu\nu,\rm
SM}+N_{i}^{\rm NP}$ includes the background and the SM $B\to K\nu\nu$ with as
well as the yield of new physics,

$\displaystyle N_{i}^{\rm
NP}=N_{B^{\pm}}\int_{q^{2}_{i}}^{q^{2}_{i+1}}\mathrm{d}x\int\mathrm{d}q^{2}\,f^{{\rm
rec}}_{q^{2}}(x)\,\epsilon(q^{2})\frac{\mathrm{d}{{\cal B}_{\rm
NP}}}{\mathrm{d}q^{2}}\,,$ (6)

where $N_{B^{\pm}}=3.99\times 10^{8}$ at Belle II. We obtain the log-
likelihood, $-2\ln L$, which is denoted as $\chi^{2}$ for simplicity. This
method accounts for only the statistical uncertainty. As a check, we obtain
${\cal B}=2.6\pm 0.4$(stat) for $B^{+}\to K^{+}\nu\nu$ for ITA, which is
reasonably consistent with ${\cal B}=2.7\pm 0.5$(stat) reported by Belle II.

Figure 4:  Normalized distributions for $B^{+}\to K^{+}\chi_{1}\chi_{2}$ via
scalar, vector, and tensor operators as functions of $q^{2}$ for massless
$\chi_{1,2}$ (solid lines) as well as massive $\chi_{1,2}$ (dotted line) with
an equal mass of $0.6\,\text{GeV}$. The distributions do not take experimental
efficiency into account.

###  II.2 Testing 2-body scenario

First, we calculate the log-likelihood $\chi^{2}$, based on the 2-body decay
$B^{+}\to K^{+}X$ scenario. We find the minimum $\chi^{2}_{\rm min}$ and
obtain the preferred parameter space by evaluating $\chi^{2}-\chi^{2}_{\rm
min}$. Since there is a significant mass preference, we perform the 2D fit in
the $m_{X}$-$\cal B$ space, and the contours with 1, 2, and 3$\sigma$ are
shown in the upper-left pane in Fig. 3.

Although our method correctly accounts for the statistical uncertainty, the
systematic uncertainty can be comparable based on the Belle II report. We
expect that including systematic uncertainty lowers the $\chi^{2}$ by a factor
of two, resulting in the allowed region being enlarged by $\sim 40\%$, but the
tendency of the preferred parameters would not be affected.

To extract the preferred mass, we use the profile log-likelihood: finding the
best $\chi^{2}$ for each mass by profiling $\cal B$, which is shown in the
bottom left panel. We find a narrow range of mass, $m_{X}=1.97\pm 0.05$(stat)
GeV. The example of distribution in the ITA is shown in Fig. 1. A similar
method is applied for the branching fraction, leading to ${\cal
B}=\left(0.79\begin{subarray}{c}+0.16\\\
-0.13\end{subarray}{\rm(stat)}\right)\times 10^{-5}$. This is in accordance
with a recent study of the 2-body scenario Altmannshofer  _et al._ (2023), but
it does not agree with the other McKeen  _et al._ (2023). We wish to emphasize
that the expected number of events in the given $q_{\rm rec}^{2}$ binning is
calculated following Eq. (5) so that we appropriately fit the signal with the
data of the ITA.

###  II.3 Testing 3-body scenarios

The same method is applied to the 3-body case $B^{+}\to K^{+}\chi_{1}\chi_{2}$
where we have multiple scenarios to consider. We assume that $\chi_{1,2}$ are
fermionic in this work. However, some operators that contribute to $B\to K$
decays also contribute to $B\to K^{*}$ decays. Belle-II has recently obtained
an upper limit on the branching fraction ${\cal B}(B^{0}\to
K^{*}\nu\nu)<1.8\times 10^{-5}$ Grygier  _et al._ (2017), but a measurement of
the branching fraction has not been obtained yet. We consider operators that
do not have significant effects on $B\to K^{*}$:

  1. 1.

Scalar operator: the interaction is given by

${\cal L}\supset{1\over\Lambda_{S}^{2}}(\bar{b}s)(\bar{\chi}_{1}\chi_{2}),$
(7)

where $\Lambda_{S}$ is a heavy scale. We strictly distinguish this case from
the pseudo-scalar case with $\bar{b}\gamma_{5}s$ as that would contribute only
to $B\to K^{*}$ and not to $B\to K$, while the scalar case contributes only to
$B\to K$ and not to $B\to K^{*}$. In our work, we consider the case where
$\chi_{1}=\nu_{L}$, which is massless for our purpose, and $\chi_{2}$ is a new
fermion, possibly massive. Note that in Ref. Deppisch  _et al._ (2020) it was
shown that a scalar current distribution for the related kaon decay mode
$K^{+}\to\pi^{+}\nu\nu$ implies lepton number violation (LNV) for SM-invariant
operators, since such a current can only be generated via $\Delta L=2$ odd
mass-dimension operators. We consider the possibility of $\chi_{2}$ being the
right-handed neutrino $N$, for which case there exists non-LNV operators at
dimension-6 which generate a scalar current, e.g. the operator
$\mathcal{O}_{LNQd}=\epsilon_{ij}\overline{L^{i}}N\overline{Q^{j}}d$ Liao and
Ma (2017).

  2. 2.

Vector operator: the interaction is given by

${\cal
L}\supset{1\over\Lambda_{V}^{2}}(\bar{b}\gamma^{\mu}s)(\bar{\chi}_{1}\gamma_{\mu}\chi_{2})\,,$
(8)

where $\Lambda_{V}$ is a heavy scale. $\chi_{1}$ and $\chi_{2}$ may be a pair
of SM $\nu_{L}\nu_{L}$ or that of new (massive) neutral fermions. The quark
part of this operator is a pure vector, which is crucially different from a
SM-like $V$$-$$A$ current. As shown in Fig. 2 of Bause  _et al._ (2023), an
additional $V$$-$$A$ contribution that fits the $B\to K$ excess would be
severely excluded by the absence of a corresponding excess in $B\to K^{*}$,
while the pure vector case has the least impact on $B\to K^{*}$ and is still
allowed by data. Therefore, for our interest in identifying the best
scenarios, we do not consider any axial-vector component in the quark
bilinear.

  3. 3.

Tensor operator: the interaction is given by

${\cal
L}\supset{1\over\Lambda_{T}^{2}}(\bar{b}\sigma^{\mu\nu}s)(\bar{\chi}_{1}\sigma_{\mu\nu}\chi_{2})\,,$
(9)

where $\Lambda_{T}$ is a heavy scale. As in the scalar case above, we consider
$\chi_{1}=\nu_{L}$ and $\chi_{2}$ to be a new neutral fermion.

In each of these three cases, we compute the $q^{2}$ spectrum of the decay
rate $\mathrm{d}\Gamma/\mathrm{d}q^{2}$, and the distributions are given in
Fig. 4. The relevant matrix elements and form factors for the $B\to K$ decay
are given in Isgur and Wise (1990); Ball and Zwicky (2005) and also in
Appendix A. In the ITA, the distributions are modified due to the efficiency
and smeared by $q_{\rm rec}^{2}$. Two representative cases are shown in Fig.
1.

We find the favored parameter space of the 3-body scenarios using the same
likelihood analysis as we perform for the 2-body one, and the results for the
cases of the vector and scalar operators are shown in Fig. 3. For the vector
case, the fit is best at $m_{\chi}=0.62\begin{subarray}{c}+0.10\\\
-0.09\end{subarray}\,\text{GeV}$ and ${\cal B}=(3.5\begin{subarray}{c}+0.8\\\
-0.6\end{subarray}(\text{stat}))\times 10^{-5}$ (see the middle panel), which
is close to the measurement made at Belle II. Since this massive $\chi$ be can
be stable, it is a candidate of dark matter.

For the scalar operator, the best-fit mass is $m_{\chi}=0$, while the
branching fraction is preferred to be high, $(7.9\begin{subarray}{c}+1.6\\\
-1.5\end{subarray}(\text{stat})\times 10^{-5})$ (see the right panel). This is
due to that the $q^{2}$-spectrum peaking at high $q^{2}\sim 20\
\text{GeV}^{2}$, as seen in Fig. 4, is multiplied by the efficiency which
drops at high $q^{2}$, see Fig. 2. Therefore, most of the signal events are
discarded by the ITA, which statistically favors a higher branching fraction
to explain the excess.

For the tensor case, since the $q^{2}$-distribution is very similar to the
scalar case, as seen in Fig. 4, the fitted results are also similar, and the
preferred branching fraction is ${\cal B}=(6.4\pm 1.2(\text{stat}))\times
10^{-5}$. We do not show the tensor case in Fig. 3 for the sake of brevity,
but we quote $\chi^{2}_{\text{min}}$ for the fit to the decay $B\to K\chi\nu$
in Table 1.

For the best-fit points in Fig. 3, we find $\Lambda_{V}\approx 6.5$ TeV for
the 3-body vector current and $\Lambda_{S}\approx 4.1$ TeV for the 3-body
scalar current. The 2-body decay occurs at tree-level, for which we find an
effective coupling $\lambda_{X}\approx 3\times 10^{-4}$ for
$\mathcal{L}\supset-\lambda_{X}\bar{b}sX$.

As already addressed, our analysis does not include systematic uncertainty.
The systematic uncertainty would enlarge the allowed parameter space by about
$40\%$ unless the uncertainty has a strong $q^{2}$ dependence.

###  II.4 Comparing All Scenarios

Based on the values $\chi^{2}_{\min}$ of different scenarios, shown in Table
1, we can distinguish the preferred ones out of them. Testing one hypothesis
$H_{1}$ against the null hypothesis $H_{0}$ can be evaluated by Workman and
Others (2022)

$\displaystyle\alpha=\frac{\max L(H_{1})}{\max L(H_{0})}=\exp\left(-[\chi_{\rm
min}^{2}(H_{1})-\chi_{\rm min}^{2}(H_{0})]/{2}\right)\,.$ (10)

Within the Belle II ITA, we find that the best scenario is the 3-body decay
with vector-operator and $m_{\chi}=0.62$ GeV, and the best 2-body scenario is
also competitive. However, the scalar-current 3-body scenario is disfavored,
and adding mass to $\chi$ makes it even worse. A SM-like shape, i.e. the
vector-current 3-body scenario with massless $\chi$, is also disfavored. Even
if unaccounted uncertainties, such as the systematics, lower $\chi^{2}$, say
by a factor of two, the scalar and vector operators with massless $\chi$ would
still have $\alpha<10\%$ against the best-fit scenario.

|  $\chi^{2}_{\rm min}-100$ |  |  2b |  V |  V′ |  S |  T |  SM
---|---|---|---|---|---|---|---|---
|  Belle II |  |  6.8 |  15.2 |  4.7 |  15.1 |  11.9 |  44.6
|  \+ BaBar SR |  |  27.6 |  30.4 |  22.1 |  31.8 |  29.8 |  61.0
|  \+ BaBar $s_{B}<0.8$ |  |  73.3 |  78.8 |  72.9 |  90.2 |  86.9 |  106.7
Table 1: The likelihood minima $\chi^{2}_{\text{min}}$ for 2-body scenario
(2b), and several 3-body scenarios: two vector cases with $m_{\chi}=0$ (V) and
$m_{\chi}\simeq 0.6\,\text{GeV}$ (V′); the scalar (S) and tensor (T) cases.
The rows correspond to the choice of analyzed datasets: first only the Belle
II ITA data Adachi  _et al._ (2023), then adding the BaBar data Lees  _et al._
(2013) within the signal region (SR), and finally including outside the SR,
$q^{2}<0.8m_{B}^{2}$. The bold ones indicate the best-fit or competitive
scenarios with $\alpha(H_{0}=\rm V^{\prime})>10\%$. The $\chi^{2}$ values for
the SM are also shown. Only statistical uncertainties are considered.

##  III Combining with Other $B\to K\nu\nu$ Measurements

Having considered the Belle II ITA, we now combine the analysis with past
BaBar measurements of $B\to h\nu\nu$ with $h=K^{0},K^{+}$ based on the
conventional hadronic tagging Lees  _et al._ (2013). The shape with the BaBar
data (after the background is subtracted) with several new physics scenarios
in Fig. 5. A major difference from the latest Belle II ITA is that the
analysis uses $q^{2}$ rather than $q_{\rm rec}^{2}$, and the resolution of
$q^{2}$ is significantly better than the adopted bin size. Therefore, we can
ignore smearing effects.

For their signal region $s_{B}\equiv q^{2}/m_{B}^{2}<0.3$, we calculate the
log-likelihood $\chi^{2}$ for this data, which we then combine with the result
of Belle II ITA. The combined fits are shown as magenta dashed lines in Fig.
3, and the values of $\chi^{2}_{\min}$ are in Table 1. While the preferred
parameter regions barely change in all scenarios, the difference of
$\chi^{2}_{\min}$ between the 2-body scenario and the best 3-body scenario
increases from 2 to 5 because the signal of 2-body decay should be localized
in the second bin where no excess is seen.

Although $s_{B}>0.3~{}(q^{2}\gtrsim 8.4\,{\rm GeV}^{2})$ is outside the signal
region and requires more careful treatment of the uncertainties, it is
interesting to examine this region for further discrimination. For some 3-body
scenarios with the scalar and tensor operators, we expect an excess in high
$q^{2}$ because the efficiency in this BaBar study stays flat until $q^{2}\sim
20\,{\rm GeV}^{2}$ while the efficiency of the Belle II ITA is very small for
$q^{2}>10\,{\rm GeV}^{2}$. Thus, increasing the branching fraction to fit the
Belle II ITA would be in tension with $B^{+}\to K^{+}\nu\nu$ result at BaBar,
as seen in Fig. 5. For this reason, the $\chi^{2}_{\min}$ values of these
scenarios are significantly worsened, see Table 1. The best scenario is still
the 3-body decay with massive $\chi$ via the vector operator, but the 2-body
one is comparable.

Figure 5:  Number of $B\to K\nu\nu$ events (black, error bars shown) for
$471\times 10^{6}$ $B\bar{B}$ pairs at BaBar Lees  _et al._ (2013) with
background subtracted, as a function of $q^{2}$. The colored lines show the
same benchmark scenarios as Fig. 1.

Other data sets that may be relevant for us, but are not included, are the
following:

  1. 1.

In Adachi  _et al._ (2023), Belle II also reported the $q^{2}$ result in the
HTA, see Fig. 6. Since the data has a small overlap with the ITA dataset, we
are unable to combine them. However, we see the trend that the third bin seems
to have a tension with the 2-body scenario with $m_{X}=2$ GeV.

  2. 2.

Belle results via hadronic Lutz  _et al._ (2013) and semileptonic Grygier  _et
al._ (2017) tagging: the number of events is measured in bins of $E_{\rm
ECL}$, which is the residual energy deposited in the electromagnetic
calorimeter (ECL). Since $E_{\rm ECL}$ is not related to $q^{2}$, we are
unable to include this data for testing the various scenarios.

  3. 3.

BaBar results via semileptonic tagging del Amo Sanchez  _et al._ (2010): The
binning is done in bins of $|p_{K}|$, the magnitude of the momentum of the
$K^{+}$ in the center-of-mass frame. While $|p_{K}|$ has a one-to-one
correspondence to $q_{\rm rec}^{2}$, the final results are reported as an
average of 20 different BDT analyses. Due to this unconventional statistical
treatment, we cannot safely add this data to our analysis.

Figure 6:  Number of $B^{+}\to K^{+}\nu\nu$ events at Belle II Adachi  _et
al._ (2023) (black dots, with error bars shown) after background subtraction
using hadronic tagging, as a function of $q^{2}$. The colored lines show the
same benchmark scenarios as in Fig. 1.

##  IV Conclusion

In this work, we perform a binned-likelihood analysis of different NP
scenarios mainly based on the recent Belle II ITA data on the decay $B^{+}\to
K^{+}\nu\nu$. The NP scenarios we consider are: (a) 2-body scenario $B\to KX$,
and (b) the 3-body scenario $B\to K\chi_{1}\chi_{2}$. For (b) we consider
several operators (scalar, vector, and tensor) as well as several masses for
the new particles.

We find it crucial to account for the fact that the Belle II ITA data is
binned not in the momentum transfer $q^{2}$ but in $q_{\rm rec}^{2}$, and for
the non-uniform efficiency. We augment our analysis with past BaBar data. Our
results are listed in Table 1 and in Fig. 3.

As seen in Table 1, the best fit NP scenario is a 3-body decay of the $B^{+}$
via a vector current into $K^{+}$ and a $\chi\chi$ pair, where $m_{\chi}\simeq
0.6\,\text{GeV}$. This $\chi$ can be stable and hence a possible dark matter
candidate. Table 1 also tells us that the 2-body scenario is not much worse,
and so it still remains competitive.

The right panel of Fig. 3, shows the preferred parameter region in the scalar-
current 3-body scenario. Massless $\chi$ fits best, for a branching fraction
of ${\cal B}(B^{+}\to K\chi_{R}\nu_{L})\sim 8\times 10^{-5}$. This is about a
factor of 3 higher than the value in Eq. (2). This is an especially clear
demonstration of the importance of considering the non-uniform signal
efficiency in inferring the branching fraction of NP from the reported excess
at Belle II.

###### Acknowledgements.

We thank Eldar Ganiev for useful discussions. This work is supported in part
by the US Department of Energy grant DE-SC0010102. The work of K.F., T.O., and
K.T. is supported in part by JSPS Grant-in-Aid for Scientific Research (Grant
No. 21H01086). K.F. acknowledges support from the JSPS KAKENHI grant
JP22K21350.

##  Appendix A Hadronic Matrix Elements and Form factors for $B\to K$

The hadronic matrix elements of some of the operators relevant to $B\to K$
decays are given below Isgur and Wise (1990); Ball and Zwicky (2005):

$\displaystyle\langle{K}|{\bar{b}s}|{B}\rangle$ $\displaystyle=$
$\displaystyle{m_{B}^{2}-m_{K}^{2}\over m_{s}-m_{b}}f_{0}(q^{2})\,,$ (11)
$\displaystyle\langle{K}|{\bar{b}\gamma^{\mu}s}|{B}\rangle$ $\displaystyle=$
$\displaystyle f_{+}(q^{2})(p_{B}+p_{K})^{\mu}$
$\displaystyle+\bigl{[}f_{0}(q^{2})-f_{+}(q^{2})\bigr{]}{m_{B}^{2}-m_{K}^{2}\over
q^{2}}q^{\mu}\,,$
$\displaystyle\langle{K}|{\bar{b}\sigma^{\mu\nu}s}|{B}\rangle$
$\displaystyle=$ $\displaystyle i{f_{T}(q^{2})\over
m_{B}+m_{K}}\bigl{[}(p_{B}+p_{K})^{\mu}q^{\nu}$ (13) $\displaystyle\hskip
55.97205pt-(p_{B}+p_{K})^{\nu}q^{\mu}\bigr{]}\,,$

where $q^{\mu}=p_{B}^{\mu}-p_{K}^{\mu}$ and the form factors $f_{0},f_{+}$ and
$f_{T}$ can be given approximate analytical expressions via lattice methods.
Here we quote the results obtained in Ref. Ball and Zwicky (2005):

$\displaystyle f_{0}(q^{2})$ $\displaystyle=$
$\displaystyle{r_{1}^{P}+r_{2}^{P}\over 1-q^{2}/m_{\text{fit}}^{2}}\,,$ (14)
$\displaystyle f_{+}(q^{2})$ $\displaystyle=$ $\displaystyle{r_{1}^{P}\over
1-q^{2}/M_{B}^{2}}+{r_{2}^{P}\over(1-q^{2}/M_{B}^{2})^{2}}\,,$ (15)
$\displaystyle f_{T}(q^{2})$ $\displaystyle=$ $\displaystyle{r_{1}^{T}\over
1-q^{2}/M_{B}^{2}}+{r_{2}^{T}\over(1-q^{2}/M_{B}^{2})^{2}}\,,$ (16)

where $r_{1}^{P}=0.162$, $r_{2}^{P}=0.173$, $r_{1}^{T}=0.161$,
$r_{2}^{T}=0.198$, $M_{B}=5.41\;\text{GeV}$ and
$m^{2}_{\text{fit}}=37.46\;\text{GeV}^{2}$ are fit parameters in the lattice.

Note that using isospin symmetry, we claim that the form factors for the
charged $B\to K$ modes are essentially identical. Using these relations we can
find expressions for the branching fraction for NP contributions to $B^{+}\to
K^{+}+$ invisible. For example, the branching fraction for an additional
vector current contribution to $B^{+}\to K^{+}\nu\nu$ is given by

${\cal B}_{\text{NP}}(B^{+}\to
K^{+}\nu\nu)\times\Gamma_{B^{+}}^{\text{tot}}=\int_{0}^{(q^{2})^{+}}\mathrm{d}q^{2}\int_{t^{-}}^{t^{+}}\mathrm{d}t\,\frac{\mathrm{d}\Gamma_{V}}{\mathrm{d}q^{2}\mathrm{d}t}$
(17)

where $(q^{2})^{+}=(p_{B^{+}}-p_{K^{+}})^{2}$ and

$\displaystyle t^{\pm}$
$\displaystyle=\frac{(m_{B^{+}}^{2}-m_{K^{+}}^{2})^{2}}{4q^{2}}$ (18)
$\displaystyle\qquad-\frac{1}{4q^{2}}\left(\lambda^{1/2}({q^{2},m_{B^{+}}^{2},m_{K^{+}}^{2}})\mp
q^{2}\right)^{2}\,,$

where $\lambda(a,b,c)$ is the Källén function. Furthermore, we have

$\displaystyle\frac{\mathrm{d}\Gamma_{V}}{\mathrm{d}q^{2}\mathrm{d}t}=$
$\displaystyle\frac{1}{(2\pi)^{3}}\frac{1}{32m_{B^{+}}^{3}}\left(\frac{1}{\Lambda_{V}^{2}}\right)^{2}$
(19)
$\displaystyle\times|\langle{K}|{\bar{b}\gamma^{\mu}s}|{B}\rangle\langle{\nu\bar{\nu}}|{\bar{\nu}\gamma_{\mu}\nu}|{0}\rangle|^{2}\,,$

where $\Gamma_{B^{+}}^{\text{tot}}\approx 4.0\times 10^{-13}$ GeV is the total
width of the $B^{+}$.

## References

  * Lutz  _et al._ (2013) O. Lutz _et al._ (Belle), Phys. Rev. D 87, 111103 (2013), arXiv:1303.3719 [hep-ex] .
  * Grygier  _et al._ (2017) J. Grygier _et al._ (Belle), Phys. Rev. D 96, 091101 (2017), [Addendum: Phys.Rev.D 97, 099902 (2018)], arXiv:1702.03224 [hep-ex] .
  * del Amo Sanchez  _et al._ (2010) P. del Amo Sanchez _et al._ (BaBar), Phys. Rev. D 82, 112002 (2010), arXiv:1009.1529 [hep-ex] .
  * Lees  _et al._ (2013) J. P. Lees _et al._ (BaBar), Phys. Rev. D 87, 112005 (2013), arXiv:1303.7465 [hep-ex] .
  * Parrott  _et al._ (2023) W. G. Parrott, C. Bouchard, and C. T. H. Davies (HPQCD), Phys. Rev. D 107, 014511 (2023), [Erratum: Phys.Rev.D 107, 119903 (2023)], arXiv:2207.13371 [hep-ph] .
  * Bečirević  _et al._ (2023) D. Bečirević, G. Piazza, and O. Sumensari, Eur. Phys. J. C 83, 252 (2023), arXiv:2301.06990 [hep-ph] .
  * Adachi  _et al._ (2023) I. Adachi _et al._ (Belle-II), (2023), arXiv:2311.14647 [hep-ex] .
  * Prim  _et al._ (2020) M. T. Prim _et al._ (Belle), Phys. Rev. D 101, 032007 (2020), arXiv:1911.03186 [hep-ex] .
  * Abudinén  _et al._ (2021) F. Abudinén _et al._ (Belle-II), Phys. Rev. Lett. 127, 181802 (2021), arXiv:2104.12624 [hep-ex] .
  * Abdughani and Reyimuaji (2023) M. Abdughani and Y. Reyimuaji, (2023), arXiv:2309.03706 [hep-ph] .
  * Datta  _et al._ (2023a) A. Datta, A. Hammad, D. Marfatia, L. Mukherjee, and A. Rashed, JHEP 03, 108 (2023a), arXiv:2210.15662 [hep-ph] .
  * Datta  _et al._ (2023b) A. Datta, D. Marfatia, and L. Mukherjee, (2023b), arXiv:2310.15136 [hep-ph] .
  * Berezhnoy and Melikhov (2023) A. Berezhnoy and D. Melikhov, (2023), arXiv:2309.17191 [hep-ph] .
  * Felkl  _et al._ (2023) T. Felkl, A. Giri, R. Mohanta, and M. A. Schmidt, (2023), arXiv:2309.02940 [hep-ph] .
  * Bause  _et al._ (2023) R. Bause, H. Gisbert, and G. Hiller, (2023), arXiv:2309.00075 [hep-ph] .
  * Dreiner  _et al._ (2023) H. K. Dreiner, J. Y. Günther, and Z. S. Wang, (2023), arXiv:2309.03727 [hep-ph] .
  * Allwicher  _et al._ (2023) L. Allwicher, D. Becirevic, G. Piazza, S. Rosauro-Alcaraz, and O. Sumensari, (2023), arXiv:2309.02246 [hep-ph] .
  * Athron  _et al._ (2023) P. Athron, R. Martinez, and C. Sierra, (2023), arXiv:2308.13426 [hep-ph] .
  * He  _et al._ (2023) X.-G. He, X.-D. Ma, and G. Valencia, (2023), arXiv:2309.12741 [hep-ph] .
  * Crivellin  _et al._ (2022) A. Crivellin, C. A. Manzari, W. Altmannshofer, G. Inguglia, P. Feichtinger, and J. Martin Camalich, Phys. Rev. D 106, L031703 (2022), arXiv:2202.12900 [hep-ph] .
  * Altmannshofer  _et al._ (2023) W. Altmannshofer, A. Crivellin, H. Haigh, G. Inguglia, and J. Martin Camalich, (2023), arXiv:2311.14629 [hep-ph] .
  * McKeen  _et al._ (2023) D. McKeen, J. N. Ng, and D. Tuckler, (2023), arXiv:2312.00982 [hep-ph] .
  * Deppisch  _et al._ (2020) F. F. Deppisch, K. Fridell, and J. Harz, JHEP 12, 186 (2020), arXiv:2009.04494 [hep-ph] .
  * Liao and Ma (2017) Y. Liao and X.-D. Ma, Phys. Rev. D 96, 015012 (2017), arXiv:1612.04527 [hep-ph] .
  * Isgur and Wise (1990) N. Isgur and M. B. Wise, Phys. Rev. D 42, 2388 (1990).
  * Ball and Zwicky (2005) P. Ball and R. Zwicky, Phys. Rev. D 71, 014015 (2005), arXiv:hep-ph/0406232 .
  * Workman and Others (2022) R. L. Workman and Others (Particle Data Group), PTEP 2022, 083C01 (2022).