# Signatures of Light New Particles in $B\to K^{(*)}E_{\rm miss}$

Patrick D. Bolton<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39,
1000 Ljubljana, Slovenia  Svjetlana Fajfer<EMAIL_ADDRESS>Jožef
Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia  Faculty of Mathematics
and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia
Jernej F. Kamenik<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39,
1000 Ljubljana, Slovenia  Faculty of Mathematics and Physics, University of
Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia  Martín Novoa-Brunet
<EMAIL_ADDRESS>Instituto de Física Corpuscular, Universitat de
Valéncia – Consejo Superior de Investigaciones Científicas, Parc Científic,
E-46980 Paterna, Valencia, Spain  Istituto Nazionale di Fisica Nucleare,
Sezione di Bari, Via Orabona 4, 70126 Bari, Italy

###### Abstract

The recent Belle II observation of $B\to KE_{\rm miss}$ challenges theoretical
interpretations in terms of Standard Model neutrino final states. Instead, we
consider new physics scenarios where up to two new light-invisible particles
of spin 0 up to 3/2 are present in the final state. We identify viable
scenarios by reconstructing the (binned) likelihoods of the relevant $B\to
K^{(*)}E_{\rm miss}$ and also $B_{s}\to E_{\rm miss}$ experimental analyses
and present preferred regions of couplings and masses. In particular, we find
that the current data prefers two-body decay kinematics involving the emission
of a single massive scalar or a vector particle, or alternatively, three-body
decays involving pairs of massive scalars or spin 1/2 fermions. When
applicable, we compare our findings with existing literature and briefly
discuss some model-building implications.

††preprint: BARI-TH/757-24

##  I Introduction

Historically, rare decays of (heavy) flavoured mesons have been important
probes and harbingers of new physics (NP). In the last decade, the LHCb and
B-factory experiments have produced several intriguing results on rare
semileptonic $B$ meson decays, including the charged current mediated $B\to
D^{(*)}\tau\nu$ and the flavour-changing neutral current (FCNC) mediated $B\to
K^{(*)}\mu^{+}\mu^{-}$ (see, e.g., Refs. Amhis  _et al._ (2023); Algueró  _et
al._ (2023)), that challenge explanations within the Standard Model (SM).

Most recently, the Belle II experiment has measured the branching ratio
$\mathcal{B}(B\to KE_{\text{miss}})=(2.3\pm 0.7)\times 10^{-5}$. Assuming that
the missing energy $E_{\text{miss}}$ is carried away by a pair of massless SM
neutrinos, the result lies $2.9\sigma$ above the SM prediction for
$\mathcal{B}(B\to K\nu\bar{\nu})$ Adachi  _et al._ (2023). The quark-level
transition $b\to s\nu\bar{\nu}$ is also probed by the complementary mode $B\to
K^{\ast}E_{\text{\rm miss}}$, with current experimental upper bounds on the
branching ratio of the order $\mathcal{B}(B\to K^{*}E_{\text{\rm
miss}})<11\times 10^{-5}$ Lees  _et al._ (2013). Finally, an upper limit on
the branching fraction of invisible $B_{s}$ decays, $\mathcal{B}(B_{s}\to
E_{\rm miss})<5.4\times 10^{-4}$ (90% CL), has been recently derived in Ref.
Alonso-Álvarez and Escudero (2023) using a recast of ALEPH data Barate  _et
al._ (2001). While innocuous in the SM, we show that this mode can put
competitive bounds on particular NP interpretations of the Belle II result.

Undetected particles (neutrinos) in the final state make these FCNC processes
experimentally more challenging compared to those producing charged leptons.
On the other hand, they are theoretically cleaner Buras (2020); Buras  _et
al._ (2015). In particular, the relevant hadronic matrix elements are well
understood within existing theoretical frameworks Bečirević  _et al._ (2023);
Gubernari  _et al._ (2023); Athron  _et al._ (2023).

A variety of NP models aiming to resolve the Belle II excess have been
proposed in the literature Bause  _et al._ (2024); Allwicher  _et al._ (2024);
Athron  _et al._ (2023); Chen and Chiang (2023); Felkl  _et al._ (2023);
Abdughani and Reyimuaji (2023); Dreiner  _et al._ (2023); He  _et al._ (2023);
Berezhnoy and Melikhov (2023); Datta  _et al._ (2023); Altmannshofer  _et al._
(2023); Fridell  _et al._ (2023); Gabrielli  _et al._ (2024); Chen  _et al._
(2024); Hou  _et al._ (2024). In some scenarios, the SM neutrinos still carry
away all of the observed missing energy Bause  _et al._ (2024); Allwicher  _et
al._ (2024); Chen and Chiang (2023); Athron  _et al._ (2023), while in others,
novel undetected decay products are also present in the final state
Altmannshofer  _et al._ (2023); Fridell  _et al._ (2023); Gabrielli  _et al._
(2024). Restricting the outgoing invisible states to SM neutrinos only, there
are two possibilities for the associated NP; either it couples universally to
all three lepton generations Athron  _et al._ (2023); Allwicher  _et al._
(2024), or it prefers some (e.g., $\nu_{\tau}$) neutrino flavours Bause  _et
al._ (2024); Allwicher  _et al._ (2024), leading to the violation of lepton
flavour universality (LFU). A general analysis of LFU and non-LFU NP, coupling
to both left- and right-handed quark current operators and contributing to the
decays $B\to K^{(\ast)}\ell^{+}\ell^{-}$ and $B\to K^{(\ast)}\nu\bar{\nu}$,
was presented in Ref. Descotes-Genon  _et al._ (2020). It showed that
enhancing the $B\to K\nu\bar{\nu}$ branching ratio while simultaneously
satisfying existing constraints on $\mathcal{B}(B\to K^{*}\nu\bar{\nu})$
necessitates a significant contribution from right-handed quark current
operators. In the LFU limit, however, these operators induce large
contributions to the $C_{9,10}^{\prime}$ Wilson coefficients, present in the
$b\to s\ell^{+}\ell^{-}$ effective Hamiltonian, which are already ruled out by
existing measurements of $B\to K^{(*)}\mu^{+}\mu^{-}$ and
$B_{s}\to\mu^{+}\mu^{-}$ decays Algueró  _et al._ (2023). Thus, the only
remaining phenomenologically viable option is LFU-violating NP, which couples
both left- and right-handed quark currents predominantly to $\tau$ neutrinos
Bause  _et al._ (2024); Allwicher  _et al._ (2024).

Motivated by the remarkable Belle II measurement of $B\to KE_{\text{\rm
miss}}$, but also the phenomenological difficulties in accommodating the
observed excess exclusively with SM neutrinos, we turn instead to the
interesting scenario that additional light invisible states are present in the
final state. We systematically consider both single scalar and vector particle
final states, as well as pairs of scalars, spin $1/2$ and $3/2$ fermions, and
vectors, following Ref. Kamenik and Smith (2012). Since several of these
possibilities significantly alter the phase space and kinematic distributions
of events in the experiments, we consider not only the total branching ratios
$\mathcal{B}(B\to K^{(*)}E_{\text{miss}})$, but also all available
distributions presented in the Belle II and BaBar analyses. Finally,
interactions producing two invisible particles in the final state of $B\to
KE_{\text{miss}}$ can induce the invisible $B_{s}$ decay $B_{s}\to
E_{\text{miss}}$ whose larger phase space can probe massive invisible states
beyond the kinematical limit of $B\to K^{*}E_{\text{miss}}$. The experimental
upper bound from ALEPH Alonso-Álvarez and Escudero (2023); Barate  _et al._
(2001) can thus be considered as complementary.

Using this data, we construct our likelihoods for the possible invisible final
states, discerning which scenarios are favoured (and, if so, what masses and
couplings are implied).

We note that partial analyses of some of the scenarios considered in this work
have already been performed in the literature Felkl  _et al._ (2023);
Abdughani and Reyimuaji (2023); Dreiner  _et al._ (2023); He  _et al._ (2023);
Berezhnoy and Melikhov (2023); Datta  _et al._ (2023); Altmannshofer  _et al._
(2023); Fridell  _et al._ (2023). However, a comprehensive study of all
possibilities taking into account all available decay distributions has not
been implemented to date. When applicable, we compare our results and findings
with the previously published results.

Figure 1: (Above) The process $B^{+}\to K^{(*)+}\sum X$, resulting in the
$B\to KE_{\text{miss}}$ signal. In the SM, $\sum X=\nu\bar{\nu}$, but
additional light invisible particles (for example, a pair of heavy Majorana
fermions, $\sum X=N\bar{N}$) may be present in the final state. (Below) The
same invisible final states with the same couplings induce the decay
$B_{s}\to\sum X$, leading to the $B_{s}\to E_{\rm miss}$ signal.

The remainder of this work is structured as follows. In Section II, we
introduce the NP fields describing the invisible final states, outlining their
different (effective) couplings to the quark currents relevant to $B$ meson
decays. We describe some features of the kinematic distributions of the
decays, while explicit expressions for the branching ratios are given in
Appendix B. In Section III, we give an overview of how Belle II and BaBar data
are used to construct likelihoods for the different NP scenarios.

In Section IV, we then discuss the results of minimising these (negative) log-
likelihoods as a function of the invisible particle masses and couplings.
Firstly, we show that certain scenarios are immediately disfavoured. Of the
scenarios that remain, we explore what masses and couplings are implied by the
Belle II excess and are at the same time compatible with existing constraints
from BaBar and ALEPH. We also compare our findings to existing results in the
literature and briefly discuss some model-building implications before
concluding in Section V.

##  II Model Considerations

As depicted in Fig. 1, we consider the contribution of additional invisible
final states, denoted as $\sum X$, to $B\to K^{(*)}E_{\text{miss}}$. This is
alongside the SM neutrinos, $\sum
X=\sum_{\alpha}\nu_{\alpha}\bar{\nu}_{\alpha}\equiv\nu\bar{\nu}$, with
$\alpha\in\\{e,\mu,\tau\\}$. Any number of invisible final states may be
present; however, to avoid phase space suppression we will consider only one
or two invisible final state particles. With this requirement, there are now
only a few possible scenarios. Considering fields
$X\in\\{\phi,\psi,V_{\mu},\Psi_{\mu}\\}$, corresponding to (massive) particles
of spin $J=\\{0,1/2,1,3/2\\}$, respectively, leads to the following possible
final states:

$\displaystyle\sum
X\in\\{\phi,V,\phi\bar{\phi},\psi\bar{\psi},V\bar{V},\Psi\bar{\Psi}\\}\,.$ (1)

Thus, two-body decays can only involve scalar and vector bosons, $B\to
K^{(*)}\phi/V$. Three-body decays can proceed to pairs of scalars, vectors,
and also spin $1/2$ and $3/2$ fermions, $B\to
K^{(*)}\phi\bar{\phi}/\psi\bar{\psi}/V\bar{V}/\Psi\bar{\Psi}$.

We assume that the invisible state is a singlet under the SM gauge group
$SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$, but leave open the possibility
that it is charged under a dark gauge or global symmetry. As a result, any
observable effect of the invisible states or hidden sector must come from
interactions involving gauge-invariant combinations of SM fields. Of these
interactions, those which are renormalisable, i.e. dimension-four operators,
are sometimes referred to as portals. For example, the generic fermion field
$\psi$ may correspond to a right-handed gauge singlet fermion $N_{R}$, with
the well-known mixing portal $-Y_{\nu}\bar{l}\tilde{H}N_{R}$ and also a
Majorana mass term $-\frac{1}{2}M_{R}\bar{N}_{R}^{c}N_{R}$ if lepton number is
not conserved. These terms result in the Type-I seesaw mechanism, generating
masses for the light neutrinos. A generic vector field $V_{\mu}$ may instead
exhibit kinetic mixing with the SM hypercharge, e.g.
$-\frac{\varepsilon}{4}B_{\mu\nu}V^{\mu\nu}$, with
$V_{\mu\nu}=\partial_{\mu}V_{\nu}-\partial_{\nu}V_{\mu}$. Recently Dror  _et
al._ (2017); Di Luzio  _et al._ (2022) it has been shown new massive vector
fields may also couple directly to SM currents as $\mathcal{J}_{\mu}V^{\mu}$,
even if $\mathcal{J}$ is not conserved. Finally, the so-called Higgs portal,
with the terms $\mu^{\prime}(H^{\dagger}H)\phi$ and
$\lambda^{\prime}(H^{\dagger}H)\phi^{\dagger}\phi$, is another possible window
to NP in the form of a scalar $\phi$.

However, heavy NP may also mediate the interactions between new light-
invisible particles and the SM. At energies below the associated NP scale
$\Lambda$, these then manifest as higher-dimensional effective operators, i.e.

$\displaystyle\mathcal{L}=\mathcal{L}_{\text{SM}+X}+\sum_{i}C_{i}^{(d)}\mathcal{O}_{i}^{(d)}\,.$
(2)

Here, $\mathcal{L}_{\text{SM}+X}$ is the dimension-four SM Lagrangian extended
with the gauge singlet field(s) $X$, while the sum denotes a tower of
dimension-$d$ SM gauge-invariant operators $\mathcal{O}_{i}^{(d)}$. In the
standard dimensional analysis, the Wilson coefficients $C_{i}^{(d)}$ are
proportional to $4-d$ powers of the heavy NP scale $\Lambda$.

At energies below the weak scale, the $X$-extended SM gauge-invariant
effective field theory is no longer appropriate. Instead, the weak effective
theory in the broken phase of the SM, i.e. operators invariant under
$SU(3)_{c}\times U(1)_{Q}$, should be used. Since in this EFT, the left- and
right-handed chiral fermion fields carry identical charges, we can consider
the parity basis for the quark fields, i.e. either $P$-even or $P$-odd quark
currents. In the context of $B$ meson decays, the expressions for the
branching ratios are simplified considerably.

Using the parity basis for $b\to s$ transition quark currents, the
interactions of the generic invisible fields
$X\in\\{\phi,\psi,V_{\mu},\Psi_{\mu}\\}$ in the weak effective theory
(including operators up to dimension $d=6$) are as follows. Starting with the
$P$-even quark currents, the coupling of the vector quark current to invisible
states is described by the effective Hamiltonian

$\displaystyle\mathcal{H}_{\text{eff}}^{V}$
$\displaystyle\supset\bar{s}\gamma_{\mu}b\bigg{[}h_{V}V^{\mu}+\frac{g_{VV}}{\Lambda^{2}}i\phi^{\dagger}\overset{\leftrightarrow\leavevmode\nobreak\
\,}{\partial^{\mu}}\phi$ $\displaystyle\hskip
36.0001pt+\frac{f_{VV}}{\Lambda^{2}}\bar{\psi}\gamma^{\mu}\psi+\frac{f_{VA}}{\Lambda^{2}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi$
$\displaystyle\hskip
36.0001pt+\frac{F_{VV}}{\Lambda^{2}}\bar{\Psi}^{\rho}\gamma^{\mu}\Psi_{\rho}+\frac{F_{VA}}{\Lambda^{2}}\bar{\Psi}^{\rho}\gamma^{\mu}\gamma_{5}\Psi_{\rho}\bigg{]}\,,$
(3)

where the corresponding hermitian conjugate terms are implicit in all our
Hamiltonians. Note that we have omitted the operator
$g_{V}(\bar{s}\gamma_{\mu}b)\partial^{\mu}\phi/\Lambda$, because it can be
rewritten using the quark field equations of motion as $g_{S}(\bar{s}b)\phi$,
with $g_{S}=im_{b}g_{V}/\Lambda$. This leads us to the general effective
Hamiltonian of the scalar quark current,

$\displaystyle\mathcal{H}_{\text{eff}}^{S}$
$\displaystyle\supset\bar{s}b\bigg{[}g_{S}\phi+\frac{g_{SS}}{\Lambda}\phi^{\dagger}\phi+\frac{h_{S}}{\Lambda}V_{\mu}^{\dagger}V^{\mu}$
$\displaystyle\hskip
27.0pt+\frac{f_{SS}}{\Lambda^{2}}\bar{\psi}\psi+\frac{f_{SP}}{\Lambda^{2}}\bar{\psi}\gamma_{5}\psi$
$\displaystyle\hskip
27.0pt+\frac{F_{SS}}{\Lambda^{2}}\bar{\Psi}^{\rho}\Psi_{\rho}+\frac{F_{SP}}{\Lambda^{2}}\bar{\Psi}^{\rho}\gamma_{5}\Psi_{\rho}\bigg{]}\,.$
(4)

Finally, the effective Hamiltonian for the tensor quark current is

$\displaystyle\mathcal{H}_{\text{eff}}^{T}$
$\displaystyle\supset\bar{s}\sigma_{\mu\nu}b\bigg{[}\frac{h_{T}}{\Lambda}V^{\mu\nu}+\frac{f_{TT}}{\Lambda^{2}}\bar{\psi}\sigma^{\mu\nu}\psi+\frac{F_{TT}}{\Lambda^{2}}\bar{\Psi}^{\rho}\sigma^{\mu\nu}\Psi_{\rho}$
$\displaystyle\hskip
40.00006pt+\frac{F_{TS}}{\Lambda^{2}}\bar{\Psi}^{[\mu}\Psi^{\nu]}+\frac{F_{TP}}{\Lambda^{2}}\bar{\Psi}^{[\mu}\gamma_{5}\Psi^{\nu]}\bigg{]}\,,$
(5)

with
$\bar{\Psi}^{[\mu}\Gamma\Psi^{\nu]}=i(\bar{\Psi}^{\mu}\Gamma\Psi^{\nu}-\bar{\Psi}^{\nu}\Gamma\Psi^{\mu})/2$.
For the $P$-odd quark currents, the effective Hamiltonians are similar to the
expressions above. In particular, for the axial vector quark current, the
couplings to invisible states can be found by replacing
$\bar{s}\gamma_{\mu}b\to\bar{s}\gamma_{\mu}\gamma_{5}b$ and $V\to A$ in Eq.
(II), e.g. $h_{V}\to h_{A}$ and $f_{VV}\to f_{AV}$. Likewise, the couplings
for the pseudoscalar and axial tensor quark currents are found with the
replacements $\bar{s}b\to\bar{s}\gamma_{5}b$ ($S\to P$) and
$\bar{s}\sigma_{\mu\nu}b\to\bar{s}\sigma_{\mu\nu}\gamma_{5}b$
($T\to\tilde{T}$) in Eqs. (II) and (II), respectively.

For completeness, in Appendix A, we give the matching relations between the
coefficients of the weak EFT operators in Eqs. (II–II), written in the parity
basis, and the coefficients of SM gauge-invariant operators, necessarily in
the chiral basis for the quark fields, as defined in Ref. Kamenik and Smith
(2012). For example, the scalar and pseudoscalar operators
$g_{S}(\bar{s}b)\phi$ and $g_{P}(\bar{s}\gamma_{5}b)\phi$ are induced by the
SM gauge-invariant operators
$C_{d\phi}^{S,L}\bar{d}_{R}H^{\dagger}q\phi/\Lambda$ and
$C_{d\phi}^{S,R}\bar{q}Hd_{R}\phi/\Lambda$, with the matching
$g_{S(P)}=v/\sqrt{2}\Lambda\times(C_{d\phi}^{S,R}\pm C_{d\phi}^{S,L})/2$.

We now comment on some interesting formal details regarding the possible
invisible final states. If the generic scalar $\phi$ or vector $V_{\mu}$ is
charged under a dark, possibly non-Abelian, gauge group, the operators
$g_{S(P)}[\bar{s}(\gamma_{5})b]\phi$ and
$h_{V(A)}[\bar{s}\gamma_{\mu}(\gamma_{5})b]V^{\mu}$ must vanish and it is
forbidden to emit a single or vector scalar boson. If the scalar field is
neutral, i.e., a real scalar field, the operator
$g_{VV}(\bar{s}\gamma_{\mu}b)i\phi^{*}\overset{\leftrightarrow\leavevmode\nobreak\
\,}{\partial^{\mu}}\phi/\Lambda^{2}$ vanishes instead.

Figure 2: Normalised differential branching fraction of the three-body decay
$B\to K\sum X$ as a function of the momentum transfer $q^{2}$ for SM neutrinos
and additional NP light-states, $\sum
X=\nu\bar{\nu},\,\phi\bar{\phi},\,\psi\bar{\psi},\,V\bar{V},\,\Psi\bar{\Psi}$.

The generic fermion field $\psi$, if massive, may either be a Dirac or
Majorana fermion. In the latter scenario ($\psi=\psi^{c}$), only viable if
$\psi$ is a singlet under a dark gauge group, the vector and tensor fermion
bilinears vanish identically,
$\bar{\psi}\gamma_{\mu}\psi=\bar{\psi}\sigma_{\mu\nu}\psi=0$. For the other
fermion bilinears, two possible contractions are possible for the Majorana
fields, so the amplitudes for Majorana fermions are a factor of two larger
than those for Dirac fermions. On the other hand, the final phase space
integration has to take two identical particles into account, reducing the
rate by a factor of two. Similar observations apply to all scenarios with two
identical particles in the final state (if $\phi=\phi^{\dagger}$,
$V=V^{\dagger}$ or $\Psi=\Psi^{c}$).

Finally, decay rates derived from the operators
$h_{V(A)}[\bar{s}\gamma_{\mu}(\gamma_{5})b]V^{\mu}$ and
$h_{S(P)}[\bar{s}(\gamma_{5})b]V_{\mu}^{\dagger}V^{\mu}/\Lambda$ diverge in
the limit where the generic vector field $V_{\mu}$ is massless, $m_{V}\to 0$.
This is an expected behaviour because the massless limit can only be
consistently defined when $V_{\mu}$ is a gauge field. It could be the case
that $V_{\mu}$ gets its mass from a dark Higgs-like mechanism, such that
$m_{V}$ is proportional to some power of the couplings $h_{V(A)}$ multiplied
by the vacuum expectation value of the dark Higgs field. Then, the massless
limit implies that the couplings $h_{V(A)}$ vanish, avoiding the divergence.
Alternatively, operators can be constructed solely from the manifestly dark
gauge-invariant field strength tensor $V_{\mu\nu}$. Then, the operators with
couplings $h_{V}$ and $h_{A}$ are generated by applying the equations of
motion for $V_{\mu}$ to such gauge-invariant operators, i.e.
$(\bar{s}\gamma_{\mu}b)\partial_{\nu}V^{\mu\nu}/\Lambda^{2}\to
m_{V}^{2}(\bar{s}\gamma_{\mu}b)V^{\mu}/\Lambda^{2}$, and consequently the
$1/m_{V}^{2}$ divergences in decay rates are tamed via the replacement
$1/m_{V}^{2}\to 1/\Lambda^{2}$. The situation is analogous for the spin $3/2$
fermions.

With these details in mind, we derive the branching ratios for the two-body
processes $B\to K^{(*)}\phi/V$ and differential branching ratios in the
momentum transfer $q^{2}$ for the three-body processes $B\to
K^{(*)}\phi\bar{\phi}/\psi\bar{\psi}/V\bar{V}/\Psi\bar{\Psi}$, given in full
in Appendix B. We do not list the differential branching ratios for $B\to
K^{*}V\bar{V}/\Psi\bar{\Psi}$, because we find that non-zero values of the
couplings $h_{S,P,T}$ and $F_{XY}$ are not favoured by the Belle II results
Adachi  _et al._ (2023). This result is foreshadowed by Fig. 2, where we
display the normalised differential branching fractions for the three-body
decays, for the SM neutrinos $\sum X=\nu\bar{\nu}$ (blue) and the light NP
states $\sum X=\phi\bar{\phi},\psi\bar{\psi},V\bar{V},\Psi\bar{\Psi}$ (orange,
green, red and purple, respectively), with the choice $m_{X}=0.6$ GeV. The
bands illustrate the theoretical uncertainties from the $B\to K$ form factors,
implemented as outlined in Appendix C. It can be seen that the differential
$B\to K\psi\bar{\psi}$ rate peaks close to the $\psi\bar{\psi}$ threshold,
while $B\to KV\bar{V}$ and $B\to K\Psi\bar{\Psi}$ peak closer to the endpoint.
The distribution of the scalar final states exhibits a relatively flat $q^{2}$
distribution. Since most of the signal observed by Belle II peaks at low
$q^{2}$, see Fig. 3, this generically disfavours the three-body decay
scenarios with spin 1 and 3/2 final states.

##  III Methodology

Next, we summarise how a possible NP contribution to the $B\to
K^{(*)}E_{\text{miss}}$ decay translates to a signal in the Belle II and BaBar
experiments. We also describe how the expected SM plus NP signal is used to
construct likelihoods for the different NP scenarios, given the distribution
of events seen in each experiment. Finally, we briefly describe the recast of
the upper limit on $B_{s}\to E_{\rm miss}$ decays from ALEPH Alonso-Álvarez
and Escudero (2023); Barate  _et al._ (2001).

Figure 3: The number of $B^{+}\to K^{+}E_{\text{miss}}$ events in
$q_{\text{rec}}^{2}$ in the $\eta(\text{BDT}_{2})>0.92$ Belle II signal
region, with the total background subtracted (black dots, with error bars
showing the total statistical uncertainty). Shown for comparison is the
predicted SM distribution, before (blue, solid) and after (blue, dashed) re-
scaling by the best-fit signal strength, and also after adding the
distribution for $B^{+}\to K^{+}\bar{\psi}\psi$ (green) for
$m_{\psi}=0.6\leavevmode\nobreak\ \text{GeV}$ and
$f_{VV}/\Lambda^{2}=1.7\times 10^{-2}\leavevmode\nobreak\ \text{TeV}^{-2}$ and
for $B^{+}\to K^{+}V$ (red) for $m_{V}=2.1\leavevmode\nobreak\ \text{GeV}$ and
$h_{V}=7.1\times 10^{-9}$.

The Belle II experiment conducted a search for the SM process $B^{+}\to
K^{+}\nu\bar{\nu}$ using $e^{+}e^{-}\to B^{+}B^{-}$ at the $\Upsilon(4S)$
resonance with an integrated luminosity of
$\mathcal{L}=362\leavevmode\nobreak\ \text{fb}^{-1}$ Adachi  _et al._ (2023).
The collaboration used two methods: an inclusive (ITA) and hadronic (HTA)
tagging analysis. The former exploits inclusive properties of the $B^{+}B^{-}$
pair, while the latter uses an explicit reconstruction of the partner $B$
meson via its hadronic decay. Thus, the former method trades a higher signal
efficiency for an increased background. The known backgrounds, separated as
$B^{+}B^{-}$, $B^{0}\bar{B}^{0}$ and continuum components, were suppressed by
utilising the kinematic properties of the decay in a multivariate classifier,
i.e. training a boosted decision tree (BDT). Cuts are then placed on the
classifier used for final event selection, $\text{BDT}_{2}$, to optimise the
signal region. Using the lower purity signal region,
$\eta(\text{BDT}_{2})>0.92$, where $\eta(\text{BDT}_{2})$ is defined in Eq.
(4) of Ref. Adachi  _et al._ (2023), Belle II observed a signal corresponding
to a branching ratio $2.9\sigma$ or a factor of four above the SM prediction.
The observed events minus the backgrounds are shown in bins of the
reconstructed momentum transfer, $q_{\text{rec}}^{2}$, in Fig. 3.

The BaBar experiment has also performed a search for the $B\to
K^{(*)}\nu\bar{\nu}$ process Lees  _et al._ (2013). Using $e^{+}e^{-}$
collisions, again at the $\Upsilon(4S)$ resonance and with an integrated
luminosity of $\mathcal{L}=429\leavevmode\nobreak\ \text{fb}^{-1}$, BaBar used
the HTA method to search for both neutral and charged invisible $B$ decays.
Using cuts to define a signal region minimises the $B^{+}B^{-}$ and continuum
backgrounds, separated as those which are correctly reconstructed (peak) and
those which are not (combinatorial). Their results were consistent with the
SM, and no significant signal was observed, resulting in upper bounds between
six and twelve times the SM rates.

For the SM or invisible final state(s), $\sum X$, we determine the
distribution of Belle II and BaBar events in the reconstructed momentum
transfer, $q_{\text{rec}}^{2}$, as

$\displaystyle\frac{dN_{\text{SM}(X)}}{dq^{2}_{\mathrm{rec}}}=N_{B}\int
dq^{2}f_{q^{2}_{\mathrm{rec}}}(q^{2})\epsilon(q^{2})\frac{d\mathcal{B}_{\text{SM}(X)}}{dq^{2}}\,,$
(6)

where $N_{B}$ is the number of $B^{+}B^{-}$ or $B^{0}\bar{B}^{0}$ pairs,
$f_{q^{2}_{\mathrm{rec}}}(q^{2})$ corresponds to the smearing of
$q_{\text{rec}}^{2}$ with respect to the true momentum transfer $q^{2}$, and
$\epsilon(q^{2})$ is the detector efficiency as a function of momentum
transfer. The kinematic range is $(\sum
m_{X})^{2}<q^{2}<(m_{B}-m_{K^{(*)}})^{2}$. We provide details on how we
obtained $f_{q^{2}_{\mathrm{rec}}}(q^{2})$ and $\epsilon(q^{2})$ for the Belle
II and BaBar analyses in Appendix D.

Figure 4: Minimised log-likelihood ratio of the SM plus $\phi/V$ scenario with
respect to the SM-only hypothesis as a function of the scalar/vector mass
$m_{\phi/V}$. The log-likelihood ratio $t_{\phi/V}$ for single effective
coupling contributions $g_{S}$, $h_{V}$ and $h_{T}/\Lambda$ to $B\to K\phi/V$
is minimised over the coupling and nuisance parameters
$(\bm{\theta}_{x},\tau_{b})$ for each $m_{\phi/V}$. As the likelihood is
independent of the nature of the light NP state and the coupling in two-body
decays, the three lines overlap. The blue dashed line corresponds to the
minimum of the log-likelihood ratio of the re-scaled SM hypothesis.

We obtain the expected SM ($X$) signal $s_{\text{SM}(X)}^{i}$ for the bin $i$,
$[q^{2}_{\mathrm{rec},i},q^{2}_{\mathrm{rec},i+1}]$, as

$\displaystyle
s^{i}_{\text{SM}(X)}=\int_{q^{2}_{\mathrm{rec},i}}^{q^{2}_{\mathrm{rec},i+1}}dq^{2}_{\mathrm{rec}}\,\frac{dN_{\text{SM}(X)}}{dq^{2}_{\mathrm{rec}}}\,,$
(7)

and the total expected event count in bin $i$ is,

$\displaystyle n_{\text{exp}}^{i}$
$\displaystyle=\mu\left(1+\theta_{\text{SM}}^{i}\right)s_{\text{SM}}^{i}+\left(1+\theta_{X}^{i}\right)s_{X}^{i}(m_{X},c_{X})$
$\displaystyle\quad+\sum_{b}\tau_{b}(1+\theta_{b}^{i})b^{i}\,,$ (8)

where $\mu$ is a signal strength parameter allowing to re-scale the SM signal
$s_{\text{SM}}^{i}$, the NP signal $s_{X}^{i}$ depends on the invisible
particle mass $m_{X}$ and coupling $c_{X}$, $b^{i}$ is the binned expected
signal for the background $b$, and $\tau_{b}$ is an overall normalisation. We
account for systematic and Monte-Carlo statistical uncertainties via the
nuisance parameters $\theta_{x}$.

Figure 5: Minimised log-likelihood ratios of the SM plus three-body $\sum X$
scenarios with respect to the SM-only hypothesis as a function of the new
particle $X$ mass $m_{X}$ for $\sum X=\psi\bar{\psi}$ (top) and $\sum
X=\phi\bar{\phi}/V\bar{V}$ (bottom). The log-likelihood ratio $t_{X}$ for
single effective coupling contributions $c_{X}$ to $B\to K\sum X$ is minimised
over the coupling and nuisance parameters $(\bm{\theta}_{x},\tau_{b})$ for
each $m_{X}$. The blue dotted line again corresponds to the re-scaled SM
hypothesis.

The combined likelihood then takes the following form

$\displaystyle L_{\text{SM}+X}$
$\displaystyle=\prod_{i}^{N_{\text{bins}}}\text{Poiss}\left[n_{\text{obs}}^{i},n_{\text{exp}}^{i}(\mu,m_{X},c_{X},\bm{\theta}_{x},\tau_{b})\right]$
$\displaystyle\quad\times\prod_{x={\rm
SM},X,b}\mathcal{N}\left(\bm{\theta}_{x};\bm{0},\Sigma_{x}\right)\prod_{b}\mathcal{N}\left(\tau_{b};0,\sigma_{b}^{2}\right)\,,$
(9)

where the event counts in each bin $n_{i}$ are Poissonian (Poiss) distributed,
while we sample the nuisance parameters $\bm{\theta}_{x}$ from a multi-normal
($\cal N$) distribution centred at $\bm{0}$ and of covariance $\Sigma_{x}$.
Finally, the overall normalisations of the backgrounds are sampled from a
univariate normal distribution, centred at $0$ and with standard deviation
$\sigma_{b}$.

The covariance for the SM signal, $\Sigma_{\text{SM}}$, is found by performing
a Monte-Carlo simulation of the SM signal, $s_{\text{SM}}^{i}$, including the
uncertainties on the efficiency and $B\to K^{(*)}$ form factors in Eq. (7).
The $q^{2}$ smearing introduces correlations among the $q^{2}_{\text{rec}}$
bins. The covariances for the background components, separated in the Belle II
and BaBar analyses as discussed above Eq. (7), are found by simply re-scaling
$\Sigma_{\text{SM}}$ to the relative size of the background. Finally, for the
covariance for the NP signal, we take
$(\Sigma_{X})_{ij}=s_{X}^{i}\delta_{ij}$, i.e. Poissonian uncertainties and
neglecting correlations between bins. This speeds up considerably the
following analysis, and we find that including the correlations has a
negligible impact on the results.

In the following, we consider three types of signal hypothesis, corresponding
to: (i) the SM-only scenario, where $\mu=1$ and $s_{X}^{i}=0$, (ii) the re-
scaled SM scenario, where $\mu$ is treated as free nuisance parameter and
$s_{X}^{i}=0$, and lastly (iii) the SM plus various NP scenarios, with $\mu=1$
and $s_{X}^{i}\neq 0$, considering separately each NP final state $\sum X$ and
its possible couplings $c_{X}$. On the one hand, the first two hypotheses
serve to cross-check the validity of our likelihood recast, see Appendix D for
details. On the other hand, they represent important likelihood benchmarks
against which we compare all of the NP scenarios.

Figure 6: Favoured regions in the parameter space $(g_{S},g_{P})$ relevant to
the two-body decay $B\to K^{(*)}\phi$, with the mass of $\phi$ fixed to a
value close to the best-fit point, $m_{\phi}=2.1$ GeV. Dashed and dotted grey
lines correspond to 1-$\sigma$ confidence regions for the Belle II and BaBar
analyses, respectively. Solid dark, medium and light orange lines correspond
to the 1-, 2- and 3-$\sigma$ combined confidence regions. Diagonal grey dashed
lines denote the chiral combinations of couplings $g_{S}=\pm g_{P}$.

For the SM + NP hypothesis, we define for convenience the log-likelihood ratio

$\displaystyle t_{X}=-2\ln\frac{L_{\text{SM}+X}}{L_{\text{SM}}}\,,$ (10)

where $L_{\text{SM}}$ is the likelihood for the SM-only hypothesis. In the
next section, we will use $t_{X}|_{\text{min}}$ (in which $L_{\text{SM}+X}$
and $L_{\text{SM}}$ are minimised with respect to the nuisance parameters
$(\bm{\theta}_{x},\tau_{b})$ and the NP couplings or $\mu$, respectively) to
see which invisible final states provide a better fit to the Belle II excess
compared to the (re-scaled) SM. We then examine the profile likelihoods
$\hat{t}_{X}=t_{X}-t_{X}|_{\text{min}}$, where $t_{X}$ is not minimised with
respect to the NP couplings, to infer what NP couplings are implied by the
data.

Lastly, we take into account the recent recast Alonso-Álvarez and Escudero
(2023), of the ALEPH search for $b\to\tau^{-}\bar{\nu}_{\tau}X$ at LEP Barate
_et al._ (2001) in terms of an upper limit on invisible $B_{s}$ decays,
$\mathcal{B}(B_{s}\to E_{\rm miss})<5.4\times 10^{-4}$ (90% CL). For $\sum
X=\phi\bar{\phi}/\psi\bar{\psi}/V\bar{V}/\Psi\bar{\Psi}$ the decays
$B_{s}\to\sum X$ depend on some of the same couplings as the $B\to K^{(*)}\sum
X$ processes. We thus include this constraint in the likelihood for these
scenarios as follows. Ref. Alonso-Álvarez and Escudero (2023) provides the
upper limit on the number of signal events at two confidence levels; from this
information, the mean and standard deviation of the signal and thus branching
fraction can be determined, which in turn can be used to construct a simple
Gaussian log-likelihood with the NP prediction for the branching ratio, given
in Appendix B.

##  IV Results

Figure 7: Favoured regions in the parameter space $(g_{PS},g_{SS})$ (left) and
$(g_{VV},g_{AV})$ (right) relevant to the three-body decay $B\to
K^{(*)}\phi\bar{\phi}$, with the mass of $\phi$ fixed to a value close to the
respective best-fit point, $m_{\phi}=0.52$ GeV and $m_{\phi}=0.38$ GeV, and
the heavy NP scale fixed at $\Lambda=10$ TeV. Colour coding is the same as
Fig. 6. The red dashed lines indicate the 1-$\sigma$ confidence regions from
ALEPH.

In the following, we present the results of minimising the binned log-
likelihood ratio $t_{X}$, defined in Eq. (10), with respect to the model and
nuisance parameters. To gain a better understanding of which NP scenarios are
favoured by the excess seen by Belle II, in each case, we find
$t_{X}|_{\text{min}}$ for different values of the light invisible particle
masses, $m_{X}$, using only the Belle II data. For the scenarios that provide
a better fit to the Belle II data compared to the re-scaled SM prediction, we
next profile over $t_{X}$ in the parameter space of effective couplings
$c_{X}$, also including the BaBar and ALEPH data. In this way, we identify
ranges of couplings compatible with both the Belle II excess as well as the
BaBar and ALEPH upper limits on $\mathcal{B}(B\to K^{*}E_{\text{miss}})$ and
$\mathcal{B}(B_{s}\to E_{\rm miss})$, respectively. In the following, all
minimisation is performed using the iminuit interface to the minuit2 package
Dembinski and et al. (2020); James and Roos (1975).

In Fig. 4, we first show the outcome of minimising the binned likelihood ratio
$t_{\phi/V}$ for the two-body decay scenarios $\sum X=\phi/V$. For different
values of the masses $m_{\phi/V}$, the minimisation is performed with respect
to the nuisance parameters ($\bm{\theta}_{x},\tau_{b})$ and the NP couplings;
$g_{S}$ for the scalar boson and $h_{V}$ or $h_{T}/\Lambda$ for the vector
boson. Firstly, we see that the minimised likelihoods in the three scenarios
overlap because the likelihood is independent of the spin and coupling of
single-particle $\sum X$; the likelihood is constructed in bins of
$q_{\text{rec}}^{2}$, while the two-body scenario only depends on
$q_{\text{rec}}^{2}$ through the smearing at the true momentum transfer
$q^{2}=m_{\phi/V}^{2}$. For masses in the range $1.7\leavevmode\nobreak\
\text{GeV}\lesssim m_{\phi/V}\lesssim 2.4\leavevmode\nobreak\ \text{GeV}$, the
SM plus an invisible scalar or vector provides a better fit than the re-scaled
SM prediction, shown as a blue dotted line. The value $m_{\phi/V}=(2.1\pm
0.1)$ GeV provides the best fit to the data, with a significance of
$4.5\sigma$ over the SM.

Figure 8: Two-dimensional favoured regions in the parameter space
$(f_{VV},f_{VA},f_{AV},f_{AA})$ relevant to the three-body decay $B\to
K^{(*)}\psi\bar{\psi}$, with the mass of $\psi$ fixed to a value close to the
best-fit point, $m_{\psi}=0.6$ GeV and the heavy NP scale fixed at
$\Lambda=10$ TeV. Colour coding is the same as Fig. 6, with green lines
denoting the 1-, 2- and 3-$\sigma$ combined confidence regions instead.
Additionally, the cyan dashed region corresponds to the 1-$\sigma$ combined
confidence region in the case that $\psi$ is a Majorana fermion.

In Fig. 5 (above), we show $t_{\psi}|_{\text{min}}$ for the scenario of the SM
plus an invisible fermion pair, $\sum X=\psi\bar{\psi}$. Again, for different
values of the fermion mass $m_{\psi}$, $t_{\psi}$ is minimised with respect to
$(\bm{\theta}_{x},\tau_{b})$ and one $f_{XY}$ coupling, with the others set to
zero. For masses $m_{\psi}\lesssim 0.85$ GeV, the vector couplings $f_{VV}$
and $f_{VA}$ provide an improved fit compared to the re-scaled SM, with
$m_{\psi}=0.60_{-0.14}^{+0.11}$ GeV giving a highest significance of
$3.7\sigma$. While still yielding an improvement over the SM-only hypothesis,
the scalar and tensor couplings are not as competitive. In Fig. 5 (below), we
show $t_{X}|_{\text{min}}$ for the other three-body decay scenarios; the SM
plus pairs of scalars, vectors or spin 3/2 fermions, $\sum
X=\phi\bar{\phi},V\bar{V},\Psi\bar{\Psi}$. For the scalar pair, both of the
couplings $g_{SS}$ (orange, solid) and $g_{VV}$ (orange, dashed) provide a
better fit than the re-scaled SM for masses in the ranges
$0.05\leavevmode\nobreak\ \text{GeV}\lesssim m_{\phi}\lesssim
0.76\leavevmode\nobreak\ \text{GeV}$ and $m_{\phi}\lesssim
0.63\leavevmode\nobreak\ \text{GeV}$, respectively. The scalar masses
$m_{\phi}=0.52^{+0.11}_{-0.14}$ GeV and $m_{\phi}=0.38^{+0.13}_{-0.15}$ GeV
provide the best fits at $3.4\sigma$. It is clear from Fig. 5 that the three-
body vector (red, solid) and spin $3/2$ (purple, solid) scenarios do not
provide a good fit to the Belle II excess, as expected from the $q^{2}$
distributions in Fig. 2.

From Figs. 4 and 5, we see that only the two-body $B\to K\phi/V$ and three-
body $B\to K\phi\bar{\phi}/\psi\bar{\psi}$ scenarios provide good fits to the
Belle II excess. For these cases, we now investigate which couplings are
favoured. To do so, we turn on two couplings at a time and calculate the
profile log-likelihood, $\hat{t}_{X}=t_{X}-t_{X}|_{\text{min}}$, where $t_{X}$
is calculated with the couplings at a particular point in the parameter space
(minimising over the nuisance parameters) and $t_{X}|_{\text{min}}$ is the
likelihood at the global minimum, i.e., minimised with respect to all
parameters. However, in both cases, we keep the SM contribution fixed,
$\mu=1$.

Starting with the $B\to K\phi$ mode, we show in Fig. 6 the constraints in the
$(g_{S},g_{P})$ parameter space, with the scalar mass fixed to $m_{\phi}=2.1$
GeV. The separate Belle II and BaBar 1-$\sigma$ regions are shown as grey
dotted and dashed lines, respectively. The 1-, 2- and 3-$\sigma$ combined
confidence regions are depicted as dark, medium and light orange-shaded
contours. Belle II favours the coupling $|g_{S}|=(1.6\pm 0.2)\times 10^{-8}$,
while BaBar imposes the upper limit $|g_{P}|<2.5\times 10^{-8}$ (90% CL). In
Fig. 7, we instead show the favoured region in the $(g_{SS},g_{PS})$ and
$(g_{VV},g_{AV})$ parameter spaces, relevant for $B\to K^{(*)}\phi\bar{\phi}$
and $B_{s}\to\phi\bar{\phi}$ processes. For $\Lambda=10$ TeV and the scalar
mass fixed to $m_{\phi}=0.52$ GeV and $0.38$ GeV for the scalar and vector
couplings, respectively, the Belle II data favours the coupling values
$|g_{SS}|=(9.3\pm 1.4)\times 10^{-4}$ and $|g_{VV}|=(3.6\pm 0.5)$, while BaBar
and ALEPH put upper bounds on the couplings $g_{PS}$ and $g_{AV}$; ALEPH
places a stronger bound on the former, $|g_{PS}|<2.1\times 10^{-3}$, and BaBar
on the latter, $|g_{AV}|<7.2$, both at 90% CL.

Figure 9: Favoured regions in the parameter space $(h_{V},h_{A})$ (left) and
$(h_{T},h_{\tilde{T}})$ (right) relevant to the two-body decay $B\to
K^{(*)}V$, with the mass of $V$ fixed to a value close to the best-fit point,
$m_{V}=2.1$ GeV and the heavy NP scale fixed at $\Lambda=10$ TeV. Colour
coding is the same as Fig. 6.

In Fig. 8, we show the constraints on the parameter space
$(f_{VV},f_{VA},f_{AV},f_{AA})$ which enter the decay rates for $B\to
K^{(*)}\psi\bar{\psi}$ and $B_{s}\to\psi\bar{\psi}$. We show planes for each
of the six combinations of the couplings, with the others set to zero, for the
choices $m_{\psi}=0.6$ GeV and $\Lambda=10$ TeV. Again, the dashed and dotted
grey lines show the 1-$\sigma$ regions from Belle II and BaBar on the
invisible Dirac fermion scenario, with the dark, medium and light green
regions the 1-, 2- and 3-$\sigma$ combined confidence regions. As seen in the
lower left plot of Fig. 8, the Belle II data suggests the values
$|f_{VV}|^{2}+|f_{VA}|^{2}=(3.0\pm 0.8)$ for $f_{AV}=f_{AA}=0$, while the
centre plot shows the BaBar upper bounds $|f_{AV}|<2.8$ and $|f_{AA}|<3.0$ for
$f_{VV}=f_{VA}=0$ (90% CL). The ALEPH bound on
$\mathcal{B}(B_{s}\to\psi\bar{\psi})$ translates to the upper limit
$\big{|}f_{AA}\big{|}<13.8$ (90% CL), and so BaBar provides more stringent
constraints. In Fig. 8, we also show the 1-$\sigma$ combined region for the
invisible Majorana fermion scenario (cyan shaded). These constraints do not
depend on the couplings $f_{VV}$ and $f_{AV}$, which vanish.

Finally, in Fig. 9, we show the favoured regions in the $(h_{V},h_{A})$ and
$(h_{T},h_{\tilde{T}})$ parameter spaces, which are probed by $B\to K^{(*)}V$,
for $m_{V}=2.1$ GeV and $\Lambda=10$ TeV. The colour coding for the Belle II
and BaBar 1-$\sigma$ regions is the same as before, and the dark, medium and
light red shaded contours again show the 1-, 2- and 3-$\sigma$ combined
confidence regions. Belle II favours the couplings values $|h_{V}|=(7.1\pm
0.7)\times 10^{-9}$ and $|h_{T}|=(4.7\pm 0.5)\times 10^{-5}$, while BaBar
enforces the upper bounds $h_{A}<7.6\times 10^{-9}$ (for $h_{V}$ at the best-
fit point) and $|h_{T}|^{2}+|h_{\tilde{T}}|^{2}<2.8\times 10^{-10}$ (90% CL).
In this case the BaBar data disfavour the $h_{T}$ value preferred by Belle II.

Several studies have previously considered the Belle II excess in terms of new
invisible particles Berezhnoy and Melikhov (2023); Felkl  _et al._ (2023);
Chen and Chiang (2023); Dreiner  _et al._ (2023); Altmannshofer  _et al._
(2023); He  _et al._ (2024). However, most derive their results solely from
$B\to K^{(*)}E_{\rm miss}$ integrated branching fractions and only one study
considered also the decay $B_{s}\to E_{\rm miss}$ Felkl  _et al._ (2023). Only
in Ref. Fridell  _et al._ (2023) was a binned likelihood analysis of $B\to
KE_{\rm miss}$ performed within different new physics scenarios relying on
both Belle II and BaBar data. We observe that the preferred mass $m_{\phi/V}$
and $m_{\psi}$ obtained from the two- and three-body decay fits respectively
are compatible with Refs. Altmannshofer  _et al._ (2023); Fridell  _et al._
(2023). Additionally, our constraint on the $h_{V}$ and $h_{A}$ couplings is
compatible with the results from Ref. Altmannshofer  _et al._ (2023)111The
normalisation of Ref. Altmannshofer  _et al._ (2023) differs by a factor of 2,
obtaining $g_{V(A)}^{(4)}=2h_{V(A)}$.. On the other hand, our estimated
significances of NP scenarios over the SM differ somewhat from those given in
Ref Fridell  _et al._ (2023). In particular, we find that two-body decay
kinematics seems to give a better fit to Belle II data than three-body decay
spectra.

##  V Conclusions

In the present work, we have explored the interesting possibility that the
signal of $B\to KE_{\rm miss}$ recently observed by Belle II arises not only
from the SM neutrinos but from new undetected particles in the final state. By
taking into account available experimental data on the total event yields, as
well as on the differential $B\to K^{(*)}E_{\rm miss}$ distributions presented
in Belle II and BaBar analyses, we were able to construct the likelihood for
different NP scenarios using experimental information from both processes. In
addition, we considered constraints on the branching ratio $B_{s}\to E_{\rm
miss}$ obtained from the ALEPH data.

We considered possible decay channels involving one new invisible scalar or
vector state as well as decays to two invisible states, either scalars,
spin-$1/2$ fermions, vectors or even spin-$3/2$ states.

The two-body $B\to K^{(*)}X$ topology with an invisible particle mass of
$m_{X}=2.1$ GeV can accommodate all data much better than the (re-scaled) SM.
In particular, the Belle II data prefer these scenarios by ($3.5\sigma$)
$4.5\sigma$, respectively. This conclusion seems to differ somewhat from the
results of Ref. Fridell  _et al._ (2023).

In the case of two invisible new particles in the final state, i.e. $B\to
K^{(*)}X\bar{X}$, we find that two scalars or two spin-$1/2$ fermions can also
accommodate all data with the most favoured masses in the range
$m_{X}=0.4-0.6$ GeV, depending on the $X$ spin and couplings considered.
Conversely, scenarios with pairs of vectors or spin $3/2$ fermions in the
final state cannot improve the SM fit of Belle II and are thus disfavoured
irrespective of their mass or couplings.

In the parity basis of effective operators, we find that scenarios with
spin-$1/2$ fermions coupled to the vector quark current are preferred over
those coupling to scalar or tensor quark currents. On the other hand,
scenarios with pairs of scalar particles coupled to vector or scalar quark
currents are equally likely.

The results of the BaBar search for $B\to K^{*}E_{\rm miss}$ are consistent
with most scenarios preferred by Belle II data. In the parameter space of
possible couplings, they yield predominantly orthogonal constraints. One
notable exception is a single massive vector field coupled to the tensor
(dipole) quark current. The Belle II and BaBar results cannot be accommodated
simultaneously in this scenario. Finally, the ALEPH data on $B_{s}\to E_{\rm
miss}$ currently provide the strongest constraint on a pair of scalars coupled
to a pseudoscalar quark bilinear operator ($\bar{s}\gamma_{5}b$).

The impact of $B\to K^{*}E_{\rm miss}$ (and $B_{s}\to E_{\rm miss}$)
constraints is somewhat more interesting in the chiral basis of quark
operators (see also Appendix A), which is often preferred for UV model
building respecting the SM chiral gauge symmetry. Interestingly, we find that
purely chiral couplings still lie within the combined $68\%$ CL region of both
Belle II and BaBar (ALEPH) results. However, a future improvement on $B\to
K^{*}E_{\rm miss}$ by Belle II could soon put most chiral scenarios under
pressure.

We also note that coupling ranges implied by our fit are mostly consistent
with NP scales even above our benchmark $\Lambda=10$ TeV. Consequently, new
degrees of freedom associated with higher dimensional operator terms in Eqs.
(II), (II), and (II) could lie above the LHC reach. Notable exceptions are
pairs of invisible scalars or fermions coupled to the quark vector current
($\bar{s}\gamma_{\mu}b$), where the fit to Belle II signal points to NP below
few TeV.

While we only considered $b\to s$ transitions in this work, in more complete
NP scenarios, other phenomenology could become relevant. In particular, NP
flavour constructions (see e.g. Badin and Petrov (2010); Bauer  _et al._
(2022); Li and Tandean (2023); He  _et al._ (2024)) would relate these decays
to $K\to\pi E_{\rm miss}$, $D\to\pi E_{\rm miss}$ and monojet/monotop
production at high $p_{T}$, leading to possibly strong constraints from NA62
Cortina Gil  _et al._ (2021), KOTO Ahn  _et al._ (2021), BESSIII Ablikim  _et
al._ (2022) and LHC Khachatryan  _et al._ (2015); Aad  _et al._ (2024)
experiments. The results of our Belle II likelihood fit however indicate that
for the preferred masses of invisible states, they could be kinematically
forbidden from being produced in rare kaon and possibly even $D$ meson decays.
Finally, in UV complete models, constraints from $B_{s}$ meson
oscillations222In scenarios with a single scalar or vector, these arise at the
tree level, but have been found in Ref. Bauer  _et al._ (2021) to be
negligible, a result which we have also verified explicitly. as well as EW
precision data could be relevant. We leave these considerations for future
work.

###### Acknowledgements.

We thank Luka Šantelj and Peter Križan for multiple illuminating discussions
regarding the details of the Belle II experimental analysis and in particular
Luka Šantelj for providing the relevant smearing and efficiency maps of the
Belle II ITA tagging analysis. We also thank Damir Bečirević and Olcyr
Sumensari for providing useful information on the relevant $B\to K^{(*)}$ form
factors. PDB, SF and JFK acknowledge financial support from the Slovenian
Research Agency (research core funding No. P1-0035, J1-3013 and N1-0321). MN
acknowledges the financial support by the Spanish Government (Agencia Estatal
de Investigación MCIN/AEI/10.13039/501100011033) and the European Union
NextGenerationEU/PRTR through the “Juan de la Cierva” program (Grant No.
JDC2022-048787-I) and Grant No. PID2020-114473GB-I00. MN also acknowledges the
support of the Generalitat Valenciana through Grant No. PROMETEO/2021/071.
This study has been partially carried out within the INFN project (Iniziativa
Specifica) QFT-HEP.

##  Appendix A Chiral Operator Basis

To convert from the parity basis (Eqs. (II), (II), and (II)) to the chiral
basis (defined in Ref. Kamenik and Smith (2012)) for the operators involving
the fields $X\in\\{\phi,\psi,V_{\mu},\Psi_{\mu}\\}$, it is convenient to
define the following matrices,

$\displaystyle P=\frac{1}{2}\begin{pmatrix}1&1\\\
-1&1\end{pmatrix}\,,\leavevmode\nobreak\ \leavevmode\nobreak\
Q=\frac{1}{4}\begin{pmatrix}1&1&1&1\\\ -1&1&-1&1\\\ -1&-1&1&1\\\ 1&-1&-1&1\\\
\end{pmatrix}\,.$ (11)

For the operators containing the dark scalar field $\phi$ in Eqs. (II) to
(II), we then have the matching relations,

$\displaystyle\begin{pmatrix}g_{S(S)}\\\ g_{P(S)}\\\
\end{pmatrix}=P\frac{v}{\sqrt{2}\Lambda}\begin{pmatrix}C_{d\phi(\phi)}^{S,L}\\\
C_{d\phi(\phi)}^{S,R}\\\ \end{pmatrix}\,,$
$\displaystyle\begin{pmatrix}g_{VV}\\\ g_{AV}\\\
\end{pmatrix}=P\begin{pmatrix}C_{d\phi\phi}^{V,L}\\\ C_{d\phi\phi}^{V,R}\\\
\end{pmatrix}\,.$ (12)

For the four-fermion operators involving the invisible fermion field $\psi$,
we have

$\displaystyle\begin{pmatrix}f_{VV}\\\ f_{VA}\\\ f_{AV}\\\ f_{AA}\\\
\end{pmatrix}=Q\begin{pmatrix}C_{d\psi}^{V,LL}\\\ C_{d\psi}^{V,LR}\\\
C_{d\psi}^{V,RL}\\\ C_{d\psi}^{V,RR}\\\ \end{pmatrix}\,,\leavevmode\nobreak\
\begin{pmatrix}f_{SS}\\\ f_{SP}\\\ f_{PS}\\\ f_{PP}\\\
\end{pmatrix}=Q\frac{v}{\sqrt{2}\Lambda}\begin{pmatrix}C_{d\psi}^{S,LL}\\\
C_{d\psi}^{S,LR}\\\ C_{d\psi}^{S,RL}\\\ C_{d\psi}^{S,RR}\\\ \end{pmatrix}\,,$
$\displaystyle\begin{pmatrix}f_{TT}\\\ f_{\tilde{T}T}\\\
\end{pmatrix}=P\frac{v}{\sqrt{2}\Lambda}\begin{pmatrix}C_{d\psi}^{T,LL}\\\
C_{d\psi}^{T,RR}\\\ \end{pmatrix}\,.$ (13)

The same relations hold for the effective couplings of the spin 3/2 field
$\Psi_{\mu}$, with the replacements $f_{XY}\to F_{XY}$ and $\psi\to\Psi$.
Additionally, the couplings $F_{TS}$, $F_{TP}$, $F_{\tilde{T}S}$ and
$F_{\tilde{T}P}$ are related to chiral basis coefficients in a similar manner
as the scalar and pseudoscalar couplings. Finally, the operators coupling the
vector field $V_{\mu}$ to the relevant quark couplings are rotated to the
chiral basis as

$\displaystyle\begin{pmatrix}h_{V}\\\ h_{A}\\\
\end{pmatrix}=P\begin{pmatrix}C_{dV}^{V,L}\\\ C_{dV}^{V,R}\\\
\end{pmatrix}\,,\quad\begin{pmatrix}h_{T}\\\ h_{\tilde{T}}\\\
\end{pmatrix}=P\frac{v}{\sqrt{2}\Lambda}\begin{pmatrix}C_{dV}^{T,L}\\\
C_{dV}^{T,R}\\\ \end{pmatrix}\,,$ $\displaystyle\begin{pmatrix}h_{S}\\\
h_{P}\\\
\end{pmatrix}=P\frac{v}{\sqrt{2}\Lambda}\begin{pmatrix}C_{dVV}^{S,L}\\\
C_{dVV}^{S,R}\\\ \end{pmatrix}\,.$ (14)

In Figs. 6, 7 and 9, the constraints on the chiral-basis couplings can be
inferred from the light grey dashed lines, which show, for example,
$g_{S}=g_{P}$ and $g_{S}=-g_{P}$ (which equivalently depict
$C_{d\phi}^{S,L}=0$ and $C_{d\phi}^{S,R}=0$, respectively).

The transformation to the chiral basis for the constraints on the fermion
couplings in Fig. 8 is slightly less trivial than a rotation because two of
couplings $(f_{VV},f_{VA},f_{AV},f_{AA})$ are taken to be zero in each of the
six plots. We instead perform an additional fit for the chiral couplings
$(C_{d\psi}^{LL},C_{d\psi}^{LR},C_{d\psi}^{RL},C_{d\psi}^{RR})$, with the
favoured regions shown in Fig. 10. The dashed and dotted grey lines again show
the 1-$\sigma$ regions from Belle II and BaBar, while the dark, medium and
light green contours show the 1-, 2- and 3-$\sigma$ combined confidence
regions.

Figure 10: Two-dimensional favoured regions in the parameter space
$(C_{d\psi}^{LL},C_{d\psi}^{LR},C_{d\psi}^{RL},C_{d\psi}^{RR})$ relevant to
the three-body decay $B\to K^{(*)}\psi\bar{\psi}$, with the mass of $\psi$
fixed to a value close to the best-fit point, $m_{\psi}=0.6$ GeV and the heavy
NP scale fixed at $\Lambda=10$ TeV. Colour coding is the same as Fig. 6, with
green lines denoting the 1-, 2- and 3-$\sigma$ combined confidence regions.

##  Appendix B Analytic Expressions of Decay Rates

In this appendix, we give explicit expressions for the relevant two-body and
three-body $B$ decays to $K^{(*)}$ and the invisible final states
$X\in\\{\phi,\psi,V,\Psi\\}$, first explored in Ref. Kamenik and Smith (2012).
We also give the decay rates for invisible two-body $B_{s}$ decays. The
results make use of the form factors and decay constants defined in Appendix
C.

###  B.1 $B\to K$ Two-Body Decays

The two-body decay rates for $B\to K^{(*)}\phi$ are

$\displaystyle\Gamma(B\to K\phi)$
$\displaystyle=\frac{|g_{S}|^{2}}{8\pi}\frac{\left|\vec{p}_{K}\right|m_{B}^{2}\delta_{K}^{2}}{(m_{b}-m_{s})^{2}}f_{0}^{2}(m_{\phi}^{2})\,,$
(15) $\displaystyle\Gamma(B\to K^{*}\phi)$
$\displaystyle=\frac{|g_{P}|^{2}}{2\pi}\frac{\left|\vec{p}_{K^{*}}\right|^{3}}{(m_{b}+m_{s})^{2}}A_{0}^{2}(m_{\phi}^{2})\,,$
(16)

where $\delta_{K^{(*)}}\equiv(1-m_{K^{(*)}}^{2}/m_{B}^{2})$, $g_{S}$ ($g_{P}$)
is the scalar (pseudoscalar) coupling of the scalar field $\phi$ to the
flavour-changing quark current $\bar{s}b$ ($\bar{s}\gamma_{5}b$), and the
three-momentum of $K^{(*)}$ is

$\displaystyle\left|\vec{p}_{K^{(*)}}\right|=\frac{\lambda^{1/2}(m_{B}^{2},m_{\phi}^{2},m_{K^{(*)}}^{2})}{2m_{B}}\,,$
(17)

with $\lambda(x,y,z)=(x-y-z)^{2}-4yz$. Note that the above rates also apply to
the operators $(\bar{s}\gamma_{\mu}b)\partial^{\mu}\phi/\Lambda$ and
$(\bar{s}\gamma_{\mu}\gamma_{5}b)\partial^{\mu}\phi/\Lambda$ if one makes the
replacement $g_{S(P)}\to im_{b}g_{V(A)}/\Lambda$.

The two-body decay rates for $B\to K^{(*)}V$ are instead given by

$\displaystyle\Gamma(B\to KV)$
$\displaystyle=\frac{|\vec{p}_{K}|}{2\pi}\Bigg{[}|h_{V}|^{2}\frac{|\vec{p}_{K}|^{2}}{m_{V}^{2}}f_{+}^{2}(m_{V}^{2})+4|h_{T}|^{2}\frac{|\vec{p}_{K}|^{2}}{\Lambda^{2}}\frac{m_{V}^{2}}{(m_{B}+m_{K})^{2}}f_{T}^{2}(m_{V}^{2})$
$\displaystyle\hskip
40.00006pt-4\,\Re[h_{V}h_{T}^{*}]\frac{|\vec{p}_{K}|^{2}}{(m_{B}+m_{K})\Lambda}f_{+}(m_{V}^{2})f_{T}(m_{V}^{2})\Bigg{]}\,,$
(18) $\displaystyle\Gamma(B\to K^{*}V)$
$\displaystyle=\frac{|\vec{p}_{K^{*}}|}{2\pi}\Bigg{[}2|h_{V}|^{2}\frac{|\vec{p}_{K^{*}}|^{2}}{(m_{B}+m_{K^{*}})^{2}}V^{2}(m_{V}^{2})+|h_{A}|^{2}\left(\frac{(m_{B}+m_{K^{*}})^{2}}{2m_{B}^{2}}A_{1}^{2}(m_{V}^{2})+\frac{16m_{K^{*}}^{2}}{m_{V}^{2}}A_{12}^{2}(m_{V}^{2})\right)$
$\displaystyle\hskip
40.00006pt+8|h_{T}|^{2}\frac{|p_{K^{*}}|^{2}}{\Lambda^{2}}T_{1}^{2}(m_{V}^{2})+2|h_{\tilde{T}}|^{2}\frac{m_{V}^{2}}{\Lambda^{2}}\bigg{(}\frac{m_{B}^{2}\delta_{K^{*}}^{2}}{m_{V}^{2}}T_{2}^{2}(m_{V}^{2})+\frac{8m_{K^{*}}^{2}}{(m_{B}+m_{K^{*}})^{2}}T_{23}^{2}(m_{V}^{2})\bigg{)}$
$\displaystyle\hskip
40.00006pt+8\,\Re[h_{V}h_{T}^{*}]\frac{|\vec{p}_{K^{*}}|^{2}}{(m_{B}+m_{K^{*}})\Lambda}V(m_{V}^{2})T_{1}(m_{V}^{2})$
$\displaystyle\hskip
40.00006pt+2\,\Re[h_{A}h_{\tilde{T}}^{*}]\frac{m_{B}+m_{K^{*}}}{\Lambda}\bigg{(}\delta_{K^{*}}A_{1}(m_{V}^{2})T_{2}(m_{V}^{2})+\frac{16m_{K^{*}}^{2}}{(m_{B}+m_{K^{*}})^{2}}A_{12}(m_{V}^{2})T_{23}(m_{V}^{2})\bigg{)}\Bigg{]}\,,$
(19)

where $h_{V}$ ($h_{A}$) is the (axial) vector coupling of the vector field
$V_{\mu}$ to the flavour-changing quark current $\bar{s}\gamma_{\mu}b$
($\bar{s}\gamma_{\mu}\gamma_{5}b$) and $h_{T}$ ($h_{\tilde{T}}$) is the
magnetic (electric) dipole-like coupling of the field strength tensor
$V^{\mu\nu}$ to the flavour-changing quark current $\bar{s}\sigma_{\mu\nu}b$
($\bar{s}\sigma_{\mu\nu}\gamma_{5}b$). The three-momentum of $K^{(*)}$ is
given by Eq. (17) with the replacement $m_{\phi}\to m_{V}$.

###  B.2 $B\to K$ Three-Body Decays

Firstly, the three-body differential decay rates in $q^{2}$ for $B\to
K^{(*)}\phi\bar{\phi}$ are given by

$\displaystyle\frac{d\Gamma(B\to K\phi\bar{\phi})}{dq^{2}}$
$\displaystyle=\frac{\beta_{\phi}}{96\pi^{3}}\frac{|\vec{p}_{K}|}{\Lambda^{2}}\Bigg{[}\frac{3}{4}|g_{SS}|^{2}\frac{m_{B}^{2}\delta_{K}^{2}}{(m_{b}-m_{s})^{2}}f_{0}^{2}(q^{2})+|g_{VV}|^{2}\frac{|\vec{p}_{K}|^{2}}{\Lambda^{2}}\beta_{\phi}^{2}f_{+}^{2}(q^{2})\Bigg{]}\,,$
(20) $\displaystyle\frac{d\Gamma(B\to K^{*}\phi\bar{\phi})}{dq^{2}}$
$\displaystyle=\frac{\beta_{\phi}}{96\pi^{3}}\frac{|\vec{p}_{K^{*}}|}{\Lambda^{2}}\Bigg{[}3|g_{PS}|^{2}\frac{|\vec{p}_{K^{*}}|^{2}}{(m_{b}+m_{s})^{2}}A_{0}^{2}(q^{2})$
$\displaystyle\hskip
75.0001pt+|g_{AV}|^{2}\frac{q^{2}}{\Lambda^{2}}\beta_{\phi}^{2}\left(\frac{(m_{B}+m_{K^{*}})^{2}}{2m_{B}^{2}}A_{1}^{2}(q^{2})+\frac{16m_{K^{*}}^{2}}{q^{2}}A_{12}^{2}(q^{2})\right)\Bigg{]}\,,$
(21)

where $\beta_{X}=\sqrt{1-4m_{X}^{2}/q^{2}}$ and the three-momentum of
$K^{(*)}$ is given by Eq. (17) with the replacement $m_{\phi}^{2}\to q^{2}$.
In the scenario of a real scalar field, where $\phi=\bar{\phi}$, the
contributions from $g_{VV}$ and $g_{AV}$ vanish, and the remaining
contributions from $g_{SS}$ and $g_{SP}$ are a factor of $2$ larger at the
amplitude level; hence, the replacements $g_{SS}\to 2g_{SS}$ and $g_{PS}\to
2g_{PS}$ can be made in Eq. (20). The rate must also be multiplied by factor
of $1/2$ to account for identical outgoing states.

The differential decay rate for $B\to K\psi\bar{\psi}$ is

$\displaystyle\frac{d\Gamma(B\to K\psi\bar{\psi})}{dq^{2}}$
$\displaystyle=\frac{\beta_{\psi}}{24\pi^{3}}\frac{|\vec{p}_{K}|q^{2}}{\Lambda^{4}}$
$\displaystyle\times\Bigg{[}\bigg{(}\Big{(}|f_{VV}|^{2}\beta_{\psi}^{\prime
2}+|f_{VA}|^{2}\beta_{\psi}^{2}\Big{)}f_{+}(q^{2})+12\,\Re\big{[}f_{VV}f_{TT}^{*}\big{]}\frac{m_{\psi}}{m_{B}+m_{K}}f_{T}(q^{2})\bigg{)}\frac{|\vec{p}_{K}|^{2}}{q^{2}}f_{+}(q^{2})$
$\displaystyle+\frac{3}{8}\bigg{(}4|f_{VA}|^{2}\frac{m_{\psi}^{2}(m_{b}-m_{s})^{2}}{q^{4}}+|f_{SS}|^{2}\beta_{\psi}^{2}+|f_{SP}|^{2}+4\,\Re\big{[}f_{VA}f_{SP}^{*}\big{]}\frac{m_{\psi}(m_{b}-m_{s})}{q^{2}}\bigg{)}\frac{m_{B}^{2}\delta_{K}^{2}}{(m_{b}-m_{s})^{2}}f_{0}^{2}(q^{2})$
$\displaystyle+2\Big{(}|f_{TT}|^{2}\beta_{\psi}^{\prime\prime
2}+|f_{\tilde{T}T}|^{2}\beta_{\psi}^{2}\Big{)}\frac{|\vec{p}_{K}|^{2}}{(m_{B}+m_{K})^{2}}f_{T}^{2}(q^{2})\Bigg{]}\,,$
(22)

where we have defined $\beta_{\psi}^{\prime
2}=(3-\beta_{\psi}^{2})/2=(1+2m_{\psi}^{2}/q^{2})$ and
$\beta_{\psi}^{\prime\prime 2}=(3-2\beta_{\psi}^{2})=(1+8m_{\psi}^{2}/q^{2})$.
Likewise, the decay rate for $B\to K^{*}\psi\bar{\psi}$ is

$\displaystyle\frac{d\Gamma(B\to K^{*}\psi\bar{\psi})}{dq^{2}}$
$\displaystyle=\frac{\beta_{\psi}}{24\pi^{3}}\frac{|\vec{p}_{K^{*}}|q^{2}}{\Lambda^{4}}$
$\displaystyle\times\Bigg{[}2\Big{(}|f_{VV}|^{2}\beta_{\psi}^{\prime
2}+|f_{VA}|^{2}\beta_{\psi}^{2}\Big{)}\frac{|\vec{p}_{K^{*}}|^{2}}{(m_{B}+m_{K^{*}})^{2}}V^{2}(q^{2})$
$\displaystyle+\Big{(}|f_{AV}|^{2}\beta_{\psi}^{\prime
2}+|f_{AA}|^{2}\beta_{\psi}^{2}\Big{)}\left(\frac{(m_{B}+m_{K^{*}})^{2}}{2m_{B}^{2}}A_{1}^{2}(q^{2})+\frac{16m_{K^{*}}^{2}}{q^{2}}A_{12}^{2}(q^{2})\right)$
$\displaystyle+\frac{3}{2}\bigg{(}4|f_{AA}|^{2}\frac{m_{\psi}^{2}(m_{b}+m_{s})^{2}}{q^{4}}+|f_{PS}|^{2}\beta_{\psi}^{2}+|f_{PP}|^{2}-4\,\Re\big{[}f_{AA}f_{PP}^{*}\big{]}\frac{m_{\psi}(m_{b}+m_{s})}{q^{2}}\bigg{)}\frac{|\vec{p}_{K^{*}}|^{2}}{(m_{b}+m_{s})^{2}}A_{0}^{2}(q^{2})$
$\displaystyle+4\bigg{(}\Big{(}|f_{TT}|^{2}\beta_{\psi}^{\prime\prime
2}+|f_{\tilde{T}T}|^{2}\beta_{\psi}^{2}\Big{)}T_{1}(q^{2})-6\,\Re\big{[}f_{VV}f_{TT}^{*}\big{]}\frac{m_{\psi}}{m_{B}+m_{K^{*}}}V(q^{2})\bigg{)}\frac{|\vec{p}_{K^{*}}|^{2}}{q^{2}}T_{1}(q^{2})$
$\displaystyle+\Big{(}|f_{TT}|^{2}\beta_{\psi}^{2}+|f_{\tilde{T}T}|^{2}\beta_{\psi}^{\prime\prime
2}\Big{)}\bigg{(}\frac{m_{B}^{2}\delta_{K^{*}}^{2}}{q^{2}}T_{2}^{2}(q^{2})+\frac{8m_{K^{*}}^{2}}{(m_{B}+m_{K^{*}})^{2}}T_{23}^{2}(q^{2})\bigg{)}$
$\displaystyle-6\,\Re\big{[}f_{AV}f_{\tilde{T}T}^{*}\big{]}\frac{m_{\psi}(m_{B}+m_{K^{*}})}{q^{2}}\left(\delta_{K^{*}}A_{1}(q^{2})T_{2}(q^{2})+\frac{16m_{K^{*}}^{2}}{(m_{B}+m_{K^{*}})^{2}}A_{12}(q^{2})T_{23}(q^{2})\right)\Bigg{]}\,.$
(23)

The differential decay rates in Eq. (B.2) and (B.2) are valid for a pair of
invisible Dirac fermions in the final state. However, the result for Majorana
fermions ($\psi=\psi^{c}$) can be obtained from Eq. (B.2) and (B.2) with the
replacements $f_{VV},f_{AV},f_{TT},f_{\tilde{T}T}\to 0$ and $f_{XY}\to
2f_{XY}$ for the remaining couplings. Additionally, a factor of $1/2$ is
required to account for identical outgoing states.

The differential decay rate for the process $B\to KV\bar{V}$ is

$\displaystyle\frac{d\Gamma(B\to KV\bar{V})}{dq^{2}}$
$\displaystyle=\frac{\beta_{V}|h_{S}|^{2}}{512\pi^{3}}\frac{|\vec{p}_{K}|q^{4}}{m_{V}^{2}\Lambda^{2}}\mathcal{J}_{V}\frac{m_{B}^{2}\delta_{K}^{2}}{(m_{b}-m_{s})^{2}}f_{0}^{2}(q^{2})\,,$
(24)

where $\mathcal{J}_{V}=(3-2\beta_{V}^{2}+3\beta_{V}^{4})/4$. Again, for
$V=\bar{V}$, the replacement $h_{S}\to 2h_{S}$ must be made and the rate
multiplied by $1/2$ for the reduction in phase space.

Finally, the differential decay rate for the process $B\to K\Psi\bar{\Psi}$ is

$\displaystyle\frac{d\Gamma(B\to K\Psi\bar{\Psi})}{dq^{2}}$
$\displaystyle=\frac{\beta_{\Psi}}{216\pi^{3}}\frac{|\vec{p}_{K}|q^{6}}{m_{\Psi}^{4}\Lambda^{4}}$
$\displaystyle\times\Bigg{[}\bigg{(}\Big{(}|F_{VV}|^{2}\mathcal{J}_{VV}+|F_{VA}|^{2}\beta_{\Psi}^{2}\mathcal{J}_{VA}\Big{)}f_{+}(q^{2})$
$\displaystyle+2\,\Re\Big{[}F_{VV}\Big{(}6F_{TT}^{*}\mathcal{J}_{VV,TT}+2F_{TS}^{*}\beta_{\Psi}^{4}+F_{\tilde{T}P}^{*}\mathcal{J}_{VV,\tilde{T}P}\Big{)}\Big{]}\frac{m_{\Psi}}{m_{B}+m_{K}}f_{T}(q^{2})\bigg{)}\frac{|\vec{p}_{K}|^{2}}{q^{2}}f_{+}(q^{2})$
$\displaystyle+\frac{3}{8}\bigg{(}|F_{VA}|^{2}(1-\beta_{\Psi}^{2})\mathcal{J}_{\Psi}\frac{(m_{b}-m_{s})^{2}}{q^{2}}$
$\displaystyle+|F_{SS}|^{2}\beta_{\Psi}^{2}\mathcal{J}_{\Psi}^{\prime}+|F_{SP}|^{2}\mathcal{J}_{\Psi}+4\,\Re\big{[}F_{VA}F_{SP}^{*}\big{]}\mathcal{J}_{\Psi}\frac{m_{\Psi}(m_{b}-m_{s})}{q^{2}}\bigg{)}\frac{m_{B}^{2}\delta_{K}^{2}}{(m_{b}-m_{s})^{2}}f_{0}^{2}(q^{2})$
$\displaystyle+2\bigg{(}|F_{TT}|^{2}\mathcal{J}_{TT}+\frac{1}{4}|F_{TS}|^{2}\beta_{\Psi}^{4}\mathcal{J}_{TS}+\frac{1}{4}|F_{TP}|^{2}\beta_{\Psi}^{2}\mathcal{J}_{TP}$
$\displaystyle+|F_{\tilde{T}T}|^{2}\beta_{\Psi}^{2}\mathcal{J}_{\Psi}+\frac{3}{16}|F_{\tilde{T}S}|^{2}\beta_{\Psi}^{2}(1-\beta_{\Psi}^{2})\mathcal{J}_{\tilde{T}S}+\frac{5}{16}|F_{\tilde{T}P}|^{2}\beta_{\Psi}^{\prime
2}(1-\beta_{\Psi}^{2})$
$\displaystyle+\Re\bigg{[}F_{TT}\bigg{(}F_{TS}^{*}\beta_{\Psi}^{4}\mathcal{J}_{TT,TS}+\frac{5}{4}F_{\tilde{T}P}^{*}\beta_{\Psi}^{\prime
2}(1-\beta_{\Psi}^{2})\bigg{)}-\frac{1}{2}F_{TP}F_{\tilde{T}S}^{*}\beta_{\Psi}^{2}(1-\beta_{\Psi}^{2})$
$\displaystyle+F_{\tilde{T}T}\bigg{(}F_{TP}^{*}\beta_{\Psi}^{2}\mathcal{J}_{\tilde{T}T,TP}-\frac{1}{4}F_{\tilde{T}S}^{*}\beta_{\Psi}^{2}(1-\beta_{\Psi}^{2})\mathcal{J}_{\tilde{T}T,\tilde{T}S}\bigg{)}\bigg{]}\bigg{)}\frac{|\vec{p}_{K}|^{2}}{(m_{B}+m_{K})^{2}}f_{T}^{2}(q^{2})\Bigg{]}\,,$
(25)

where we have defined for convenience the factors
$\mathcal{J}_{\Psi}=(9-6\beta_{\Psi}^{2}+5\beta_{\Psi}^{4})/8$,
$\mathcal{J}_{\Psi}^{\prime}=(5-6\beta_{\Psi}^{2}+9\beta_{\Psi}^{4})/8$ and,

$\displaystyle\mathcal{J}_{VV}$
$\displaystyle=(15-15\beta_{\Psi}^{2}+25\beta_{\Psi}^{4}-9\beta_{\Psi}^{6})/16\,,$
$\displaystyle\mathcal{J}_{VA}$
$\displaystyle=(5-2\beta_{\Psi}^{2}+5\beta_{\Psi}^{4})/8\,,$
$\displaystyle\mathcal{J}_{TT}$
$\displaystyle=(15-20\beta_{\Psi}^{2}+23\beta_{\Psi}^{4}-10\beta_{\Psi}^{6})/8\,,$
$\displaystyle\mathcal{J}_{TS}$
$\displaystyle=\mathcal{J}_{TT,TS}=(9-5\beta_{\Psi}^{2})/4\,,$
$\displaystyle\mathcal{J}_{TP}$
$\displaystyle=(7-3\beta_{\Psi}^{2})/4\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \mathcal{J}_{\tilde{T}S}=(11-5\beta_{\Psi}^{2})/6\,,$
$\displaystyle\mathcal{J}_{VV,TT}$
$\displaystyle=(15-10\beta_{\Psi}^{2}+19\beta_{\Psi}^{4})/24\,,$
$\displaystyle\mathcal{J}_{VV,\tilde{T}P}$
$\displaystyle=(15-10\beta_{\Psi}^{2}+3\beta_{\Psi}^{4})/8\,,$
$\displaystyle\mathcal{J}_{\tilde{T}T,TP}$
$\displaystyle=(3+\beta_{\Psi}^{2})/4\,,$
$\displaystyle\mathcal{J}_{\tilde{T}T,\tilde{T}S}$
$\displaystyle=(-3+5\beta_{\Psi}^{2})/2\,.$ (26)

For $\Psi=\Psi^{c}$, we get the decay rate from Eq. (B.2) with
$F_{VV},F_{AV},F_{TT},F_{\tilde{T}T}\to 0$ and $F_{XY}\to 2F_{XY}$ for the
remaining couplings. Eq. (B.2) must also be multiplied by $1/2$.

###  B.3 Invisible $B_{s}$ Decays

In the presence of light hidden states $X\in\\{\phi,\psi,V,\Psi\\}$, invisible
$B_{s}$ decays, $B_{s}\to\sum X$, can also occur via the same operators that
trigger $B\to K^{*}\sum X$. Here, we give expressions for the induced decay
rates, which are used to recast the upper bound on
$\mathcal{B}(B_{s}\to\text{inv})$ from ALEPH onto the effective coupling
parameter space. In the following, we make use of the $B_{s}$ decay matrix
elements defined in Appendix C.3.

The decay rate for a pair of outgoing invisible scalars is

$\displaystyle\Gamma(B_{s}\to\phi\bar{\phi})=\frac{\beta_{\phi}}{16\pi}\frac{f_{B_{s}}^{2}m_{B_{s}}}{\Lambda^{2}}\bigg{[}|g_{PS}|^{2}\frac{m_{B_{s}}^{2}}{(m_{b}+m_{s})^{2}}+\frac{1}{4}|g_{AV}|^{2}\frac{m_{B_{s}}^{2}}{\Lambda^{2}}\bigg{]}\,,$
(27)

where $\beta_{X}=\sqrt{1-4m_{X}^{2}/m_{B_{s}}^{2}}$. Again, $g_{AV}$ vanishes
if $\phi=\bar{\phi}$.

The decay rate for outgoing invisible fermions is given by,

$\displaystyle\Gamma(B_{s}\to\psi\bar{\psi})=\frac{\beta_{\psi}}{8\pi}\frac{f_{B_{s}}^{2}m_{B_{s}}^{3}}{\Lambda^{4}}\bigg{[}$
$\displaystyle
4|f_{AA}|^{2}\frac{m_{\psi}^{2}(m_{b}+m_{s})^{2}}{m_{B_{s}}^{4}}+|f_{PS}|^{2}\beta_{\psi}^{2}+|f_{PP}|^{2}$
$\displaystyle-4\,\Re\big{[}f_{AA}f_{PP}^{*}\big{]}\frac{m_{\psi}(m_{b}+m_{s})}{m_{B_{s}}^{2}}\bigg{]}\frac{m_{B_{s}}^{2}}{(m_{b}+m_{s})^{2}}\,.$
(28)

Likewise, for a pair of outgoing vector bosons the decay rate is

$\displaystyle\Gamma(B_{s}\to
V\bar{V})=\frac{\beta_{V}|h_{P}|^{2}}{64\pi}\frac{f_{B_{s}}^{2}m_{B_{s}}^{5}}{m_{V}^{4}\Lambda^{2}}\mathcal{J}_{V}\frac{m_{B_{s}}^{2}}{(m_{b}+m_{s})^{2}}\,,$
(29)

with $\mathcal{J}_{V}=(3-2\beta_{V}^{2}+3\beta_{V}^{4})/4$.

Finally, the decay rate for a pair of final state spin 3/2 fermions is

$\displaystyle\Gamma(B_{s}\to\Psi\bar{\Psi})=\frac{\beta_{\Psi}}{72\pi}\frac{f_{B_{s}}^{2}m_{B_{s}}^{7}}{m_{\Psi}^{4}\Lambda^{4}}\bigg{[}$
$\displaystyle
4|F_{AA}|^{2}\mathcal{J}_{\Psi}\frac{m_{\Psi}^{2}(m_{b}+m_{s})^{2}}{m_{B_{s}}^{4}}+|F_{PS}|^{2}\beta_{\Psi}^{2}\mathcal{J}_{\Psi}^{\prime}+|F_{PP}|^{2}\mathcal{J}_{\Psi}$
$\displaystyle-4\,\Re\big{[}F_{AA}F_{PP}^{*}\big{]}\mathcal{J}_{\Psi}\frac{m_{\Psi}(m_{b}+m_{s})}{m_{B_{s}}^{2}}\bigg{]}\frac{m_{B_{s}}^{2}}{(m_{b}+m_{s})^{2}}\,,$
(30)

where the $\mathcal{J}_{\Psi}$ and $\mathcal{J}_{\Psi}^{\prime}$ factors are
the same as those below Eq. (B.2). Again, all of the rates above must be
multiplied by a factor $1/2$ in the case of identical outgoing states
($\phi=\phi^{\dagger}$, $\psi=\psi^{c}$, $V=V^{\dagger}$, $\Psi=\Psi^{c}$),
and the couplings re-scaled in the same way as the previous sub-section.

###  B.4 Standard Model Rates

In the SM, the signal $B\to K^{(*)}E_{\text{miss}}$ is induced by the process
$B\to K^{(*)}\nu\bar{\nu}$, with the effective Hamiltonian,

$\displaystyle\mathcal{H}_{\text{eff}}^{\text{SM}}=-\frac{4G_{F}}{\sqrt{2}}\lambda_{t}C_{L}^{\text{SM}}\mathcal{O}_{L}+\text{h.c.}\,,$
(31)

where $\lambda_{t}=V_{tb}V_{ts}^{*}$ and $C_{L}^{\text{SM}}=-X_{t}/s_{w}^{2}$,
with the factor $X_{t}=1.469\pm 0.017$ accounting for NLO QCD Buchalla and
Buras (1994); Misiak and Urban (1999); Buchalla and Buras (1999) and two-loop
electroweak Brod  _et al._ (2011) corrections, and the left-handed operator is

$\displaystyle\mathcal{O}_{L}=\frac{e^{2}}{16\pi^{2}}(\bar{b}\gamma_{\mu}P_{L}s)(\bar{\nu}_{\alpha}\gamma^{\mu}(1-\gamma_{5})\nu_{\alpha})\,.$
(32)

In the operator basis of this work, this corresponds to

$\displaystyle
C^{V,LL}_{\underset{sb\alpha\beta}{d\nu}}=\frac{4G_{F}\lambda_{t}}{\sqrt{2}}\frac{\alpha}{2\pi}\frac{X_{t}}{s_{w}^{2}}\delta_{\alpha\beta}\,,$
(33)

which can be inserted into Eqs. (13), (B.2) and (B.2) to obtain the
differential rates

$\displaystyle\frac{d\Gamma_{\text{SM}}}{dq^{2}}=\frac{G_{F}^{2}\alpha^{2}\left|\lambda_{t}\right|^{2}X_{t}^{2}}{32\pi^{5}s_{w}^{4}}\left|\vec{p}_{K^{(*)}}\right|q^{2}f_{K^{(*)}}^{2}(q^{2})\,,$
(34)

where the $K^{(*)}$ three-momentum is given by Eq. (17) with $m_{\phi}^{2}\to
q^{2}$. In Eq. (34), we have defined the relevant hadronic form factor
combinations

$\displaystyle f_{K}^{2}(q^{2})$
$\displaystyle\equiv\frac{\left|\vec{p}_{K}\right|^{2}}{q^{2}}f_{+}^{2}(q^{2})\,,$
(35) $\displaystyle f_{K^{*}}^{2}(q^{2})$
$\displaystyle\equiv\frac{2|\vec{p}_{K^{*}}|^{2}}{(m_{B}+m_{K^{*}})^{2}}V^{2}(q^{2})+\frac{(m_{B}+m_{K^{*}})^{2}}{2m_{B}^{2}}A_{1}^{2}(q^{2})$
$\displaystyle\quad+\frac{16m_{K^{*}}^{2}}{q^{2}}A_{12}^{2}(q^{2})\,,$ (36)

where the form factors $f_{+}(q^{2})$, $V(q^{2})$, $A_{1}(q^{2})$ and
$A_{12}(q^{2})$ are defined in Appendix C.

On the other hand, the invisible $B_{s}$ decay width is negligible in the SM.
The two-body $B_{s}\to\nu\bar{\nu}$ rate is suppressed by tiny neutrino
masses, so the total invisible width is dominated by the four-body decay
$B_{s}\to\nu\bar{\nu}\nu\bar{\nu}$, yielding $\mathcal{B}(B_{s}\to E_{\rm
miss})\simeq 5\times 10^{-15}$ Bhattacharya  _et al._ (2019).

##  Appendix C Hadronic Matrix Elements

The latest form factors for both the $B\to K$ and $B\to K^{*}$ transitions can
be found in Ref. Gubernari  _et al._ (2023), where the authors perform a
dispersive analysis of both sets of form factors. The most up-to-date results
for the $B\to K$ form factors coming from the lattice are in Ref. Parrott  _et
al._ (2023), while in the case of $B\to K^{*}$ a determination on the lattice
of the full set of form factors can be found in Refs. Horgan  _et al._ (2014,
2015), which can be complemented by the LCSR determination in Ref. Gubernari
_et al._ (2019).

###  C.1 $B\to P$ Form Factors

The relevant hadronic matrix elements for the decay of a $B$ meson to a
pseudoscalar meson $P$, i.e. $B(p)\to P(k)$, can be written as follows:

$\displaystyle\bra{P}\bar{q}b\ket{B}$
$\displaystyle=\frac{m_{B}^{2}-m_{P}^{2}}{m_{b}-m_{q}}f_{0}(q^{2})\,,$ (37)
$\displaystyle\bra{P}\bar{q}\gamma_{\mu}b\ket{B}$
$\displaystyle=\bigg{[}P_{\mu}-\frac{m_{B}^{2}-m_{P}^{2}}{q^{2}}q_{\mu}\bigg{]}f_{+}(q^{2})+\frac{m_{B}^{2}-m_{P}^{2}}{q^{2}}q_{\mu}f_{0}(q^{2})\,,$
(38) $\displaystyle\bra{P}\bar{q}\sigma_{\mu\nu}b\ket{B}$
$\displaystyle=\frac{i(P_{\mu}q_{\nu}-P_{\nu}q_{\mu})}{m_{B}+m_{P}}f_{T}(q^{2})\,,$
(39)

where $P_{\mu}=p_{\mu}+k_{\mu}$, $q_{\mu}=p_{\mu}-k_{\mu}$, and the relation
$f_{+}(0)=f_{0}(0)$ holds. The form factors for the quark currents
$\bar{q}\gamma_{5}b$ and $\bar{q}\gamma_{\mu}\gamma_{5}b$ vanish identically,
while the form factor for $\bar{q}\sigma_{\mu\nu}\gamma_{5}b$ can be obtained
from the tensor form factor in Eq. (39) via the identity
$\sigma_{\mu\nu}\gamma_{5}=\frac{i}{2}\varepsilon_{\mu\nu\alpha\beta}\sigma^{\alpha\beta}$.

###  C.2 $B\to V$ Form Factors

The relevant hadronic matrix elements for $B$ meson decay to a vector meson
$V$, i.e. $B(p)\to V(k,\epsilon^{*})$, can be expressed as

$\displaystyle\bra{V}\bar{q}\gamma_{5}b\ket{B}$
$\displaystyle=-\frac{2im_{V}\epsilon^{*}_{\mu}q^{\mu}}{m_{b}+m_{q}}A_{0}(q^{2})\,,$
(40) $\displaystyle\bra{V}\bar{q}\gamma_{\mu}b\ket{B}$
$\displaystyle=\frac{2\varepsilon_{\mu\nu\alpha\beta}\epsilon^{*\nu}p^{\alpha}k^{\beta}}{m_{B}+m_{V}}V(q^{2})\,,$
(41) $\displaystyle\bra{V}\bar{q}\gamma_{\mu}\gamma_{5}b\ket{B}$
$\displaystyle=i\epsilon^{*\nu}\bigg{[}\frac{2m_{V}q_{\mu}q_{\nu}}{q^{2}}\big{(}A_{0}(q^{2})-A_{3}(q^{2})\big{)}+(m_{B}+m_{V})g_{\mu\nu}A_{1}(q^{2})-\frac{P_{\mu}q_{\nu}}{m_{B}+m_{V}}A_{2}(q^{2})\bigg{]}$
(42) $\displaystyle\bra{V}\bar{q}\sigma_{\mu\nu}b\ket{B}$
$\displaystyle=-i\varepsilon_{\mu\nu\alpha\beta}\epsilon^{*}_{\rho}\bigg{[}\bigg{(}P^{\alpha}-\frac{m_{B}^{2}-m_{V}^{2}}{q^{2}}q^{\alpha}\bigg{)}g^{\beta\rho}T_{1}(q^{2})+\frac{m_{B}^{2}-m_{V}^{2}}{q^{2}}q^{\alpha}g^{\beta\rho}T_{2}(q^{2})$
$\displaystyle\hskip
60.00009pt+\frac{q^{\alpha}P^{\beta}}{q^{2}}q^{\rho}\big{(}T_{1}(q^{2})-\widetilde{T}_{3}(q^{2})\big{)}\bigg{]}\,,$
(43)

where $A_{0}(0)=A_{3}(0)$, $T_{1}(0)=T_{2}(0)$, and

$\displaystyle\widetilde{T}_{3}(q^{2})\equiv
T_{2}(q^{2})+\frac{q^{2}}{m_{B}^{2}-m_{V}^{2}}T_{3}(q^{2})\,.$ (44)

The form factor for the scalar quark current $\bar{b}q$ vanishes identically,
while the form factor for $\bar{b}\sigma_{\mu\nu}\gamma_{5}q$ can again be
obtained from Eq. (43) via
$\sigma_{\mu\nu}\gamma_{5}=\frac{i}{2}\varepsilon_{\mu\nu\alpha\beta}\sigma^{\alpha\beta}$.

The form factor $A_{3}(q^{2})$ can be eliminated using the relation

$\displaystyle
A_{3}(q^{2})=\frac{m_{B}+m_{V}}{2m_{V}}A_{1}(q^{2})-\frac{m_{B}-m_{V}}{2m_{V}}A_{2}(q^{2})\,.$
(45)

The form factors $A_{2}(q^{2})$ and $T_{3}(q^{2})$ can likewise be eliminated
in favour of the so-called helicity form factors,

$\displaystyle A_{12}(q^{2})$
$\displaystyle\equiv\frac{(m_{B}+m_{V})^{2}(m_{B}^{2}-m_{V}^{2}-q^{2})A_{1}(q^{2})-\lambda(m_{B}^{2},q^{2},m_{V}^{2})A_{2}(q^{2})}{16m_{B}m_{V}^{2}(m_{B}+m_{V})}\,,$
(46) $\displaystyle T_{23}(q^{2})$
$\displaystyle\equiv\frac{(m_{B}^{2}-m_{V}^{2})(m_{B}^{2}+3m_{V}^{2}-q^{2})T_{2}(q^{2})-\lambda(m_{B}^{2},q^{2},m_{V}^{2})T_{3}(q^{2})}{8m_{B}m_{V}^{2}(m_{B}-m_{V})}\,,$
(47)

which appear in decay rates for the longitudinal polarisation of the vector
meson.

The required form factors for $B\to K^{(*)}$ decays are therefore $f_{+}$,
$f_{0}$, $f_{T}$, $V$, $A_{0}$, $A_{1}$, $A_{12}$, $T_{1}$, $T_{2}$ and
$T_{23}$. For these, we use the BSZ parametrisation Bharucha  _et al._ (2016)
results of Ref. Gubernari  _et al._ (2023).

###  C.3 $B_{s}$ Decay Matrix Elements

The relevant matrix elements for $B_{s}$ decay, which are used to derive the
decay rates $B_{s}\to\phi\bar{\phi}/\psi\bar{\psi}/V\bar{V}/\Psi\bar{\Psi}$ in
Appendix B, are

$\displaystyle\bra{0}\bar{s}\gamma_{\mu}\gamma_{5}b\ket{B_{s}}=if_{B_{s}}P_{\mu}\,,$
$\displaystyle\bra{0}\bar{s}\gamma_{5}b\ket{B_{s}}=-i\frac{m_{B_{s}}^{2}f_{B_{s}}}{m_{b}+m_{s}}\,,$
(48)

where $P_{\mu}$ is the four-momentum of $B_{s}$ and we take the value
$f_{B_{s}}=230.3(1.3)\leavevmode\nobreak\ \text{MeV}$ for the $B_{s}$ decay
constant Aoki  _et al._ (2022).

##  Appendix D Experimental Likelihood Reconstruction

In this appendix, we give further details on the reconstruction of the
likelihoods for the Belle II and BaBar analyses Adachi  _et al._ (2023); Lees
_et al._ (2013).

The Belle II Adachi  _et al._ (2023) analysis measured the differential decay
width for the charged channel $B^{+}\to K^{+}E_{\rm miss}$ with both inclusive
(ITA) and hadronic (HTA) tag methods. In order to calculate the expected
signal events using Eq. (6), we have been provided with the smearing
$f_{q^{2}_{\rm rec}}(q^{2})$ of the reconstructed momentum transfer
$q^{2}_{\rm rec}$ for the ITA analysis, and with the efficiencies for both ITA
and HTA analyses Adachi  _et al._ . We computed the Monte-Carlo statistical
errors for each background component ($B^{+}B^{-}$, $B^{0}\bar{B}^{0}$ and
continuum) as $\sigma^{i}_{b,{\rm stat}}=\sqrt{b^{i}/F_{\rm MC}}$, where
$b^{i}$ corresponds to the number of background events in the bin $i$ and
$F_{\rm MC}$ is a bin independent Monte-Carlo factor that reproduces the
simulation statistical uncertainties in Fig. 17 of Ref. Adachi  _et al._
(2023). As discussed in the main text, we performed a Monte-Carlo simulation
of the expected SM signal including uncertainties on the efficiencies and form
factors, to find the covariance $\Sigma_{\text{SM}}$, including correlations
between bins. The covariances for the background components were then found by
re-scaling $\Sigma_{\text{SM}}$ according to the background statistical errors
$\sigma^{i}_{b,{\rm stat}}$. We furthermore introduce a systematic uncertainty
for the normalisation of the backgrounds, which are fitted to reproduce the
profile log-likelihood ratio as a function of the signal strength given in
Fig. 16 of Ref. Adachi  _et al._ (2023).

The BaBar Lees  _et al._ (2013) analysis performed a search for both the
neutral and charged channels $B^{0(+)}\to K^{*0(+)}E_{\rm miss}$, combining
two different final states for each channel ($K^{+}\pi^{-}$, $K_{S}\pi^{0}$,
$K^{+}\pi^{0}$ and $K_{S}\pi^{+}$) and using a hadronic tag method. We extract
the different background contributions and their $q^{2}$-dependence from Fig.
5 of Ref. Lees  _et al._ (2013), while we extract the efficiencies
$\epsilon(q^{2})$ from Fig. 6 in their binned forms. We compute the
statistical errors similarly to the Belle II analysis, fitting the Monte-Carlo
factor to reproduce the values in Table IV of Ref. Lees  _et al._ (2013). We
obtained the systematic uncertainties from Table II and III of Ref. Lees  _et
al._ (2013), taking into account the correlation between the systematic
uncertainties in the background and efficiency.

Since no information is available on the smearing $f_{q^{2}_{\rm rec}}(q^{2})$
of the reconstructed momentum transfer $q^{2}_{\rm rec}$, we exclude the
constraints obtained from $B^{+}\to K^{*+}E_{\rm miss}$, as both final states
contain neutral particles leading to non-negligible and asymmetrical smearing.
In the case $B^{0}\to K^{*0}E_{\rm miss}$ we assume that the $K^{*0}$
dominantly decays to $K^{+}\pi^{-}$ (as there is a factor of 10 in the
efficiency) and we therefore neglect the smearing.

Lastly, we have verified that our likelihoods reproduce correctly the
confidence intervals and limits given in both Refs. Adachi  _et al._ (2023);
Lees  _et al._ (2013).

## References

  * Amhis  _et al._ (2023) Y. S. Amhis _et al._ (HFLAV), Phys. Rev. D 107, 052008 (2023), arXiv:2206.07501 [hep-ex] .
  * Algueró  _et al._ (2023) M. Algueró, A. Biswas, B. Capdevila, S. Descotes-Genon, J. Matias, and M. Novoa-Brunet, Eur. Phys. J. C 83, 648 (2023), arXiv:2304.07330 [hep-ph] .
  * Adachi  _et al._ (2023) I. Adachi _et al._ (Belle-II), (2023), arXiv:2311.14647 [hep-ex] .
  * Lees  _et al._ (2013) J. P. Lees _et al._ (BaBar), Phys. Rev. D 87, 112005 (2013), arXiv:1303.7465 [hep-ex] .
  * Alonso-Álvarez and Escudero (2023) G. Alonso-Álvarez and M. Escudero, (2023), arXiv:2310.13043 [hep-ph] .
  * Barate  _et al._ (2001) R. Barate _et al._ (ALEPH), Eur. Phys. J. C 19, 213 (2001), arXiv:hep-ex/0010022 .
  * Buras (2020) A. Buras,  _Gauge Theory of Weak Decays_ (Cambridge University Press, 2020).
  * Buras  _et al._ (2015) A. J. Buras, J. Girrbach-Noe, C. Niehoff, and D. M. Straub, JHEP 02, 184 (2015), arXiv:1409.4557 [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] .
  * Gubernari  _et al._ (2023) N. Gubernari, M. Reboud, D. van Dyk, and J. Virto, JHEP 12, 153 (2023), arXiv:2305.06301 [hep-ph] .
  * Athron  _et al._ (2023) P. Athron, R. Martinez, and C. Sierra, (2023), arXiv:2308.13426 [hep-ph] .
  * Bause  _et al._ (2024) R. Bause, H. Gisbert, and G. Hiller, Phys. Rev. D 109, 015006 (2024), arXiv:2309.00075 [hep-ph] .
  * Allwicher  _et al._ (2024) L. Allwicher, D. Becirevic, G. Piazza, S. Rosauro-Alcaraz, and O. Sumensari, Phys. Lett. B 848, 138411 (2024), arXiv:2309.02246 [hep-ph] .
  * Chen and Chiang (2023) C.-H. Chen and C.-W. Chiang, (2023), arXiv:2309.12904 [hep-ph] .
  * Felkl  _et al._ (2023) T. Felkl, A. Giri, R. Mohanta, and M. A. Schmidt, Eur. Phys. J. C 83, 1135 (2023), arXiv:2309.02940 [hep-ph] .
  * Abdughani and Reyimuaji (2023) M. Abdughani and Y. Reyimuaji, (2023), arXiv:2309.03706 [hep-ph] .
  * Dreiner  _et al._ (2023) H. K. Dreiner, J. Y. Günther, and Z. S. Wang, (2023), arXiv:2309.03727 [hep-ph] .
  * He  _et al._ (2023) X.-G. He, X.-D. Ma, and G. Valencia, (2023), arXiv:2309.12741 [hep-ph] .
  * Berezhnoy and Melikhov (2023) A. Berezhnoy and D. Melikhov, (2023), arXiv:2309.17191 [hep-ph] .
  * Datta  _et al._ (2023) A. Datta, D. Marfatia, and L. Mukherjee, (2023), arXiv:2310.15136 [hep-ph] .
  * Altmannshofer  _et al._ (2023) W. Altmannshofer, A. Crivellin, H. Haigh, G. Inguglia, and J. Martin Camalich, (2023), arXiv:2311.14629 [hep-ph] .
  * Fridell  _et al._ (2023) K. Fridell, M. Ghosh, T. Okui, and K. Tobioka, (2023), arXiv:2312.12507 [hep-ph] .
  * Gabrielli  _et al._ (2024) E. Gabrielli, L. Marzola, K. Müürsepp, and M. Raidal, (2024), arXiv:2402.05901 [hep-ph] .
  * Chen  _et al._ (2024) F.-Z. Chen, Q. Wen, and F. Xu, (2024), arXiv:2401.11552 [hep-ph] .
  * Hou  _et al._ (2024) B.-F. Hou, X.-Q. Li, M. Shen, Y.-D. Yang, and X.-B. Yuan, (2024), arXiv:2402.19208 [hep-ph] .
  * Descotes-Genon  _et al._ (2020) S. Descotes-Genon, S. Fajfer, J. F. Kamenik, and M. Novoa-Brunet, Phys. Lett. B 809, 135769 (2020), [Addendum: Phys.Lett.B 840, 137830 (2023)], arXiv:2005.03734 [hep-ph] .
  * Kamenik and Smith (2012) J. F. Kamenik and C. Smith, JHEP 03, 090 (2012), arXiv:1111.6402 [hep-ph] .
  * Dror  _et al._ (2017) J. A. Dror, R. Lasenby, and M. Pospelov, Phys. Rev. D 96, 075036 (2017), arXiv:1707.01503 [hep-ph] .
  * Di Luzio  _et al._ (2022) L. Di Luzio, M. Nardecchia, and C. Toni, Phys. Rev. D 105, 115042 (2022), arXiv:2204.05945 [hep-ph] .
  * Dembinski and et al. (2020) H. Dembinski and P. O. et al.,  (2020), 10.5281/zenodo.3949207.
  * James and Roos (1975) F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
  * He  _et al._ (2024) X.-G. He, X.-D. Ma, M. A. Schmidt, G. Valencia, and R. R. Volkas, (2024), arXiv:2403.12485 [hep-ph] .
  * Badin and Petrov (2010) A. Badin and A. A. Petrov, Phys. Rev. D 82, 034005 (2010), arXiv:1005.1277 [hep-ph] .
  * Bauer  _et al._ (2022) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 09, 056 (2022), arXiv:2110.10698 [hep-ph] .
  * Li and Tandean (2023) G. Li and J. Tandean, JHEP 11, 205 (2023), arXiv:2306.05333 [hep-ph] .
  * Cortina Gil  _et al._ (2021) E. Cortina Gil _et al._ (NA62), JHEP 06, 093 (2021), arXiv:2103.15389 [hep-ex] .
  * Ahn  _et al._ (2021) J. K. Ahn _et al._ (KOTO), Phys. Rev. Lett. 126, 121801 (2021), arXiv:2012.07571 [hep-ex] .
  * Ablikim  _et al._ (2022) M. Ablikim _et al._ (BESIII), Phys. Rev. D 105, L071102 (2022), arXiv:2112.14236 [hep-ex] .
  * Khachatryan  _et al._ (2015) V. Khachatryan _et al._ (CMS), Phys. Rev. Lett. 114, 101801 (2015), arXiv:1410.1149 [hep-ex] .
  * Aad  _et al._ (2024) G. Aad _et al._ (ATLAS), (2024), arXiv:2403.02793 [hep-ex] .
  * Bauer  _et al._ (2021) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 04, 063 (2021), arXiv:2012.12272 [hep-ph] .
  * Buchalla and Buras (1994) G. Buchalla and A. J. Buras, Nucl. Phys. B 412, 106 (1994), arXiv:hep-ph/9308272 .
  * Misiak and Urban (1999) M. Misiak and J. Urban, Phys. Lett. B 451, 161 (1999), arXiv:hep-ph/9901278 .
  * Buchalla and Buras (1999) G. Buchalla and A. J. Buras, Nucl. Phys. B 548, 309 (1999), arXiv:hep-ph/9901288 .
  * Brod  _et al._ (2011) J. Brod, M. Gorbahn, and E. Stamou, Phys. Rev. D 83, 034030 (2011), arXiv:1009.0947 [hep-ph] .
  * Bhattacharya  _et al._ (2019) B. Bhattacharya, C. M. Grant, and A. A. Petrov, Phys. Rev. D 99, 093010 (2019), arXiv:1809.04606 [hep-ph] .
  * Parrott  _et al._ (2023) W. G. Parrott, C. Bouchard, and C. T. H. Davies ((HPQCD collaboration)§, HPQCD), Phys. Rev. D 107, 014510 (2023), arXiv:2207.12468 [hep-lat] .
  * Horgan  _et al._ (2014) R. R. Horgan, Z. Liu, S. Meinel, and M. Wingate, Phys. Rev. D 89, 094501 (2014), arXiv:1310.3722 [hep-lat] .
  * Horgan  _et al._ (2015) R. R. Horgan, Z. Liu, S. Meinel, and M. Wingate, PoS LATTICE2014, 372 (2015), arXiv:1501.00367 [hep-lat] .
  * Gubernari  _et al._ (2019) N. Gubernari, A. Kokulu, and D. van Dyk, JHEP 01, 150 (2019), arXiv:1811.00983 [hep-ph] .
  * Bharucha  _et al._ (2016) A. Bharucha, D. M. Straub, and R. Zwicky, JHEP 08, 098 (2016), arXiv:1503.05534 [hep-ph] .
  * Aoki  _et al._ (2022) Y. Aoki _et al._ (Flavour Lattice Averaging Group (FLAG)), Eur. Phys. J. C 82, 869 (2022), arXiv:2111.09849 [hep-lat] .
  * (53) I. Adachi _et al._ (Belle-II), Private Communication.