aainstitutetext: Institute of Particle Physics and Key Laboratory of Quark and
Lepton Physics (MOE),
Central China Normal University, Wuhan, Hubei 430079, Chinabbinstitutetext:
Center for High Energy Physics, Peking University, Beijing 100871,
Chinaccinstitutetext: Institute of Particle and Nuclear Physics, Henan Normal
University, Xinxiang 453007, China

# Deciphering the Belle II data on $\boldsymbol{B\to K\nu\bar{\nu}}$ decay in
the (dark) SMEFT with minimal flavour violation

Biao-Feng Hou<EMAIL_ADDRESS>a,b Xin-Qiang Li<EMAIL_ADDRESS>a Meng
Shen<EMAIL_ADDRESS>a,c Ya-Dong Yang<EMAIL_ADDRESS>a and
Xing-Bo Yuan<EMAIL_ADDRESS>
###### Abstract

Recently, the Belle II collaboration announced the first measurement of the
branching ratio $\mathcal{B}(B^{+}\to K^{+}\nu\bar{\nu})$, which is found to
be about $2.7\sigma$ higher than the Standard Model (SM) prediction. We
decipher the data with two new physics scenarios: the underlying quark-level
$b\to s\nu\bar{\nu}$ transition is, besides the SM contribution, further
affected by heavy new mediators that are much heavier than the electroweak
scale, or amended by an additional decay channel with undetected light final
states like dark matter or axion-like particles. These two scenarios can be
most conveniently analyzed in the SM effective field theory (SMEFT) and the
dark SMEFT (DSMEFT) framework, respectively. We consider the flavour
structures of the resulting effective operators to be either generic or
satisfy the minimal flavour violation (MFV) hypothesis, both for the quark and
lepton sectors. In the first scenario, once the MFV assumption is made, only
one SM-like low-energy effective operator induced by the SMEFT dimension-six
operators can account for the Belle II excess, the parameter space of which
is, however, excluded by the Belle upper bound of the branching ratio
$\mathcal{B}(B^{0}\to K^{*0}\nu\bar{\nu})$. In the second scenario, it is
found that the Belle II excess can be accommodated by 22 of the DSMEFT
operators involving one or two scalar, fermionic, or vector dark matters as
well as axion-like particles. These operators also receive dominant
constraints from the $B^{0}\to K^{*0}+\rm inv$ and $B_{s}\to\rm inv$ decays.
Once the MFV hypothesis is assumed, the number of viable operators is reduced
to 14, and the $B^{+}\to\pi^{+}+\rm inv$ and $K^{+}\to\pi^{+}+\rm inv$ decays
start to put further constraints on them. Within the parameter space allowed
by all the current experimental data, the $q^{2}$ distributions of the $B\to
K^{(*)}+\rm inv$ decays are then studied for each viable operator. We find
that the resulting prediction of the operator
$\mathcal{Q}_{q\chi}=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}\big{(}\bar{\chi}\gamma^{\mu}\chi\big{)}$
with a fermionic dark matter mass $m_{\chi}\approx 700\,\rm MeV$ can closely
match the Belle II event distribution in the bins $2\leq q^{2}\leq 7\,{\rm
GeV}^{2}$. In addition, the future prospects at Belle II, CEPC and FCC-ee are
also discussed for some of these FCNC processes.

##  1 Introduction

The flavour-changing neutral current (FCNC) processes are highly suppressed by
the Glashow-Iliopoulos-Maiani (GIM) mechanism in the Standard Model (SM)
Glashow:1970gm , making them sensitive probes of new physics (NP) beyond the
SM. In particular, the $B\to K^{(\ast)}\nu\bar{\nu}$ decays are one of the
cleanest channels within the SM and most suitable for indirect NP searches
Altmannshofer:2009ma ; Buras:2014fpa ; Becirevic:2023aov  (see also refs.
Felkl:2021uxi ; Bause:2021cna ; He:2021yoz ; Browder:2021hbl  for recent
studies on these as well as other related decays in the context of NP, but
before the recent Belle II announcement Belle-II:2023esi ).

Recently, using the dataset with an integrated luminosity of $362{\,\rm
fb}^{-1}$, the Belle II collaboration reported the first evidence of the
$B^{+}\to K^{+}+\nu\bar{\nu}$ decay with a branching ratio Belle-II:2023esi

$\displaystyle{\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})_{\text{Belle II}}$
$\displaystyle=\left[23\pm 5(\text{stat})^{+5}_{-4}(\text{syst})\right]\times
10^{-6},$ (1)

which is obtained by combining the inclusive and hadronic tagging results.
Interestingly, this measurement shows about $2.7\sigma$ higher than the SM
prediction $(4.16\pm 0.57)\times 10^{-6}$, which is based on the calculations
in refs. Altmannshofer:2009ma ; Buras:2014fpa ; Becirevic:2023aov  but with
updated input parameters. In the future, the Belle II with just $5{\,\rm
ab}^{-1}$ integrated luminosity is expected to measure this decay with $20\sim
30\%$ uncertainty relative to the SM value Belle-II:2018jsg ; Belle-II:2022cgf
, and could therefore confirm or exclude such an excess. Nevertheless, it has
already motivated numerous studies of possible NP explanations Athron:2023hmz
; Bause:2023mfe ; Allwicher:2023xba ; Felkl:2023ayn ; Abdughani:2023dlr ;
Dreiner:2023cms ; He:2023bnk ; Berezhnoy:2023rxx ; Datta:2023iln ;
Altmannshofer:2023hkn ; McKeen:2023uzo ; Fridell:2023ssf ; Ho:2024cwk ;
Chen:2024jlj ; Gabrielli:2024wys ; Li:2024thq .

As an interesting feature, the Belle II result, once the $2\sigma$ range in
eq. (1) is considered, would be higher than the SM prediction by a factor of
$3\sim 8$. Such a large excess motivates us to ask the following questions:
why has such a large NP effect not shown up or will it be found in other $b\to
s$ processes, such as the $B\to K^{\ast}\nu\bar{\nu}$ and
$B_{s}\to\nu\bar{\nu}$ decays Workman:2022ynf ; Alonso-Alvarez:2023mgc ? In
this paper, we aim to answer these questions by performing a detailed model-
independent analysis of possible NP effects on these and related FCNC
processes in the Effective Field Theory (EFT) framework.

The $B^{+}\to K^{+}\nu\bar{\nu}$ excess reported by Belle II could be
explained by introducing heavy NP particles that are much heavier than the
electroweak scale and hence contribute to the $b\to s\nu\bar{\nu}$ transition
as virtual mediators. In this case, the Standard Model Effective Field Theory
(SMEFT) provides a well-suited model-independent framework to study such kinds
of NP effects Brivio:2017vri ; Isidori:2023pyp . Here, all the NP effects are
encoded in the Wilson coefficients of the higher-dimensional operators that
are invariant under the SM gauge group. Below the electroweak scale, the heavy
SM particles like the massive bosons and the top quark also decouple, and the
dynamics is then described by the Low-energy Effective Field Theory (LEFT)
Buchalla:1995vs ; Jenkins:2017jig , which consists of, besides the QCD and QED
Lagrangians for light SM fermions, a set of higher-dimensional operators
compatible with the QCD and QED gauge symmetries. The Wilson coefficients of
these higher-dimensional operators encode all the physics related to the heavy
SM and NP degrees of freedom. Starting with the SMEFT, one can also fix these
low-energy Wilson coefficients by perform a matching between these two EFTs at
the electroweak scale Jenkins:2017jig ; Dekens:2019ept ; Aebischer:2015fzz .
Such a procedure has been extensively applied to various $B$-meson and kaon
decays Buchalla:1995vs ; Jenkins:2017jig . In this paper, we will interpret
the Belle II excess in terms of the SMEFT operators and consider the LEFT as
the low-energy limit of SMEFT to calculate the relevant $B$-meson and kaon
decays.

The Belle II excess could also be accommodated by amending the $B^{+}\to
K^{+}\nu\bar{\nu}$ decay by an additional decay channel with undetected light
final states like dark matter (DM) or axion-like particles (ALP). To mimic the
signal, these new light states must couple to the SM quark sector and be
sufficiently long-lived to escape the detector. The potential candidates are
new light particles that are singlet under the SM gauge group. In order to
study the interactions of such new states with the SM quark sector model-
independently, it is also convenient to work in an EFT framework, such as the
DM EFT DelNobile:2011uf ; Baumgart:2022vwr  without any stabilizing symmetry
on the DM fields.111In this work, we do not impose any symmetry to stabilize
the DM fields in the EFT framework. Thus, the DM fields we are considering are
actually just new singlets under the SM gauge group. In this respect, the
SMEFT can be extended by including the scalar, fermionic or vector DM fields,
resulting in the so-called Dark SMEFT (DSMEFT) Duch:2014xda ; Kamenik:2011vy ;
Brod:2017bsw ; Criado:2021trs ; Aebischer:2022wnl ; Song:2023jqm . In the
DSMEFT, the new operators involving the DM fields are invariant under the SM
$SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}$ gauge group and the electroweak
symmetry is broken by the usual Higgs mechanism. The general non-redundant
operator basis can be found in ref. Aebischer:2022wnl  up to dim$=6$, and in
ref. Song:2023jqm  up to dim$=8$. In order to describe the low-energy
processes below the electroweak scale, such as $d_{i}\to d_{j}+{\rm DM}$
decays, a general framework is the generalization of the LEFT to include the
scalar, fermionic and vector DM fields, referred to as the Dark LEFT (DLEFT)
Aebischer:2022wnl . In the DLEFT, the DM particles should be lighter than the
electroweak scale and the effective operators are invariant under the
$SU(3)_{C}\otimes U(1)_{\rm em}$ gauge group. A complete basis of the DLEFT
operators can be found in refs. Aebischer:2022wnl ; He:2022ljo  for dim $\leq
6$, and in ref. Liang:2023yta  for dim $\leq 7$. In this work, starting from
the DSMEFT and considering the DLEFT as the low-energy limit of DSMEFT, we
investigate the parameter space of each DSMEFT operator required to explain
the Belle II excess, together under the constraints from the other $b\to
s+\text{inv}$ decays. The future prospects at Belle II, CEPC and FCC-ee will
also be discussed for some of these FCNC processes. In addition, for
completeness, our analyses also include the ALP case Peccei:1977hh ;
Peccei:1977ur ; Weinberg:1977ma ; Wilczek:1977pj  in the EFT framework
Georgi:1986df ; Brivio:2017ije ; Chala:2020wvs ; Bauer:2020jbp ; Galda:2021hbr
; Song:2023lxf .

The above two NP scenarios proposed to account for the Belle II excess could
also affect the quark-level $b\to d$ and $s\to d$ processes. In order to
correlate the $b\to s$ to the $b\to d$ and $s\to d$ transitions, we have to
specify the underlying flavour structures of the EFT considered Isidori:2010kg
; Altmannshofer:2022aml . To this end, we assume the minimal flavour violation
(MFV) hypothesis Chivukula:1987py ; Buras:2000dm ; DAmbrosio:2002vsn  and
implement it into the SMEFT and DSMEFT frameworks. According to the
hypothesis, all the flavour violating currents are controlled by the SM Yukawa
couplings, so that all the FCNC interactions in the quark sector are naturally
suppressed by the CKM factors. Therefore, any potential large FCNCs in the
SMEFT and DSMEFT can be avoided. Considering the current experimental data on
the $b\to s+\text{inv}$, $b\to d+\text{inv}$ and $s\to d+\text{inv}$
processes, such as the $B\to\pi+\text{inv}$ and $K\to\pi+\text{inv}$ decays,
we derive the constraints on each operator in the SMEFT and DSMEFT with the
MFV hypothesis. Implications of the resulting parameter space of each viable
operator for the differential distributions of these decays are also
investigated in detail, which can be tested at the Belle II with more
statistics.

The remainder of this paper is organized as follows. In section 2, we
introduce the MFV hypothesis. In sections 3 and 4, we investigate the various
$d_{i}\to d_{j}\nu\bar{\nu}$ and $d_{i}\to d_{j}+{\rm DM}$ decays in the SMEFT
and DSMEFT, respectively. In these EFT frameworks, the flavour structures of
the effective operators are taken to be either of the most general form or
satisfy the MFV hypothesis. We present our detailed numerical analyses and
discussions in section 5. Our conclusion is finally made in section 6. The
hadronic matrix elements involved throughout this work, such as the transition
form factors and decay constants, are summarized in appendix A.

##  2 Minimal Flavour Violation

In this section, we discuss the idea of MFV hypothesis and its application in
the SMEFT and DSMEFT.

###  2.1 Quark sector

In the interaction eigenbasis, the SM Yukawa interactions in the quark sector
are described by the Lagrangian

$\displaystyle-{\mathcal{L}}_{Y}=\bar{q}\,Y_{d}Hd+\bar{q}\,Y_{u}\tilde{H}u+\mathrm{h.c.},$
(2)

where $\tilde{H}\equiv i\sigma_{2}H^{*}$ with $H$ being the SM Higgs doublet.
The fields $q$ and $u,d$ denote the left-handed quark doublet and the right-
handed quark singlets in the SM, respectively. The Yukawa coupling matrices
$Y_{u,d}$ are $3\times 3$ complex matrices in flavour space. In the SM, these
Yukawa interactions violate the global flavour symmetry Chivukula:1987py

$\displaystyle G_{\rm QF}=SU(3)_{q}\otimes SU(3)_{u}\otimes SU(3)_{d}.$ (3)

We can formally recover this flavour symmetry by promoting the Yukawa matrices
$Y_{u,d}$ to spurion fields DAmbrosio:2002vsn . The Yukawa spurions should
transform under $G_{\rm QF}$ as

$\displaystyle
Y_{u}\sim\big{(}\boldsymbol{3},\boldsymbol{\bar{3}},\boldsymbol{1}\big{)},\qquad
Y_{d}\sim\big{(}\boldsymbol{3},\boldsymbol{1},\boldsymbol{\bar{3}}\big{)}.$
(4)

Then, we can construct two useful basic building blocks ${\mathsf{A}}\equiv
Y_{u}Y_{u}^{\dagger}$ and ${\mathsf{B}}\equiv Y_{d}Y_{d}^{\dagger}$, which
transform as
$(\boldsymbol{1}\oplus\boldsymbol{8},\boldsymbol{1},\boldsymbol{1})$ under the
flavour group $G_{\rm QF}$. As will be seen later, these building blocks are
useful to parameterize the quark flavour transitions in effective theories.

In the SMEFT and DSMEFT, there are four types of quark currents relevant for
the down-type FCNC transitions, which read

$\displaystyle\bar{q}\gamma^{\mu}{\mathcal{C}}q,\qquad\bar{d}\gamma^{\mu}{\mathcal{C}}d,\qquad\bar{q}{\mathcal{C}}d,\qquad\bar{q}\sigma^{\mu\nu}{\mathcal{C}}d,$
(5)

where ${\mathcal{C}}$ denotes the Wilson coefficient of the operator
${\mathcal{Q}}$ and is a $3\times 3$ matrix in quark flavour space. Here and
in the following, we suppress the operator index $i$, i.e., ${\mathcal{C}}$
denoting ${\mathcal{C}}_{i}$. In the MFV hypothesis, it is assumed that at low
energy the Yukawa couplings are the only source of flavour violation in the SM
and beyond DAmbrosio:2002vsn . Technically, this hypothesis implies that the
quark currents in eq. (5) should be invariant under the flavour group $G_{\rm
QF}$. Therefore, the Wilson coefficient ${\mathcal{C}}$ can be written in the
following form:

$\displaystyle{\mathcal{C}}^{\rm
MFV}=\begin{cases}f({\mathsf{A}},{\mathsf{B}})&\text{for
\;}\bar{q}\gamma^{\mu}{\mathcal{C}}q,\\\
f({\mathsf{A}},{\mathsf{B}})Y_{d}&\text{for
\;}\bar{q}{\mathcal{C}}d,\,\bar{q}\sigma^{\mu\nu}{\mathcal{C}}d,\\\
\epsilon_{0}\mathds{1}+Y_{d}^{\dagger}g({\mathsf{A}},{\mathsf{B}})Y_{d}&\text{for
\;}\bar{d}\gamma^{\mu}{\mathcal{C}}d,\end{cases}$ (6)

where $\epsilon_{0}$ is a real constant. The functions
$f({\mathsf{A}},{\mathsf{B}})$ and $g({\mathsf{A}},{\mathsf{B}})$ are infinite
series of the building blocks ${\mathsf{A}}$ and ${\mathsf{B}}$, e.g.,
explicitly
$f({\mathsf{A}},{\mathsf{B}})\equiv\xi_{ijk\dotsm}{\mathsf{A}}^{i}{\mathsf{B}}^{j}{\mathsf{A}}^{k}\dotsb$
with real coefficients $\xi_{ijk\dotsm}$ to ensure no new sources of CP
violation beyond the Yukawa couplings. By application of the Cayley-Hamilton
identity, the infinite series in $f({\mathsf{A}},{\mathsf{B}})$ and
$g({\mathsf{A}},{\mathsf{B}})$ can be resummed into a finite number of terms.
Taking $f({\mathsf{A}},{\mathsf{B}})$ as an example, its explicit summed
expression can be written as Colangelo:2008qp ; Mercolli:2009ns ;
Grinstein:2023njq

$\displaystyle f({\mathsf{A}},{\mathsf{B}})=\epsilon_{0}\mathds{1}+$
$\displaystyle\epsilon_{1}{\mathsf{A}}+\epsilon_{3}{\mathsf{A}}^{2}+\epsilon_{5}{\mathsf{A}}{\mathsf{B}}+\epsilon_{7}{\mathsf{A}}{\mathsf{B}}{\mathsf{A}}+\epsilon_{10}{\mathsf{A}}{\mathsf{B}}^{2}+\epsilon_{12}{\mathsf{A}}^{2}{\mathsf{B}}^{2}+\epsilon_{14}{\mathsf{B}}^{2}{\mathsf{A}}{\mathsf{B}}+\epsilon_{15}{\mathsf{A}}{\mathsf{B}}^{2}{\mathsf{A}}^{2}$
$\displaystyle+$
$\displaystyle\epsilon_{2}{\mathsf{B}}+\epsilon_{4}{\mathsf{B}}^{2}+\epsilon_{6}{\mathsf{B}}{\mathsf{A}}+\epsilon_{9}{\mathsf{B}}{\mathsf{A}}{\mathsf{B}}+\epsilon_{8}{\mathsf{B}}{\mathsf{A}}^{2}+\epsilon_{13}{\mathsf{B}}^{2}{\mathsf{A}}^{2}+\epsilon_{11}{\mathsf{A}}{\mathsf{B}}{\mathsf{A}}^{2}+\epsilon_{16}{\mathsf{B}}^{2}{\mathsf{A}}^{2}{\mathsf{B}}.$

During the resummation, higher-order polynomials of ${\mathsf{A}}$ and
${\mathsf{B}}$ contribute to the coefficients $\epsilon_{i}$, and make them
complex generally. However, it is found that the imaginary parts of the
coefficients are tiny (e.g.,
$\text{Im}\epsilon_{i}\propto|\text{Tr}({\mathsf{A}}^{2}{\mathsf{B}}{\mathsf{A}}{\mathsf{B}}^{2})|\ll
1$) and can be therefore neglected in the numerical analysis Colangelo:2008qp
; Mercolli:2009ns ; He:2014uya ; He:2014fva ; He:2014efa . In the following,
we assume the coefficients $\epsilon_{i}$ are of the same orders of magnitude.
Therefore, only the terms built by ${\mathsf{A}}$ are kept, while all the
terms involving ${\mathsf{B}}$ are neglected, as the spurion ${\mathsf{B}}$ is
suppressed by ${\mathcal{O}}(\hat{\lambda}_{d}^{2})$, where
$\hat{\lambda}_{d}$ denotes the diagonal Yukawa coupling matrix for down-type
quarks defined by eq. (9). Finally, we obtain the approximation Chiang:2017hlj

$\displaystyle
f({\mathsf{A}},{\mathsf{B}})\approx\epsilon_{0}\mathds{1}+\epsilon_{1}{\mathsf{A}}+\epsilon_{2}{\mathsf{A}}^{2},$
(7)

where the coefficients $\epsilon_{0,1,2}$ are free real parameters. Similarly,
being suppressed by ${\mathcal{O}}(\hat{\lambda}_{d}^{2}$), the approximation
$Y_{d}^{\dagger}g({\mathsf{A}},{\mathsf{B}})Y_{d}\approx 0$ can be taken.

In the MFV framework discussed above, the quark currents are all given in the
interaction eigenbasis. In order to rotate to the fermion mass eigenbasis, the
following transforms should be performed:

$\displaystyle u_{L}\to U_{u}u_{L},\qquad u_{R}\to W_{u}u_{R},\qquad d_{L}\to
U_{d}d_{L},\qquad d_{R}\to W_{d}d_{R},$ (8)

where $U_{u,d}$ and $W_{u,d}$ are $3\times 3$ unitary matrices introduced to
diagonalize the Yukawa matrices

$\displaystyle U_{d}^{\dagger}Y_{d}W_{d}=\hat{\lambda}_{d},\qquad
U_{u}^{\dagger}Y_{u}W_{u}=\hat{\lambda}_{u}.$ (9)

Therefore, the CKM matrix is given by $V=U_{u}^{\dagger}U_{d}$. In the
following, we turn to the fermion mass eigenbasis and derive explicit
expressions of the MFV couplings.

The current $\bar{q}\gamma^{\mu}{\mathcal{C}}q$ involves the following down-
type quark interactions in the mass eigenbasis:

$\displaystyle\bar{d}_{L}U_{d}^{\dagger}\gamma_{\mu}{\mathcal{C}}^{\rm
MFV}U_{d}d_{L}\approx\bar{d}_{L}\gamma^{\mu}V^{\dagger}\big{(}\epsilon_{0}+\epsilon_{1}\hat{\lambda}_{u}^{2}+\epsilon_{2}\hat{\lambda}_{u}^{4}\big{)}Vd_{L},$
(10)

where $d_{L}$ denotes now the left-handed down-type quark in the mass
eigenbasis, and $\hat{\lambda}_{u}$ denotes the diagonal Yukawa coupling
matrix for up-type quarks, whose diagonal elements are given by
$\lambda_{u_{i}}=\sqrt{2}m_{u_{i}}/v$ with the vacuum expectation value
$v=246{\,\rm GeV}$. It is noted that, due to the large hierarchy among the
diagonal elements $\lambda_{u_{i}}$, the matrices
$\epsilon_{1}\hat{\lambda}_{u}^{2}$ and $\epsilon_{2}\hat{\lambda}_{u}^{4}$
have almost the same structures. Therefore, one can obtain the following
approximation Zhang:2018nmy :

$\displaystyle
V^{\dagger}(\epsilon_{1}\hat{\lambda}_{u}^{2}+\epsilon_{2}\hat{\lambda}_{u}^{4})V\approx
V^{\dagger}(\epsilon_{1}\hat{\lambda}_{u}^{2})V\equiv\epsilon_{1}\Delta_{q},$
(11)

which is equivalent to taking the approximation
$\lambda_{t}^{4}-\lambda_{u,c}^{4}\approx\lambda_{t}^{2}(\lambda_{t}^{2}-\lambda_{u,c}^{2})$
and redefining
$\epsilon_{1}+\epsilon_{2}\lambda_{t}^{2}\to\epsilon_{1}$.222Using the
unitarity of the CKM matrix, one can obtain
$[V^{\dagger}(\epsilon_{1}\hat{\lambda}_{u}^{2}+\epsilon_{2}\hat{\lambda}_{u}^{4})V]_{ij}=\lambda_{t}^{2}(\epsilon_{1}+\epsilon_{2}\lambda_{t}^{2})\delta_{ij}-(\lambda_{t}^{2}-\lambda_{u}^{2})[\epsilon_{1}+\epsilon_{2}(\lambda_{t}^{2}+\lambda_{u}^{2})]V_{ui}^{*}V_{uj}-(\lambda_{t}^{2}-\lambda_{c}^{2})[\epsilon_{1}+\epsilon_{2}(\lambda_{t}^{2}+\lambda_{c}^{2})]V_{ci}^{*}V_{cj}\approx\lambda_{t}^{2}(\epsilon_{1}+\epsilon_{2}\lambda_{t}^{2})\delta_{ij}-(\lambda_{t}^{2}-\lambda_{u}^{2})[\epsilon_{1}+\epsilon_{2}\lambda_{t}^{2}]V_{ui}^{*}V_{uj}-(\lambda_{t}^{2}-\lambda_{c}^{2})[\epsilon_{1}+\epsilon_{2}\lambda_{t}^{2}]V_{ci}^{*}V_{cj}=(\epsilon_{1}+\epsilon_{2}\lambda_{t}^{2})[V^{\dagger}\hat{\lambda}_{u}^{2}V]_{ij}$,
where the approximation
$\lambda_{t}^{2}+\lambda_{u,c}^{2}\approx\lambda_{t}^{2}$ have been used. For
the observables investigated in this paper, we have verified that the
numerical differences caused by this approximation are negligible. Finally, by
defining the basic FCNC couplings

$\displaystyle\Delta_{q}=V^{\dagger}\hat{\lambda}_{u}^{2}V,$ (12)

we arrive at

$\displaystyle\bar{d}_{L}U_{d}^{\dagger}\gamma_{\mu}{\mathcal{C}}^{\rm
MFV}U_{d}d_{L}\approx\bar{d}_{L}\gamma^{\mu}\big{(}\epsilon_{0}\mathds{1}+\epsilon_{1}\Delta_{q})d_{L}.$
(13)

Similarly, from the current $\bar{q}\,{\mathcal{C}}d$, the following scalar
interactions among the down-type quarks are obtained in the mass eigenbasis:

$\displaystyle\bar{d}_{L}U_{d}^{\dagger}{\mathcal{C}}^{\rm
MFV}W_{d}d_{R}\approx\bar{d}_{L}V^{\dagger}\big{(}\epsilon_{0}\mathds{1}+\epsilon_{1}\hat{\lambda}_{u}^{2}+\epsilon_{2}\hat{\lambda}_{u}^{4}\big{)}V\hat{\lambda}_{d}d_{R}\approx\bar{d}_{L}\big{(}\epsilon_{0}\hat{\lambda}_{d}+\epsilon_{1}\Delta_{q}\hat{\lambda}_{d}\big{)}d_{R},$
(14)

where $d_{R}$ is the right-handed down-type quark in the mass eigenbasis, and
the diagonal elements of the diagonal Yukawa coupling matrix for down-type
quarks are given by $\lambda_{d}^{i}=\sqrt{2}m_{d_{i}}/v$. Here the
approximation in eq. (11) has been used. In addition, due to the same flavour
structure as of the scalar interactions, the tensor interactions
$\bar{d}_{L}U_{d}^{\dagger}\sigma^{\mu\nu}{\mathcal{C}}^{\rm MFV}W_{d}d_{R}$
resulting from the current $\bar{q}\sigma^{\mu\nu}{\mathcal{C}}^{\rm MFV}d$
have the same MFV coupling as in eq. (14).

For the vector interactions originating from the current
$\bar{d}\gamma^{\mu}{\mathcal{C}}d$, the MFV coupling in the mass eigenbasis
is given by

$\displaystyle\bar{d}_{R}W_{d}^{\dagger}\gamma^{\mu}{\mathcal{C}}^{\rm
MFV}W_{d}d_{R}\approx\bar{d}_{R}\gamma^{\mu}(\epsilon_{0}\mathds{1})d_{R}.$
(15)

This implies that the vector interactions among the right-handed down-type
quarks are approximately flavour diagonal and universal. Thus, the FCNC vector
interactions among the right-handed down-type quarks are forbidden in the MFV
scenario.

As a summary, the MFV couplings for various down-type quark currents are given
in the mass eigenbasis by

$\displaystyle{\mathcal{C}}_{i}^{\rm
MFV}=\begin{cases}\epsilon_{0}^{i}\mathds{1}+\epsilon_{1}^{i}\Delta_{q}&\text{for
}\bar{d}_{L}\gamma^{\mu}{\mathcal{C}}_{i}d_{L},\\\
\epsilon_{0}^{i}\hat{\lambda}_{d}+\epsilon_{1}^{i}\Delta_{q}\hat{\lambda}_{d}&\text{for
}\bar{d}_{L}{\mathcal{C}}_{i}d_{R},\,\bar{d}_{L}\sigma^{\mu\nu}{\mathcal{C}}_{i}d_{R},\\\
\epsilon_{0}^{i}\mathds{1}&\text{for
}\bar{d}_{R}\gamma^{\mu}{\mathcal{C}}_{i}d_{R}.\end{cases}$ (16)

It is then clear that only two types of MFV couplings, $\Delta_{q}$ and
$\Delta_{q}\hat{\lambda}_{d}$, can generate FCNC interactions among the down-
type quarks and their numerical values are given, respectively, by

$\displaystyle\Delta_{q}=\begin{pmatrix}0.8&-3.3-1.5i&79.3+35.4i\\\
-3.3+1.5i&16.6&-397.5+8.1i\\\ 79.3-35.4i&-397.5-8.1i&9839.0\end{pmatrix}\times
10^{-4},$ (17)

and

$\displaystyle\Delta_{q}\hat{\lambda}_{d}=\begin{pmatrix}0.0021&-0.18-0.08i&191.3+85.4i\\\
-0.009+0.004i&0.88&-958.7+19.6i\\\ 0.21-0.10i&-21.1-0.4i&23728.1\\\
\end{pmatrix}\times 10^{-6}.$ (18)

One can see that a large hierarchy exists among the down-type FCNC
interactions, i.e., ${\mathcal{C}}_{bs}^{\rm MFV}\gg{\mathcal{C}}_{bd}^{\rm
MFV}\gg{\mathcal{C}}_{sd}^{\rm MFV}$. It is also seen that the flavour-
conserving couplings ${\mathcal{C}}_{bb}^{\rm MFV}$ is more than one order of
magnitude larger than all the other couplings.

###  2.2 Lepton sector

One can also apply the MFV hypothesis to the lepton sector. However, since the
mechanism of neutrino mass generation is still unknown, there are different
approaches to formulate the leptonic MFV Cirigliano:2005ck ; Cirigliano:2006su
; Davidson:2006bd ; Branco:2006hz ; Gavela:2009cd ; Alonso:2011jd ;
AristizabalSierra:2012myh . Here, we consider the realization of leptonic MFV
within the so-called minimal field content Cirigliano:2005ck ;
Cirigliano:2006su , in which the neutrino masses are generated by the Weinberg
operator. In this case, the Yukawa interactions in the lepton sector can be
written as

$\displaystyle-\Delta{\mathcal{L}}=\bar{e}\,Y_{e}H^{\dagger}l+\frac{1}{2\Lambda_{\rm
LN}}\big{(}\bar{l}^{c}\tau_{2}H\big{)}Y_{\nu}\big{(}H^{T}\tau_{2}l\big{)}+\mathrm{h.c.},$
(19)

where $l$ denotes the left-handed lepton doublet with the charge conjugated
field given by $l^{c}=-i\gamma_{2}l^{\ast}$, and $e$ is the right-handed
charged lepton singlet. $\Lambda_{\rm LN}$ denotes the breaking scale of the
lepton number symmetry $U(1)_{\rm LN}$. $Y_{e}$ and $Y_{\nu}$ stand for the
$3\times 3$ Yukawa coupling matrices in flavour space. In the absence of these
Yukawa couplings, the lepton sector respects the flavour symmetry

$\displaystyle G_{\rm LF}=SU(3)_{l}\otimes SU(3)_{e}.$ (20)

Analogous to the quark sector, this flavour symmetry can be restored by
promoting the Yukawa couplings $Y_{e,\nu}$ to spurion fields. Then, the
relevant building blocks are ${\mathsf{A}}_{\ell}=Y_{\nu}^{\dagger}Y_{\nu}$
and ${\mathsf{B}}_{\ell}=Y_{e}^{\dagger}Y_{e}$, which transform as
$(\textbf{1}\oplus\textbf{8},\textbf{1})$ under the lepton flavour group
$G_{\rm LF}$.

For the three SMEFT operators, ${\mathcal{Q}}_{lq}^{(1,3)}$ and
${\mathcal{Q}}_{ld}$, relevant to the $B^{+}\to K^{+}\nu\bar{\nu}$ decay, they
contain the lepton current

$\displaystyle\bar{l}\gamma^{\mu}{\mathcal{C}}l,$ (21)

where ${\mathcal{C}}$ denotes the Wilson coefficient of the operator
${\mathcal{Q}}$ and is a $3\times 3$ matrix in lepton flavour space. In the
MFV hypothesis, analogous to the quark current, the Wilson coefficient
${\mathcal{C}}$ should take the form ${\mathcal{C}}_{\rm
MFV}=h({\mathsf{A}}_{\ell},{\mathsf{B}}_{\ell})$, where
$h({\mathsf{A}}_{\ell},{\mathsf{B}}_{\ell})$ is a finite polynomial of
${\mathsf{A}}_{\ell}$ and ${\mathsf{B}}_{\ell}$. After neglecting all the
terms involving ${\mathsf{B}}_{\ell}$, which are suppressed by the small
lepton Yukawa couplings $Y_{e}$, we obtain

$\displaystyle{\mathcal{C}}_{\rm
MFV}\approx\kappa_{0}+\kappa_{1}{\mathsf{A}}_{\ell}+\kappa_{2}{\mathsf{A}}_{\ell}^{2},$
(22)

where the coefficients $\kappa_{0,1,2}$ are free real parameters. In the
numerical analysis, we keep only the leading lepton flavour violation term
${\mathsf{A}}_{\ell}$ for simplicity, i.e., $\kappa_{2}=0$. Turning to the
lepton mass eigenbasis, the current $\bar{l}\gamma^{\mu}{\mathcal{C}}l$ gives
in the MFV hypothesis the following interactions:

$\displaystyle\bar{e}_{L}\gamma^{\mu}(\kappa_{0}\mathds{1}+\kappa_{0}\Delta_{\ell})e_{L}+\bar{\nu}_{L}\gamma^{\mu}(\kappa_{0}\mathds{1}+\kappa_{0}\hat{\lambda}_{\nu}^{2})\nu_{L},$
(23)

where the basic LFV coupling $\Delta_{\ell}$ can be obtained from
${\mathsf{A}}_{\ell}$ and takes the form

$\displaystyle\Delta_{\ell}=U\hat{\lambda}_{\nu}^{2}U^{\dagger},$ (24)

where $U$ is the PMNS matrix. Here $e_{L}$ and $\nu_{L}$ denote the left-
handed charged lepton and neutrino in the mass eigenbasis, respectively.
$\hat{\lambda}_{\nu}\equiv 2m_{\nu}\Lambda_{\rm LN}/v^{2}$ stands for the
diagonal effective neutrino Yukawa coupling and $m_{\nu}$ the diagonal
neutrino mass matrix. As can be seen from eq. (23), the lepton vector current
does not induce neutral LFV in the leptonic MFV with minimal field
content.333This is, however, not true for the most generic leptonic MFV. For
example, in the leptonic MFV with type-I seesaw mechanism, the basic LFV
coupling $\Delta_{\ell}=U\hat{\lambda}_{\nu}U^{\dagger}$ should be replaced by
$U\hat{\lambda}_{\nu}^{1/2}OO^{\dagger}\hat{\lambda}_{\nu}^{1/2}U^{\dagger}$,
where $O$ is a general complex orthogonal matrix satisfying
$OO^{T}=\mathds{1}$ Chiang:2017hlj . Therefor, any non-diagonal matrix
$OO^{\dagger}$ can induce neutral LFV.

As discussed above, the charged LFV interactions are governed by the MFV
coupling $\Delta_{\ell}$. Numerically, we obtain

$\displaystyle\Delta_{\ell}^{\rm
NO}=\begin{pmatrix}-0.19-0.01i\;&-0.25-0.02i\;&\hphantom{-}0.31-0.04i\\\
\hphantom{-}0.12+0.01i\;&\hphantom{-}0.28-0.00i\;&\hphantom{-}0.29+0.04i\\\
-0.37-0.01i\;&\hphantom{-}0.21-0.05i\;&-0.03+0.01i\end{pmatrix},$ (25)

for normal ordering (NO), and

$\displaystyle\Delta_{\ell}^{\rm
IO}=\begin{pmatrix}0.21+0.09i\;&-0.34+0.05i\;&\hphantom{-}0.03+0.11i\\\
0.31+0.12i\;&\hphantom{-}0.19+0.00i\;&-0.15-0.14i\\\
0.12-0.02i\;&\hphantom{-}0.04-0.19i\;&\hphantom{-}0.34-0.10i\end{pmatrix},$
(26)

for inverted ordering (IO) of neutrino masses. Here $\Lambda_{\rm
LN}=10^{14}{\,\rm GeV}$ and the lightest neutrino mass $m_{\nu_{1}}=0.2\,{\rm
eV}$ have been taken. One can see that the magnitudes of all the matrix
elements of $\Delta_{\ell}$ are of similar size except for
$(\Delta_{\ell}^{\rm NO})_{33}$, which has much smaller magnitude than the
other couplings.

Let us finally make a comment on the DSMEFT with MFV. In general, similar to
the SM fermions, the DM particles can also appear in a form of multiple
generations. In this case, the SM flavour symmetry could be extended to
include the global symmetry associated with the DM multiplet and the MFV
hypothesis can be generalized to the DM sector Batell:2011tc ; Lopez-
Honorez:2013wla ; Agrawal:2014aoa ; Chen:2015jkt ; Acaroglu:2021qae . In this
paper, we focus on a minimal DM sector, which contains only one or two DM
generations, and leave a detained study of the DM with MFV for future work.

##  3 $\boldsymbol{M_{1}\to M_{2}+\nu\bar{\nu}}$

In this section, we first summarize the SMEFT framework with the most general
flavour structure as well as in the MFV hypothesis. Then, the NP effects on
various rare FCNC decays are discussed.

###  3.1 SMEFT

If new particles are much heavier than the electroweak scale and the
electroweak symmetry breaking is realized linearly, their effects can be
parameterized by a series of higher dimensional gauge-invariant operators in
the SMEFT (see refs. Brivio:2017vri ; Isidori:2023pyp  for recent reviews).
The SMEFT Lagrangian at dim-6 takes the form Buchmuller:1985jz ;
Grzadkowski:2010es

$\displaystyle\mathcal{L}_{\rm
SMEFT}=\sum_{i}\frac{{\mathcal{C}}_{i}}{\Lambda^{2}}{\mathcal{Q}}_{i},$ (27)

where $\Lambda$ denotes the NP scale. The dim-6 effective operators
${\mathcal{Q}}_{i}$ are built out of the SM fields and respect the SM gauge
symmetry $SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}$, with
${\mathcal{C}}_{i}$ the corresponding dimensionless Wilson coefficients.

For the $b\to s\nu\bar{\nu}$ transition, the relevant operators are given by

$\displaystyle{\mathcal{Q}}_{lq}^{(1)}$
$\displaystyle=\big{(}\bar{l}_{p}\gamma^{\mu}l_{r}\big{)}\big{(}\bar{q}_{s}\gamma_{\mu}q_{t}\big{)}\,,$
$\displaystyle{\mathcal{Q}}_{Hq}^{(1)}$
$\displaystyle=\big{(}H^{\dagger}i\overleftrightarrow{D}_{\mu}H\big{)}\big{(}\bar{q}_{p}\gamma^{\mu}q_{r}\big{)}\,,$
(28) $\displaystyle{\mathcal{Q}}_{lq}^{(3)}$
$\displaystyle=\big{(}\bar{l}_{p}\gamma^{\mu}\tau^{I}l_{r}\big{)}\big{(}\bar{q}_{s}\tau^{I}\gamma_{\mu}q_{t}\big{)}\,,$
$\displaystyle{\mathcal{Q}}_{Hq}^{(3)}$
$\displaystyle=\big{(}H^{\dagger}i\overleftrightarrow{D}^{I}_{\mu}H\big{)}\big{(}\bar{q}_{p}\tau^{I}\gamma^{\mu}q_{r}\big{)},$
$\displaystyle{\mathcal{Q}}_{ld}$
$\displaystyle=\big{(}\bar{l}_{p}\gamma^{\mu}l_{r}\big{)}\big{(}\bar{d}_{s}\gamma_{\mu}d_{t}\big{)},$
$\displaystyle{\mathcal{Q}}_{Hd}$
$\displaystyle=\big{(}H^{\dagger}i\overleftrightarrow{D}_{\mu}H\big{)}\big{(}\bar{d}_{p}\gamma^{\mu}d_{r}\big{)},$

where
$D_{\mu}=\partial_{\mu}+ig_{s}T^{a}G^{a}_{\mu}+ig_{2}t^{I}W_{\mu}^{I}+ig_{1}YB_{\mu}$
denotes the covariant derivative under the SM gauge symmetry. Here
$T^{a}=\lambda^{a}/2$ and $t^{I}=\tau^{I}/2$ are the $SU(3)$ and $SU(2)$
generators, where $\lambda^{a}$ and $\tau^{I}$ stand for the Gell-Mann and
Pauli matrices, respectively. $l$ and $q$ denote the left-handed lepton and
quark doublets respectively, while $d$ is the right-handed down-type quark
singlet, all being given in the mass eigenbasis. The generation indices of the
SM fields are characterized by $p$, $r$, $s$, $t$.

The operators ${\mathcal{Q}}_{Hq}^{(1)}$, ${\mathcal{Q}}_{Hq}^{(3)}$ and
${\mathcal{Q}}_{Hd}$ can induce a tree-level $\bar{s}\gamma^{\mu}bZ_{\mu}$
coupling, and thus contribute to the $b\to s\nu\bar{\nu}$ transition. However,
the same $\bar{s}\gamma^{\mu}bZ_{\mu}$ interaction also affects the $b\to
s\ell^{+}\ell^{-}$ processes Altmannshofer:2009ma ; Buras:2014fpa . Due to the
$SU(2)_{L}$ gauge invariance, its contributions to these two processes are of
equal magnitudes and lepton flavour universal. Thus, it is impossible to
explain the Belle II data on the branching ratio ${\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})$ while satisfying the stringent constraints from the
precisely measured $b\to s\ell^{+}\ell^{-}$ observables, e.g., the branching
ratio ${\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$. Explicitly, it is found that
${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ can only be enhanced by about
$20\%$ after considering other relevant experimental constraints
Allwicher:2023xba . Therefore, we will not consider these operators and focus
on the remaining three four-fermion operators in the following.

In the MFV hypothesis, by using eqs. (13), (15) and (23), the relevant four-
fermion operators take the following forms:

$\displaystyle{\mathcal{C}}_{ld}{\mathcal{Q}}_{ld}$
$\displaystyle=\Big{[}\bar{e}_{L}\gamma_{\mu}\left(\kappa_{0}\mathds{1}+\kappa_{1}\Delta_{\ell}\right)e_{L}+\bar{\nu}_{L}\gamma_{\mu}\big{(}\kappa_{0}\mathds{1}+\kappa_{1}\hat{\lambda}_{\nu}^{2}\big{)}\nu_{L}\Big{]}\Big{[}\bar{d}_{R}\gamma^{\mu}(\epsilon_{0}\mathds{1})d_{R}\Big{]},$
(29) $\displaystyle{\mathcal{C}}_{lq}^{(1,3)}{\mathcal{Q}}_{lq}^{(1,3)}$
$\displaystyle\supset\Big{[}\bar{e}_{L}\gamma_{\mu}\left(\kappa_{0}\mathds{1}+\kappa_{1}\Delta_{\ell}\right)e_{L}\pm\bar{\nu}_{L}\gamma_{\mu}\big{(}\kappa_{0}\mathds{1}+\kappa_{1}\hat{\lambda}_{\nu}^{2}\big{)}\nu_{L}\Big{]}\Big{[}\bar{d}_{L}\gamma^{\mu}\left(\epsilon_{0}\mathds{1}+\epsilon_{1}\Delta_{q}\right)d_{L}\Big{]},$

in the mass eigenbasis. Here we have suppressed the operator indices for the
MFV parameters $\epsilon_{0,1}$ and $\kappa_{0,1}$ are suppressed. It is noted
that ${\mathcal{Q}}_{ld}$ does not involve any quark FCNC interactions, and
thus has no contribution to the $b\to s\nu\bar{\nu}$ transition.

###  3.2 Observables

####  3.2.1 $b\to s(d)\nu\bar{\nu}$ decay

The $B\to K^{(*)}\nu\bar{\nu}$ decays are induced by the $b\to
s\nu_{i}\bar{\nu}_{j}$ transitions. In the most general case, the effective
weak Hamiltonian contains only two operators and can be written as
Altmannshofer:2009ma ; Buras:2014fpa

$\displaystyle\mathcal{H}_{\rm eff}^{b\to
s\nu\nu}=-\dfrac{4G_{F}}{\sqrt{2}}\frac{\alpha_{e}}{4\pi}V_{tb}V_{ts}^{*}\sum_{i,j=e,\mu,\tau}\left(C^{\nu_{i}\nu_{j}}_{L}{\mathcal{O}}_{L}^{\nu_{i}\nu_{j}}\,+C^{\nu_{i}\nu_{j}}_{R}{\mathcal{O}}_{R}^{\nu_{i}\nu_{j}}\right)+\mathrm{h.c.},$
(30)

with the effective operators defined by

$\displaystyle{\mathcal{O}}_{L}^{\nu_{i}\nu_{j}}=\big{(}\bar{s}\gamma_{\mu}P_{L}b\big{)}\big{(}\bar{\nu}_{i}\gamma^{\mu}P_{L}\nu_{j}\big{)},\qquad{\mathcal{O}}_{R}^{\nu_{i}\nu_{j}}=\big{(}\bar{s}\gamma_{\mu}P_{R}b\big{)}\big{(}\bar{\nu}_{i}\gamma^{\mu}P_{L}\nu_{j}\big{)}.$
(31)

In the SM, $C_{R}^{\nu_{i}\nu_{j}}$ is negligible and
$C_{L}^{\nu_{i}\nu_{j}}=C_{L}^{\rm SM}\delta_{ij}$ with $C_{L}^{\rm
SM}=-12.64\pm 0.14$ Buras:2014fpa  after including the NLO QCD corrections
Buchalla:1993bv ; Misiak:1999yg ; Buchalla:1998ba  and the two-loop
electroweak contributions Brod:2010hi . In the SMEFT, by comparing eqs. (27)
and (30), one can derive that the NP contributions to these two operators are
given, respectively, by

$\displaystyle C_{L,{\rm
NP}}^{\nu_{i}\nu_{j}}=c_{\nu}\big{(}\big{[}\mathcal{C}_{lq}^{(1)}\big{]}_{ij23}-\big{[}\mathcal{C}_{lq}^{(3)}\big{]}_{ij23}\big{)},\qquad
C_{R,{\rm NP}}^{\nu_{i}\nu_{j}}=c_{\nu}\big{[}\mathcal{C}_{ld}\big{]}_{ij23},$
(32)

with the normalization factor
$c_{\nu}\equiv(2\pi/(\alpha_{e}V_{tb}V_{ts}^{*}))\cdot(v^{2}/\Lambda^{2})$. In
the MFV framework, by using eq. (29), we find that

$\displaystyle C_{L,{\rm MFV}}^{\nu_{i}\nu_{j}}$
$\displaystyle=c_{\nu}\Delta_{q}^{23}\big{[}\epsilon_{1}^{(1)}\big{(}\kappa_{0}^{(1)}+\kappa_{1}^{(1)}\lambda_{\nu_{i}}^{2}\big{)}-\big{(}(1)\to(3)\big{)}\big{]}\delta_{ij},$
(33) $\displaystyle C_{R,{\rm MFV}}^{\nu_{i}\nu_{j}}$ $\displaystyle=0,$

where $\epsilon_{1}^{(n)}$ and $\kappa_{0,1}^{(n)}$ denote the MFV parameters
associated with the operator ${\mathcal{Q}}_{lq}^{(n)}$, with $n=1,3$. Thus,
one can see that the transition $b\to s\nu\bar{\nu}$ in the MFV hypothesis is
only induced by the effective operator ${\mathcal{O}}_{L}^{\nu_{i}\nu_{i}}$,
which also conserves the lepton flavour. In addition, it is straightforward to
obtain the operators and Wilson coefficients for the $b\to
d\nu_{i}\bar{\nu}_{j}$ transitions by changing the corresponding flavour
indices.

Starting with the effective weak Hamiltonian in eq. (30), we can write the
dineutrino invariant mass spectrum of the $B\to K^{*}\nu\bar{\nu}$ decay as
Altmannshofer:2009ma ; Buras:2014fpa

$\displaystyle\frac{d\Gamma(B\to
K^{*}\nu\bar{\nu})}{ds_{B}}=m_{B}^{2}\sum_{i,j=e,\mu,\tau}\Big{(}\left|A^{ij}_{\perp}(s_{B})\right|^{2}+\left|A^{ij}_{\|}(s_{B})\right|^{2}+\left|A^{ij}_{0}(s_{B})\right|^{2}\Big{)},$
(34)

where $s_{B}=q^{2}/m_{B}^{2}$, and $q^{2}$ denotes the dineutrino invariant
mass squared in the $B$-meson rest frame. The three polarization amplitudes
$A_{\perp,\|,0}^{ij}$ are functions of the $B\to K^{*}$ transition form
factors and the Wilson coefficients, and take the form

$\displaystyle A^{ij}_{\perp}(s_{B})$
$\displaystyle=2\sqrt{2}N\,\lambda^{1/2}(1,\tilde{m}_{K^{*}}^{2},s_{B})\,(C^{\nu_{i}\nu_{j}}_{L}+C^{\nu_{i}\nu_{j}}_{R})\frac{V(s_{B})}{1+\tilde{m}_{K^{*}}},$
$\displaystyle A^{ij}_{\|}(s_{B})$
$\displaystyle=-2\sqrt{2}N\,(1+\tilde{m}_{K^{*}})\,(C^{\nu_{i}\nu_{j}}_{L}-C^{\nu_{i}\nu_{j}}_{R})\,A_{1}(s_{B}),$
$\displaystyle A^{ij}_{0}(s_{B})$
$\displaystyle=-\frac{N(C^{\nu_{i}\nu_{j}}_{L}-C^{\nu_{i}\nu_{j}}_{R})}{\widetilde{m}_{K^{*}}\sqrt{s_{B}}}\biggl{[}(1-\tilde{m}_{K^{*}}^{2}-s_{B})(1+\tilde{m}_{K^{*}})A_{1}(s_{B})$
$\displaystyle\hskip
184.9429pt-\lambda(1,\tilde{m}_{K^{*}}^{2},s_{B})\frac{A_{2}(s_{B})}{1+\tilde{m}_{K^{*}}}\biggr{]},$
(35)

with $\tilde{m}_{K^{*}}=m_{K^{*}}/m_{B}$ and
$N=1/2\,V_{tb}V_{ts}^{*}\left[G_{F}^{2}\alpha_{e}^{2}m_{B}^{3}s_{B}\lambda^{1/2}(1,\tilde{m}_{K^{*}}^{2},s_{B})/(3\cdot
2^{10}\pi^{5})\right]^{1/2}$. Similarly, the dineutrino invariant mass
distribution of the $B\to K\nu\bar{\nu}$ decay can be written as
Colangelo:1996ay ; Altmannshofer:2009ma ; Buras:2014fpa

$\displaystyle\frac{d\Gamma(B\to K\nu\bar{\nu})}{ds_{B}}$
$\displaystyle=\frac{G_{F}^{2}\alpha_{e}^{2}}{3\cdot
2^{10}\pi^{5}}|V_{tb}V_{ts}^{*}|^{2}m_{B}^{5}\lambda^{3/2}\left(1,\tilde{m}_{K}^{2},s_{B}\right)\left[f_{+}^{B\to
K}\left(s_{B}\right)\right]^{2}$ $\displaystyle\hskip
128.0374pt\times\sum_{i,j=e,\mu,\tau}\big{\lvert}C^{\nu_{i}\nu_{j}}_{L}+C^{\nu_{i}\nu_{j}}_{R}\big{\rvert}^{2}\,.$
(36)

For convenience, all the relevant form factors present in the above
expressions are given in appendix A.2. In addition, it is straightforward to
obtain the distributions of other $b\to s(d)\nu\bar{\nu}$ decays, such as
$B\to\pi\nu\bar{\nu}$ and $B_{s}\to\phi\nu\bar{\nu}$, from the above
expressions with some replacements. Finally, one should notice that for
charged $B$-meson decays, there is a long-distance contribution arising from
the tree-level weak annihilation mediated by the on-shell $\tau$ lepton
Kamenik:2009kc . For example, the $B^{+}\to\tau^{+}(\to K^{+}\bar{\nu})\nu$
process provides an additional ${\mathcal{O}}(10\%)$ contribution to the
$B^{+}\to K^{+}\nu\bar{\nu}$ decay, and constitutes a background in the
experimental extraction of the branching ratio Belle-II:2023esi . Here we do
not include the tree-level contribution, which has already been subtracted
away during the experimental data analysis Belle-II:2023esi .

####  3.2.2 $b\to s\ell\ell$ decay

The quark-level $b\to s\ell^{+}\ell^{-}$ transitions, such as the rare
$B_{s}\to\ell^{+}\ell^{-}$, $B\to K^{(*)}\ell^{+}\ell^{-}$ and
$B_{s}\to\phi\ell^{+}\ell^{-}$ decays, can provide promising probes of various
NP effects (see refs. London:2021lfn ; Capdevila:2023yhq  for a recent
review). For a general $b\to s\ell_{i}^{-}\ell_{j}^{+}$ transition, the
effective weak Hamiltonian can be written as Buchalla:1995vs

$\displaystyle\mathcal{H}_{\rm eff}^{b\to
s\ell\ell}=-\frac{4G_{F}}{\sqrt{2}}\frac{e^{2}}{16\pi^{2}}V_{tb}V_{ts}^{*}\sum_{k}C_{k}{\mathcal{O}}_{k}+\mathrm{h.c.},$
(37)

where the most relevant operators are ${\mathcal{O}}_{9}^{(\prime)}$ and
${\mathcal{O}}_{10}^{(\prime)}$. They are defined, respectively, by

$\displaystyle{\mathcal{O}}_{9}^{ij}$
$\displaystyle=\big{(}\bar{s}\gamma_{\mu}P_{L}b\big{)}\big{(}\bar{\ell}_{i}\gamma^{\mu}\ell_{j}\big{)},$
$\displaystyle{\mathcal{O}}_{10}^{ij}$
$\displaystyle=\big{(}\bar{s}\gamma_{\mu}P_{L}b\big{)}\big{(}\bar{\ell}_{i}\gamma^{\mu}\gamma_{5}\ell_{j}\big{)},$
(38) $\displaystyle{\mathcal{O}}_{9}^{\prime\,ij}$
$\displaystyle=\big{(}\bar{s}\gamma_{\mu}P_{R}b\big{)}\big{(}\bar{\ell}_{i}\gamma^{\mu}\ell_{j}\big{)},$
$\displaystyle{\mathcal{O}}_{10}^{\prime\,ij}$
$\displaystyle=\big{(}\bar{s}\gamma_{\mu}P_{R}b\big{)}\big{(}\bar{\ell}_{i}\gamma^{\mu}\gamma_{5}\ell_{j}\big{)}.$

In the SMEFT, the NP contributions to their Wilson coefficients are given by

$\displaystyle C_{9,{\rm NP}}^{ij}$
$\displaystyle=+c_{\ell}\big{(}\big{[}{\mathcal{C}}_{lq}^{(1)}\big{]}_{ij23}+\big{[}{\mathcal{C}}_{lq}^{(3)}\big{]}_{ij23}\big{)},$
$\displaystyle C_{9,{\rm NP}}^{\prime\,ij}$
$\displaystyle=+c_{\ell}\big{[}{\mathcal{C}}_{ld}\big{]}_{ij23},$ (39)
$\displaystyle C_{10,{\rm NP}}^{ij}$
$\displaystyle=-c_{\ell}\big{(}\big{[}{\mathcal{C}}_{lq}^{(1)}\big{]}_{ij23}+\big{[}{\mathcal{C}}_{lq}^{(3)}\big{]}_{ij23}\big{)},$
$\displaystyle C_{10,{\rm NP}}^{\prime\,ij}$
$\displaystyle=-c_{\ell}\big{[}{\mathcal{C}}_{ld}\big{]}_{ij23},$

with the normalization constant
$c_{\ell}\equiv(\pi/(\alpha_{e}V_{tb}V_{ts}^{*}))\cdot(v^{2}/\Lambda^{2})$. In
the MFV framework, from eq. (29), we can further obtain

$\displaystyle C_{9,{\rm MFV}}^{ij}$ $\displaystyle=-C_{10,{\rm
MFV}}^{ij}=c_{\ell}\Delta_{q}^{23}\big{[}\epsilon_{1}^{(1)}\big{(}\kappa_{0}^{(1)}\delta^{ij}+\kappa_{1}^{(1)}\Delta_{\ell}^{ij}\big{)}+\big{(}(1)\to(3)\big{)}\big{]},$
$\displaystyle C_{9,{\rm MFV}}^{\prime\,ij}$ $\displaystyle=C_{10,{\rm
MFV}}^{\prime\,ij}=0,$ (40)

where $\epsilon_{1}^{(n)}$ and $\kappa_{0,1}^{(n)}$ denote the MFV parameters
for the operator ${\mathcal{Q}}_{lq}^{(n)}$ ($n=1,3$). For the
$B\to\pi\ell_{i}\ell_{j}$, $B\to K^{(*)}\ell_{i}\ell_{j}$ and
$B_{s}\to\ell_{i}\ell_{j}$ decays, which are all induced by the quark-level
$b\to s(d)\ell_{i}\ell_{j}$ transitions and will be included in our numerical
analysis, the explicit expressions of their branching ratios can be found in,
e.g., refs. Crivellin:2015era ; Lee:2015qra ; Zhang:2018nmy .

####  3.2.3 $s\to d\nu\bar{\nu}$ decay

Both the $K^{+}\to\pi^{+}\nu\bar{\nu}$ and $K_{L}\to\pi^{0}\nu\bar{\nu}$
decays are induced by the $s\to d\nu_{i}\bar{\nu_{j}}$ transition. The
corresponding effective weak Hamiltonian, including the effective operators
and their Wilson coefficients, can be obtained from the ones of the $b\to
s\nu\bar{\nu}$ transition in section 3.2.1 by replacing the related flavour
indices. In addition, for the $K^{+}\to\pi^{+}\nu\bar{\nu}$ decay, a non-
negligible charm contribution must be taken into account Buchalla:1998ba ;
Misiak:1999yg .

In absence of the operator $O_{R}^{\nu_{i}\nu_{j}}$, the branching ratio of
the $K^{+}\to\pi^{+}\nu\bar{\nu}$ decay is given by Buras:2004uu ;
Mescia:2007kn ; Cirigliano:2011ny

$\displaystyle{\mathcal{B}}(K^{+}\to\pi^{+}\nu\bar{\nu})$
$\displaystyle=\kappa_{+}(1+\Delta_{\rm
EM})\sum_{i,j=e,\mu,\tau}\frac{1}{3}\bigg{[}\bigg{(}\sin^{2}\theta_{W}\frac{{\rm
Im}[\lambda_{t}C_{L}^{\nu_{i}\nu_{j}}]}{2\lambda^{5}}\bigg{)}^{2}$
$\displaystyle\hskip 56.9055pt+\bigg{(}\frac{{\rm
Re}\lambda_{c}}{\lambda}P_{c}(X)\delta^{ij}-\sin^{2}\theta_{W}\frac{{\rm
Re}[\lambda_{t}C_{L}^{\nu_{i}\nu_{j}}]}{2\lambda^{5}}\bigg{)}^{2}\,\bigg{]}\,,$
(41)

where $\lambda_{i}=V_{is}^{*}V_{id}$, $\lambda=|V_{us}|$, and
$\sin^{2}\theta_{W}=0.23116$ is the weak mixing angle. The parameter
$\kappa_{+}=(5.173\pm 0.025)\times 10^{-11}(\lambda/0.225)^{8}$ contains,
besides some other factors, the hadronic matrix element of the weak current
$\bar{s}\gamma_{\mu}P_{L}d$ that can be extracted from the $K_{\ell 3}$ data
with the help of isospin symmetry Buras:2015qea . $\Delta_{\rm EM}=-0.003$
accounts for the isospin-breaking electromagnetic corrections Mescia:2007kn ;
Cirigliano:2011ny . $P_{c}(X)=0.404\pm 0.024$ comprises the short- and long-
distance charm contributions Buras:2015qea ; Buras:2006gb ; Brod:2008ss ;
Isidori:2005xm .

The rare decay $K_{L}\to\pi^{0}\nu\bar{\nu}$ proceeds almost entirely through
direct CP violation in the SM Littenberg:1989ix ; Buchalla:1998ux , and the
charm contribution can be fully neglected Buras:2004uu . The branching ratio
of this decay can be written as Buras:2004uu ; Buchalla:1996fp

$\displaystyle{\mathcal{B}}(K_{L}\to\pi^{0}\nu\bar{\nu})=\kappa_{L}(1-\delta_{\epsilon})\sum_{i,j=e,\mu,\tau}\frac{1}{3}\left(\frac{\sin^{2}\theta_{W}}{2\lambda^{5}}\text{Im}\big{[}\lambda_{t}C_{L}^{\nu_{i}\nu_{j}}\big{]}\right)^{2},$
(42)

where $\kappa_{L}=(2.231\pm 0.013)\times 10^{-10}(\lambda/0.225)^{8}$ encodes,
besides some other factors, the hadronic matrix element that can be again
extracted from the $K_{\ell 3}$ decay with the help of isospin symmetry
Marciano:1996wy ; Mescia:2007kn . The parameter $\delta_{\epsilon}$ encodes,
on the other hand, the highly suppressed indirect CP-violating contribution
Buchalla:1998ux .

##  4 $\boldsymbol{M_{1}\to M_{2}+\text{DM}}$

In this section, we first recapitulate the theoretical framework of the DSMEFT
and DLEFT, and then detail the calculation of $M_{1}\to M_{2}+\text{DM}$
decays.

###  4.1 DSMEFT

It is known that EFT can provide a systematic framework for parameterizing
various effects induced by the DM particles. The DSMEFT is a generalization of
the SMEFT Buchmuller:1985jz ; Grzadkowski:2010es  to include the DM particles,
which can be either lighter or heavier than the electroweak scale and are
gauge singlets under the SM gauge group Aebischer:2022wnl ; Liang:2023yta ;
Song:2023jqm . It is also assumed that electroweak symmetry breaking must be
implemented via the usual Higgs mechanism. Following the convention of ref.
Aebischer:2022wnl , the DSMEFT Lagrangian takes the form

$\displaystyle{\mathcal{L}}_{\mathrm{DSMEFT}}$
$\displaystyle\supset\frac{1}{\Lambda}\sum_{i}{\mathcal{C}}_{i}{\mathcal{Q}}_{i}^{(5)}+\frac{1}{\Lambda^{2}}\sum_{j}{\mathcal{C}}_{j}{\mathcal{Q}}_{j}^{(6)},$
(43)

where $\Lambda$ is the NP scale above which new degrees of freedom appear and
interact with the SM and DM particles. The dim-$d$ effective operators
${\mathcal{Q}}_{i}^{(d)}$ involve at least one DM field and respect the SM
$SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}$ gauge symmetry.
${\mathcal{C}}_{i}$ denote the corresponding dimensionless Wilson
coefficients. We will define the DSMEFT at the electroweak scale $\mu_{\rm
EW}$.

In the DSMEFT framework, three types of DM particles can be added: spin-0
scalar $\phi$, spin-1/2 fermion $\chi$ and spin-1 vector $X$. Generally, each
type of the DM particles can have multiple generations. For simplicity, the
case of one generation DM will be assumed in this paper. For the operators
vanishing in the case of one generation, however, two generations of DM
particles will be considered. Furthermore, in order to affect the rare
$B$-meson and kaon FCNC decays, only the operators that can induce tree-level
FCNC interactions among the down-type quarks are included in our analysis.
Finally, motivated by the large excess of ${\mathcal{B}}(B^{+}\to
K^{+}+\nu\bar{\nu})$ reported by Belle II Belle-II:2023esi , we will focus on
the dim-5 and dim-6 operators, which are less suppressed than the ones with
dim $\geq 7$. Consequently, the DSMEFT operators we are considering must
involve two down-type quarks and at least one DM field. Furthermore, in order
to mimic the rare $B$-meson and kaon FCNC decays involving a dineutrino final
state, the DM particles must be light enough to be produced in these decays.
For each type of the DM particles, the relevant effective operators are given
in the following.

For the scalar DM particles, there are four operators,

$\displaystyle{\mathcal{Q}}_{d\phi}$
$\displaystyle=\big{(}\bar{q}_{p}d_{r}H\big{)}\phi+\mathrm{h.c.},$
$\displaystyle{\mathcal{Q}}_{d\phi^{2}}$
$\displaystyle=\big{(}\bar{q}_{p}d_{r}H\big{)}\phi^{2}+\mathrm{h.c.},$
$\displaystyle{\mathcal{Q}}_{\phi q}$
$\displaystyle=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}\big{(}i\phi_{1}\overleftrightarrow{\partial^{\mu}}\phi_{2}\big{)},$
$\displaystyle{\mathcal{Q}}_{\phi d}$
$\displaystyle=\big{(}\bar{d}_{p}\gamma_{\mu}d_{r}\big{)}\big{(}i\phi_{1}\overleftrightarrow{\partial^{\mu}}\phi_{2}\big{)},$
(44)

where
$\phi_{1}\overleftrightarrow{\partial^{\mu}}\phi_{2}=\phi_{1}({\partial^{\mu}}\phi_{2})-({\partial^{\mu}}\phi_{1})\phi_{2}$,
and $\phi_{(i)}$ is the real scalar DM field with the generation index
characterized by $i$. For the operators which are not hermitian,
“$+\mathrm{h.c.}$” should be added to them to indicate that the hermitian
conjugated operators should be included in the DSMEFT Lagrangian eq. (43). In
addition, constructions of the operators ${\mathcal{Q}}_{\phi q}$ and
${\mathcal{Q}}_{\phi d}$ requires at least two generations of scalar DM
fields.

For the fermionic DM particles, there are only two operators,

$\displaystyle{\mathcal{Q}}_{q\chi}=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}(\bar{\chi}\gamma^{\mu}\chi\big{)},\qquad\qquad{\mathcal{Q}}_{d\chi}=\big{(}\bar{d}_{p}\gamma_{\mu}d_{r}\big{)}(\bar{\chi}\gamma^{\mu}\chi\big{)},$
(45)

where $\chi$ denotes the fermionic DM field. In this paper, all the
observables considered are insensitive to the chirality of the fermionic DM.
Without loss of generality, $\chi$ will be therefore chosen as the right-
handed Dirac field.

For the vector DM particles, by using the field strength tensor as the
building block, there is only one operator,

$\displaystyle{\mathcal{Q}}_{dHX}=\big{(}\bar{q}_{p}\sigma_{\mu\nu}d_{r}\big{)}HX^{\mu\nu}+\mathrm{h.c.},$
(46)

where $X_{\mu\nu}=\partial_{\mu}X_{\nu}-\partial_{\nu}X_{\mu}$, and $X_{\mu}$
denotes the real vector DM field. Taking $X_{\mu}$ an additional building
block, the following operators also appear Song:2023jqm :

$\displaystyle{\mathcal{Q}}_{dX}$
$\displaystyle=\big{(}\bar{d}_{p}\gamma_{\mu}d_{r}\big{)}X^{\mu},$
$\displaystyle{\mathcal{Q}}_{HdX}$
$\displaystyle=\big{(}H^{\dagger}H\big{)}\big{(}\bar{d}_{p}\gamma^{\mu}d_{r}\big{)}X_{\mu},$
$\displaystyle{\mathcal{Q}}_{qX}$
$\displaystyle=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}X^{\mu},$
$\displaystyle{\mathcal{Q}}_{HqX}^{(1)}$
$\displaystyle=\big{(}H^{\dagger}H\big{)}\big{(}\bar{q}_{p}\gamma^{\mu}q_{r}\big{)}X_{\mu},$
$\displaystyle{\mathcal{Q}}_{dX^{2}}$
$\displaystyle=\big{(}\bar{q}_{p}d_{r}H\big{)}X_{\mu}X^{\mu}+\mathrm{h.c.},$
$\displaystyle{\mathcal{Q}}_{HqX}^{(3)}$
$\displaystyle=\big{(}H^{\dagger}\tau^{I}H\big{)}\big{(}\bar{q}_{p}\tau^{I}\gamma^{\mu}q_{r}\big{)}X_{\mu},$
$\displaystyle{\mathcal{Q}}_{qXX}$
$\displaystyle=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}X^{\mu\nu}X_{\nu}\,,$
$\displaystyle{\mathcal{Q}}_{dXX}$
$\displaystyle=\big{(}\bar{d}_{p}\gamma_{\mu}d_{r}\big{)}X^{\mu\nu}X_{\nu},$
$\displaystyle{\mathcal{Q}}_{q\tilde{X}X}$
$\displaystyle=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}\tilde{X}^{\mu\nu}X_{\nu}\,,$
$\displaystyle{\mathcal{Q}}_{d\tilde{X}X}$
$\displaystyle=\big{(}\bar{d}_{p}\gamma_{\mu}d_{r}\big{)}\tilde{X}^{\mu\nu}X_{\nu}\,,$
$\displaystyle{\mathcal{Q}}_{DqX^{2}}$
$\displaystyle=\big{(}\bar{q}_{p}\gamma^{\mu}D^{\nu}q_{r}\big{)}X_{\mu}X_{\nu}+\mathrm{h.c.},$
$\displaystyle{\mathcal{Q}}_{DdX^{2}}$
$\displaystyle=\big{(}\bar{d}_{p}\gamma^{\mu}D^{\nu}d_{r}\big{)}X_{\mu}X_{\nu}+\mathrm{h.c.},$
$\displaystyle{\mathcal{Q}}_{dHX^{2}}$
$\displaystyle=(\bar{q}_{p}\sigma_{\mu\nu}d_{r}H)X_{1}^{\mu}X_{2}^{\nu}+\mathrm{h.c.},$
(47)

where $\tilde{X}^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}X_{\rho\sigma}/2$. It is
noted that the rare $B$-meson and kaon decay rates induced by the operators
involving the vector DM field $X_{\mu}$ are divergent in the $m_{X}\to 0$
limit. Such a divergence is caused by the longitudinal polarization of the
vector field and related to the phenomenon that it is impossible to
consistently define a massless limit for a vector field without an active
gauge symmetry Kamenik:2011vy . As in refs. Kamenik:2011vy ; Williams:2011qb ;
He:2022ljo , this singularity is treated by assuming some kind of Higgs
mechanism in the dark sector. As a result, the Wilson coefficients of the
operators in eq. (4.1) should be proportional to some powers of
$(m_{X}/\Lambda)$ and can be redefined as

$\displaystyle{\mathcal{C}}_{i}=\tilde{\mathcal{C}}_{i}\cdot\begin{cases}(m_{X}/\Lambda)^{2}&\text{
for
}{\mathcal{Q}}_{i}={\mathcal{Q}}_{dX^{2}},{\mathcal{Q}}_{DdX^{2}},{\mathcal{Q}}_{DqX^{2}},{\mathcal{Q}}_{dHX^{2}},\\\
(m_{X}/\Lambda)&\text{ for }{\mathcal{Q}}_{i}=\text{others}.\end{cases}$ (48)

These redefined Wilson coefficients $\tilde{\mathcal{C}}_{i}$ will be used in
the following analysis.

We can also include the axion or axion-like particles (ALPs) Peccei:1977hh ;
Peccei:1977ur ; Weinberg:1977ma ; Wilczek:1977pj  in the EFT approach
Georgi:1986df ; Brivio:2017ije ; Chala:2020wvs ; Bauer:2020jbp ; Galda:2021hbr
; Song:2023lxf . They are pseudo Nambu-Goldstone bosons with non-trivial
properties, and can emerge from spontaneous breaking of some global
symmetries. Up to dim-5, there are only two operators,

$\displaystyle{\mathcal{Q}}_{qa}=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}\partial^{\mu}a,\qquad\qquad{\mathcal{Q}}_{da}=\big{(}\bar{d}_{p}\gamma_{\mu}d_{r}\big{)}\partial^{\mu}a,$
(49)

where $a$ denotes the ALP field.

In the MFV hypothesis discussed in section 2, the Wilson coefficients of the
DSMEFT operators should take special flavour structures. Explicitly, by using
the quark currents in eqs. (13)-(15), we obtain

$\displaystyle{\mathcal{C}}^{\rm
MFV}_{i}=\begin{cases}\epsilon_{0}^{i}\hat{\lambda}_{d}+\epsilon_{1}^{i}\Delta_{q}\hat{\lambda}_{d}&\text{for
}{\mathcal{Q}}_{i}={\mathcal{Q}}_{d\phi},{\mathcal{Q}}_{d\phi^{2}},{\mathcal{Q}}_{dHX},{\mathcal{Q}}_{dHX^{2}},{\mathcal{Q}}_{dX^{2}},\\\
\epsilon_{0}^{i}\mathds{1}+\epsilon_{1}^{i}\Delta_{q}&\text{for
}{\mathcal{Q}}_{i}={\mathcal{Q}}_{\phi
q},{\mathcal{Q}}_{q\chi},{\mathcal{Q}}_{qXX},{\mathcal{Q}}_{q\tilde{X}X},{\mathcal{Q}}_{DqX^{2}},{\mathcal{Q}}_{qX},{\mathcal{Q}}_{HqX}^{(1,3)},{\mathcal{Q}}_{qa},\\\
\epsilon_{0}^{i}\mathds{1}&\text{for }{\mathcal{Q}}_{i}={\mathcal{Q}}_{\phi
d},{\mathcal{Q}}_{d\chi},{\mathcal{Q}}_{dXX},{\mathcal{Q}}_{d\tilde{X}X},{\mathcal{Q}}_{DdX^{2}},{\mathcal{Q}}_{dX},{\mathcal{Q}}_{HdX},{\mathcal{Q}}_{da},\end{cases}$
(50)

for the DSMEFT operators discussed above. Here the Wilson coefficient
${\mathcal{C}}^{\rm MFV}_{i}$ are given for the down-type quark currents
expressed in the mass eigenbasis, and they are $3\times 3$ matrices in flavour
space. Generally, the real parameters $\epsilon_{0,1}^{i}$ are different from
each other for different operators ${\mathcal{Q}}_{i}$. We can see that the
operators listed in the last line of the above equation do not generate tree-
level FCNC interactions among the down-type quarks.

###  4.2 DLEFT

Below the electroweak scale, the NP effects induced by the DM particles can be
described by the DLEFT Aebischer:2022wnl . This EFT is obtained by integrating
out all the heavy particles with masses at or above the electroweak scale and
can be considered as an extension of the LEFT Buchalla:1995vs ;
Jenkins:2017jig  by including besides the light SM fields also the light DM
particles. The DLEFT respects the SM $SU(3)_{C}\otimes U(1)_{\rm em}$ gauge
symmetry, and does not contain the Higgs boson. Thus, it can be considered
without reference to the DSMEFT for the light DM and light SM particles
interacting at energies below the electroweak scale. Following the same
convention as in ref. Aebischer:2022wnl , we can write the DLEFT Lagrangian as

$\displaystyle{\mathcal{L}}_{\mathrm{DLEFT}}\supset\frac{1}{\Lambda}\sum_{i}L_{i}{\mathcal{O}}_{i}^{(5)}+\frac{1}{\Lambda^{2}}\sum_{j}L_{j}{\mathcal{O}}_{j}^{(6)},$
(51)

up to dim-6, where ${\mathcal{O}}_{i}^{(d)}$ denote the dim-$d$ effective
operators and $L_{i}$ the corresponding Wilson coefficients.

Analogous to the DSMEFT case, we can write down all the relevant DLEFT
operators involving the scalar, fermionic and vector DM fields. For the scalar
DM, we have

$\displaystyle{\mathcal{O}}_{d\phi}$
$\displaystyle=(\bar{d}_{Lp}d_{Rr})\phi+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{\phi d}^{L}$
$\displaystyle=(\bar{d}_{Lp}\gamma_{\mu}d_{Lr})\big{(}i\phi_{1}\overleftrightarrow{\partial^{\mu}}\phi_{2}\big{)},$
$\displaystyle{\mathcal{O}}_{d\phi^{2}}$
$\displaystyle=(\bar{d}_{Lp}d_{Rr})\phi^{2}+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{\phi d}^{R}$
$\displaystyle=(\bar{d}_{Rp}\gamma_{\mu}d_{Rr})\big{(}i\phi_{1}\overleftrightarrow{\partial^{\mu}}\phi_{2}\big{)},$
(52)

where $d_{L}$ and $d_{R}$ denote the left-handed and right-handed down-type
quarks, respectively. For the fermionic DM, the relevant operators are given
by

$\displaystyle{\mathcal{O}}_{d\chi}^{V,\,LR}$
$\displaystyle=(\bar{d}_{Lp}\gamma_{\mu}d_{Lr})(\bar{\chi}_{a}\gamma^{\mu}\chi_{b}),$
$\displaystyle{\mathcal{O}}_{d\chi}^{V,RR}$
$\displaystyle=(\bar{d}_{Rp}\gamma_{\mu}d_{Rr})(\bar{\chi}_{a}\gamma^{\mu}\chi_{b}),$
$\displaystyle{\mathcal{O}}_{d\chi}^{S,LR}$
$\displaystyle=(\bar{d}_{Rp}d_{Lr})(\chi_{a}^{T}C\chi_{b})+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{d\chi}^{S,RR}$
$\displaystyle=(\bar{d}_{Lp}d_{Rr})(\chi_{a}^{T}C\chi_{b})+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{d\chi}^{T,RR}$
$\displaystyle=(\bar{d}_{Lp}\sigma^{\mu\nu}d_{Rr})(\chi_{a}^{T}C\sigma_{\mu\nu}\chi_{b})+\mathrm{h.c.},$
(53)

with the charge conjugate operator defined by $C=i\gamma^{2}\gamma^{0}$. For
the vector DM, there is only one operator,

$\displaystyle{\mathcal{O}}_{dX}^{T}$
$\displaystyle=(\bar{d}_{Lp}\sigma_{\mu\nu}d_{Rr})X_{a}^{\mu\nu}+\mathrm{h.c.},$
(54)

by using the field strength $X_{\mu\nu}$ as the building block. However,
including $X_{\mu}$ as additional building block, the following operators are
obtained:

$\displaystyle{\mathcal{O}}_{dX}^{L}$
$\displaystyle=\big{(}\bar{d}_{Lp}\gamma_{\mu}d_{Lr}\big{)}X^{\mu},$
$\displaystyle{\mathcal{O}}_{dX}^{R}$
$\displaystyle=\big{(}\bar{d}_{Rp}\gamma_{\mu}d_{Rr}\big{)}X^{\mu},$
$\displaystyle{\mathcal{O}}_{dXX}^{L}$
$\displaystyle=\big{(}\bar{d}_{Lp}\gamma_{\mu}d_{Lr}\big{)}X^{\mu\nu}X_{\nu},$
$\displaystyle{\mathcal{O}}_{dXX}^{R}$
$\displaystyle=\big{(}\bar{d}_{Rp}\gamma_{\mu}d_{Rr}\big{)}X^{\mu\nu}X_{\nu},$
$\displaystyle{\mathcal{O}}_{d\tilde{X}X}^{L}$
$\displaystyle=\big{(}\bar{d}_{Lp}\gamma_{\mu}d_{Lr}\big{)}\tilde{X}^{\mu\nu}X_{\nu},$
$\displaystyle{\mathcal{O}}_{d\tilde{X}X}^{R}$
$\displaystyle=\big{(}\bar{d}_{Rp}\gamma_{\mu}d_{Rr}\big{)}\tilde{X}^{\mu\nu}X_{\nu},$
$\displaystyle{\mathcal{O}}_{DdX^{2}}^{L}$
$\displaystyle=\big{(}\bar{d}_{Lp}\gamma^{\nu}D^{\mu}d_{Lr}\big{)}X_{\mu}X_{\nu}+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{DdX^{2}}^{R}$
$\displaystyle=\big{(}\bar{d}_{Rp}\gamma^{\nu}D^{\mu}d_{Rr}\big{)}X_{\mu}X_{\nu}+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{dX^{2}}$
$\displaystyle=\big{(}\bar{d}_{Lp}d_{Rr}\big{)}X_{\mu}X^{\mu}+\mathrm{h.c.},$
$\displaystyle{\mathcal{O}}_{dX^{2}}^{T}$
$\displaystyle=\big{(}\bar{d}_{Lp}\sigma_{\mu\nu}d_{Rr}\big{)}X_{1}^{\mu}X_{2}^{\nu}+\mathrm{h.c.},$
(55)

where $D_{\mu}=\partial_{\mu}+ig_{s}T^{a}G_{\mu}^{a}+ieQA_{\mu}$ denotes the
covariant derivative of the $SU(3)_{C}\otimes U(1)_{\rm em}$ gauge symmetry.
Here we have chosen the basis where the combinations of $X_{\mu\nu}$ and
$X_{\mu}$ are similar as in the DSMEFT operators of eq. (4.1). One can prove
that this basis is equivalent to the one adopted in ref. He:2022ljo . In
addition, there are two operators involving the ALP field, which take the form

$\displaystyle{\mathcal{O}}_{da}^{L}=\big{(}\bar{d}_{Lp}\gamma_{\mu}d_{Lr}\big{)}\partial^{\mu}a,\qquad\qquad{\mathcal{O}}_{da}^{R}=\big{(}\bar{d}_{Rp}\gamma_{\mu}d_{Rr}\big{)}\partial^{\mu}a.$
(56)

Although the DLEFT can be considered without reference to the DSMEFT, the
former should arise as the low-energy limit of the latter, if we start from
the DSMEFT. At the electroweak scale $\mu_{\rm EW}$, the tree-level matching
conditions between these two EFTs are summarized below:

  * •

Scalar DM:

$\displaystyle L_{\phi d}^{L}$ $\displaystyle={\mathcal{C}}_{\phi q},$
$\displaystyle L_{\phi d}^{R}$ $\displaystyle={\mathcal{C}}_{\phi d},$
$\displaystyle L_{d\phi^{2}}$
$\displaystyle=\frac{v}{\sqrt{2}\Lambda}{\mathcal{C}}_{d\phi^{2}},$
$\displaystyle L_{d\phi}$
$\displaystyle=\frac{v}{\sqrt{2}\Lambda}{\mathcal{C}}_{d\phi}.$ (57)

  * •

Fermionic DM:

$\displaystyle L_{d\chi}^{V,\,LR}$ $\displaystyle={\mathcal{C}}_{q\chi},$
$\displaystyle L_{d\chi}^{V,\,RR}$ $\displaystyle={\mathcal{C}}_{d\chi},$
$\displaystyle L_{d\chi}^{S,\,LR}$ $\displaystyle=0,$ $\displaystyle
L_{d\chi}^{S,\,RR}$ $\displaystyle=0,$ $\displaystyle L_{d\chi}^{T,\,RR}$
$\displaystyle=0.$ (58)

  * •

Vector DM:

$\displaystyle
L_{dX}^{T}=\frac{v}{\sqrt{2}\Lambda}{\mathcal{C}}_{dHX},L_{dX}^{L}={\mathcal{C}}_{qX}+\frac{v^{2}}{2\Lambda^{2}}\big{(}{\mathcal{C}}_{HqX}^{(1)}+{\mathcal{C}}_{HqX}^{(3)}\big{)},L_{dX}^{R}={\mathcal{C}}_{dX}+\frac{v^{2}}{2\Lambda^{2}}{\mathcal{C}}_{HdX},$
$\displaystyle L_{dXX}^{L}$ $\displaystyle={\mathcal{C}}_{qXX},$
$\displaystyle L_{d\tilde{X}X}^{L}$
$\displaystyle={\mathcal{C}}_{q\tilde{X}X},$ $\displaystyle L_{DdX^{2}}^{L}$
$\displaystyle={\mathcal{C}}_{DqX^{2}},$ $\displaystyle L_{dX^{2}}$
$\displaystyle=\frac{v}{\sqrt{2}\Lambda}{\mathcal{C}}_{dX^{2}},$
$\displaystyle L_{dXX}^{R}$ $\displaystyle={\mathcal{C}}_{dXX},$
$\displaystyle L_{d\tilde{X}X}^{R}$
$\displaystyle={\mathcal{C}}_{d\tilde{X}X},$ $\displaystyle L_{DdX^{2}}^{R}$
$\displaystyle={\mathcal{C}}_{DdX^{2}},$ $\displaystyle L_{dX^{2}}^{T}$
$\displaystyle=\frac{v}{\sqrt{2}\Lambda}{\mathcal{C}}_{dHX^{2}}.$ (59)

  * •

ALP:

$\displaystyle L_{da}^{L}={\mathcal{C}}_{qa},\qquad\qquad
L_{da}^{R}={\mathcal{C}}_{da}.$ (60)

Here all the Wilson coefficients are given in the mass eigenbasis. In the MFV
hypothesis, the Wilson coefficients in DSMEFT take the form of eq. (50).

Since the above matching conditions are given at the electroweak scale, i.e.,
$L_{i}\equiv L_{i}(\mu_{\rm EW})$ and
${\mathcal{C}}_{i}\equiv{\mathcal{C}}_{i}(\mu_{\rm EW})$, the renormalization
group (RG) evolution should be performed in the DLEFT to obtain the
corresponding Wilson coefficients at the typical scales of $B$-meson and kaon
decays. In the DLEFT, as the DM particles are color singlet, we need only
consider the QCD RG evolution of the quark current in each DLEFT operator.
Furthermore, the different DLEFT operators do not mix under the QCD RG
running. Accordingly, the RG equations take the form

$\displaystyle\frac{dL_{i}(\mu)}{d\ln\mu}=\frac{\alpha_{s}}{4\pi}\gamma_{i}\,L_{i}(\mu),$
(61)

where $L_{i}(\mu)$ denote the Wilson coefficients of the operators
${\mathcal{O}}_{i}$ at the scale $\mu$. $\gamma_{i}$ denote the anomalous
dimension matrices (ADM) for the operators ${\mathcal{O}}_{i}$. They can be
obtained from the ADM of the semi-leptonic operators in the LEFT
Buchalla:1995vs ; Jenkins:2017dyc . For the DLEFT operators with a non-zero
matching condition, the values of these ADM are given by

$\displaystyle\gamma_{i}=\begin{cases}-6C_{F}&\text{for
}{\mathcal{O}}_{i}={\mathcal{O}}_{d\phi},{\mathcal{O}}_{d\phi^{2}},{\mathcal{O}}_{dX^{2}},\\\
+2C_{F}&\text{for
}{\mathcal{O}}_{i}={\mathcal{O}}_{dX}^{T},{\mathcal{O}}_{dX^{2}}^{T},\\\
0&\text{for }{\mathcal{O}}_{i}={\mathcal{O}}_{\phi
d}^{L,R},{\mathcal{O}}_{d\chi}^{V,LR},{\mathcal{O}}_{d\chi}^{V,RR},{\mathcal{O}}_{dX}^{L,R},{\mathcal{O}}_{dXX}^{L,R},\,{\mathcal{O}}_{d\tilde{X}X}^{L,R},{\mathcal{O}}_{DdX^{2}}^{L,R},\end{cases}$
(62)

with $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and $N_{c}=3$. Numerically, we find

$\displaystyle L_{i}(4.8{\,\rm GeV})$
$\displaystyle=L_{i}(m_{Z})\cdot\begin{cases}1.37&\text{for
}{\mathcal{O}}_{i}={\mathcal{O}}_{d\phi},{\mathcal{O}}_{d\phi^{2}},{\mathcal{O}}_{dX^{2}},\\\
0.90&\text{for
}{\mathcal{O}}_{i}={\mathcal{O}}_{dX}^{T},{\mathcal{O}}_{dX^{2}}^{T},\end{cases}$
(63) $\displaystyle L_{i}(2.0{\,\rm GeV})$
$\displaystyle=L_{i}(m_{Z})\cdot\begin{cases}1.59&\text{for
}{\mathcal{O}}_{i}={\mathcal{O}}_{d\phi},{\mathcal{O}}_{d\phi^{2}},{\mathcal{O}}_{dX^{2}},\\\
0.86&\text{for }{\mathcal{O}}_{i}={\mathcal{O}}_{dX},\end{cases}$

which can be used in the $B$-meson and kaon decays, respectively.

###  4.3 Observables

Within the DLEFT, by using the similar treatment as for the $d_{i}\to
d_{j}\nu\bar{\nu}$ transitions discussed in section 3.2, one can calculate the
various $d_{i}\to d_{j}+\text{DM}$ processes. Generally, the amplitudes of
these decays can be factorized into the product of the Wilson coefficient and
the hadronic matrix element of the corresponding effective operator that can
be parameterized form factors or decay constants. In appendix A, we summarize
the relevant hadronic matrix elements used in this paper. For the explicit
expressions of the decay rates, we refer the readers to refs. Bird:2004ts ;
Li:2019cbk ; Kamenik:2011vy ; He:2020jly ; Geng:2020seh ; Li:2021sqe ;
Kling:2022uzy  for the $d_{i}\to d_{j}+\phi\,(+\phi)$ decays, refs. Su:2019tjn
; Kamenik:2011vy ; Li:2020dpc ; Felkl:2021uxi  for the $d_{i}\to
d_{j}+\chi+\chi$ decay, refs. Kamenik:2011vy ; Li:2021sqe ; He:2022ljo  for
the $d_{i}\to d_{j}+X\,(+X)$ decays, and refs. MartinCamalich:2020dfe ;
Bauer:2021mvw ; Guerrera:2022ykl  for the $d_{i}\to d_{i}+a$ decays. All the
relevant analytical expressions have been implemented in our recently
developed package HadronToNP HadronToNP:2024 , which will be used in the
following numerical analysis.

##  5 Numerical analysis

In this section, we proceed to present our numerical results and discussions.
In table 1, we list the main input parameters used in our numerical analysis.
Table 2 summaries our SM predictions and the up-to-date experimental
measurements of the relevant observables. We take $\Lambda_{\rm NP}=1{\,\rm
TeV}$ in the SMEFT and DSMEFT. In the following, we investigate two
explanations to the recent Belle II measurement of ${\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})$ Belle-II:2023esi . In section 5.1, we consider the $b\to
s+\nu\bar{\nu}$ transition in the SMEFT. In section 5.2, the missing energy
signals in the Belle II measurement are interpreted in the DSMEFT as
contributions from DM particles.

Input | Value | Unit | Reference
---|---|---|---
$m_{t}^{\text{pole}}$ | $172.69\pm 0.30$ | GeV | Workman:2022ynf
$|V_{cb}|$(semi-leptonic) | $41.15\pm 0.34\pm 0.45$ | $10^{-3}$ | CKMfitter
$|V_{ub}|$(semi-leptonic) | $3.88\pm 0.08\pm 0.21$ | $10^{-3}$ | CKMfitter
$|V_{us}|f^{K\to\pi}_{+}(0)$ | $0.2165\pm 0.0004$ |  | CKMfitter
$\gamma$ | $72.1^{+5.4}_{-5.7}$ | ∘ | CKMfitter
$f^{K\to\pi}_{+}(0)$ | $0.9675\pm 0.0009\pm 0.0023$ |  | CKMfitter
$\Delta m_{21}^{2}$ |  $7.41^{+0.21}_{-0.20}$ ($7.41^{+0.21}_{-0.20}$ ) | $10^{-5}\cdot\text{eV}^{2}$ | Workman:2022ynf
$\Delta m_{32}^{2}$ |  $2.437^{+0.028}_{-0.027}$ ($-2.498^{+0.032}_{-0.025}$) | $10^{-3}\cdot\text{eV}^{2}$ | Workman:2022ynf
$\theta_{12}$ |  $33.41^{+0.75}_{-0.72}$ ($33.41^{+0.75}_{-0.72}$) | ∘ | Workman:2022ynf
$\theta_{23}$ |  $49.1^{+1.0}_{-1.3}$ ($49.5^{+0.9}_{-1.2}$) | ∘ | Workman:2022ynf
$\theta_{13}$ |  $8.54^{+0.11}_{-0.12}$ ($8.57^{+0.12}_{-0.11}$) | ∘ | Workman:2022ynf
$\delta_{\rm CP}$ |  $197^{+42}_{-25}$ ($286^{+27}_{-32}$) | ∘ | Workman:2022ynf
Table 1: Main input parameters used in our numerical analysis. The values (in brackets) of the parameters in the lepton sector correspond to the NO (IO) of neutrino masses. Observable | SM | Exp | Unit | Sensitivity
---|---|---|---|---
${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ | $4.16\pm 0.57$ |  $23\pm 5^{+5}_{-4}$ belleIItalk2023  | $10^{-6}$ | 19% Belle-II:2022cgf
${\mathcal{B}}(B^{0}\to K^{0}\nu\bar{\nu})$ | $3.85\pm 0.52$ |  $<26$ Belle:2017oht  | $10^{-6}$ | 87% Belle-II:2022cgf
${\mathcal{B}}(B^{+}\to K^{*+}\nu\bar{\nu})$ | $9.70\pm 0.94$ |  $<61$ Belle:2017oht  | $10^{-6}$ | 75% Belle-II:2022cgf
${\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})$ | $9.00\pm 0.87$ |  $<18$ Belle:2017oht  | $10^{-6}$ | 40% Belle-II:2022cgf
${\mathcal{B}}(B_{s}\to\phi\nu\bar{\nu})$ | $9.93\pm 0.72$ |  $<5400$ DELPHI:1996ohp  | $10^{-6}$ |  2% Li:2022tov
${\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})$ | $1.40\pm 0.18$ |  $<140$ Belle:2017oht  | $10^{-7}$ |
${\mathcal{B}}(B^{0}\to\pi^{0}\nu\bar{\nu})$ | $6.52\pm 0.85$ |  $<900$ Belle:2017oht  | $10^{-8}$ |
${\mathcal{B}}(B^{+}\to\rho^{+}\nu\bar{\nu})$ | $4.06\pm 0.79$ |  $<300$ Belle:2017oht  | $10^{-7}$ |
${\mathcal{B}}(B^{0}\to\rho^{0}\nu\bar{\nu})$ | $1.89\pm 0.36$ |  $<400$ Belle:2017oht  | $10^{-7}$ |
${\mathcal{B}}(K^{+}\to\pi^{+}\nu\bar{\nu})$ | $8.42\pm 0.61$ |  $10.6^{+4.0}_{-3.4}\pm 0.9$ NA62:2021zjw  | $10^{-11}$ |
${\mathcal{B}}(K_{L}\to\pi^{0}\nu\bar{\nu})$ | $3.41\pm 0.45$ |  $<300$ KOTO:2018dsc  | $10^{-11}$ |
${\mathcal{B}}(B_{s}\to\nu\bar{\nu})$ | $\approx 0$ |  $<5.9$ Alonso-Alvarez:2023mgc  | $10^{-4}$ |  $<1.1\times 10^{-5}$ Belle-II:2018jsg
${\mathcal{B}}(B^{0}\to\nu\bar{\nu})$ | $\approx 0$ |  $<1.4$ Alonso-Alvarez:2023mgc  | $10^{-4}$ |  $<5.0\times 10^{-6}$ Belle-II:2018jsg
Table 2: Our SM predictions and the experimental measurements of the relevant
observables. Upper limits are all given at $90\%$ confidence level (CL). The
last column shows the expected sensitivities at the Belle II ($5{\,\rm
ab}^{-1}$) or at the CEPC (Tera-$Z$ phase). The sensitivities at the FCC-ee
are very similar to that at the CEPC, and hence are not given here.

###  5.1 $\boldsymbol{M_{1}\to M_{2}+\nu\bar{\nu}}$

The recent Belle II measurement of the $B^{+}\to K^{+}+\text{inv}$ decay could
be explained by NP contributions to the $b\to s\nu\bar{\nu}$ transition. In
this case, the theoretical prediction of the related $M_{1}\to
M_{2}+\text{inv}$ processes takes the form

$\displaystyle{\mathcal{B}}_{\rm th}(M_{1}\to
M_{2}+\text{inv})=\sum_{i,j=e,\mu,\tau}{\mathcal{B}}(M_{1}\to
M_{2}\nu_{i}\bar{\nu}_{j}).$ (64)

The NP contributions to these processes can be generally described in the
framework of SMEFT. We refer to refs. Allwicher:2023xba ; Bause:2023mfe ;
Chen:2024jlj  for the recent relevant studies. In the following, we also focus
on the SMEFT but with the MFV hypothesis.

In the SMEFT, as mentioned in section 3.2.1, one interesting implication of
the MFV hypothesis is that the $d_{i}\to d_{j}\nu\bar{\nu}$ transitions are
all governed by only one effective operator at low energy. As a consequence,
the ratio of the branching fractions of any two $b\to s\nu\bar{\nu}$ decays
should be independent of the NP contribution. Considering as an example the
$B^{+}\to K^{+}\nu\bar{\nu}$ and $B^{0}\to K^{*0}\nu\bar{\nu}$ decays, we
obtain

$\displaystyle\frac{{\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})}{{\mathcal{B}}(B^{0}\to
K^{*0}\nu\bar{\nu})}=\frac{{\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})_{\rm
SM}}{{\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})_{\rm SM}}=0.46\pm 0.07,$ (65)

where summation over the neutrino flavours has been taken into account. In
this ratio, all the relevant SMEFT operators, i.e. ${\mathcal{Q}}_{lq}^{(1)}$
and ${\mathcal{Q}}_{lq}^{(3)}$, have been included simultaneously.444The ratio
still holds after including the operators ${\mathcal{Q}}_{Hq}^{(1,3)}$,
although they are not considered in our analysis as mentioned in section 3.1.
When deriving this ratio, it is noted that only the MFV hypothesis in the
quark sector is required to eliminate the right-handed operator
${\mathcal{O}}_{R}$. Therefore, the ratio still holds even by assuming generic
couplings in the lepton sector, i.e., without the leptonic MFV. In addition,
as both the CKM and Wilson coefficients involved are cancelled in this ratio,
the theoretical uncertainty is reduced compared to that of the individual
branching ratios. If future experimental measurements find any deviation from
the theoretical value, it could imply new sources of flavour violation beyond
the Yukawa couplings. By using the above ratio, the Belle $90\%$ CL bound
${\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})<18\times 10^{-6}$ would imply that

$\displaystyle{\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})_{\rm MFV}<10.5\times
10^{-6},$ (66)

which is obviously lower that the $90\%$ CL lower bound of the Belle II
measurement, ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})>12.5\times 10^{-6}$.
Or inversely, the Belle II measurement ${\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})=(23_{-6}^{+7})\times 10^{-6}$ would indicate that

$\displaystyle{\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})_{\rm
MFV}=\big{(}50_{-16}^{+17}\big{)}\times 10^{-6},$ (67)

which is incompatible with the Belle $90\%$ CL bound ${\mathcal{B}}(B^{0}\to
K^{*0}\nu\bar{\nu})<18\times 10^{-6}$. As a conclusion, the recent Belle II
measurement of the $B^{+}\to K^{+}\nu\bar{\nu}$ decay cannot be explained in
the SMEFT with MFV hypothesis. However, the MFV prediction of
${\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})$ is well within the reach of Belle
II at $5{\,\rm ab}^{-1}$ (cf. table 2), which can be used as a further test of
the MFV hypothesis considered here.

Similarly, the ratio of the branching fractions of a $b\to s\nu\bar{\nu}$ and
a $b\to d\nu\bar{\nu}$ decay is also independent of the NP contribution in the
SMEFT with MFV hypothesis. As an illustration, let us consider the $B^{+}\to
K^{+}\nu\bar{\nu}$ and $B^{+}\to\pi^{+}\nu\bar{\nu}$ decays, and obtain

$\displaystyle\frac{{\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})}{{\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})}=\frac{{\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})_{\rm SM}}{{\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})_{\rm
SM}}=29.7\pm 5.6,$ (68)

where summation over the neutrino flavours has been included. Taking as input
the Belle II measurement in eq. (1), this ratio implies that

$\displaystyle{\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})_{\rm
MFV}=\big{(}7.8_{-2.6}^{+2.8}\big{)}\times 10^{-7},$ (69)

which shows a $2.3\sigma$ deviation from the SM prediction
${\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})_{\rm SM}=(1.40\pm 0.18)\times
10^{-7}$, but is still compatible with the current $90\%$ CL experimental
bound ${\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})<140\times 10^{-7}$. Of
course, similar ratios can also be constructed for other processes to test the
framework used here.

In the MFV framework, all the quark FCNC or the LFV processes are related with
each other. Explicitly, for the operators ${\mathcal{Q}}_{lq}^{(i)}$
($i=1,3$), all the quark FCNC transitions are governed by the two parameter
products $\epsilon_{1}^{(i)}\kappa_{0}^{(i)}$ and
$\epsilon_{1}^{(i)}\kappa_{1}^{(i)}$. In order to constrain them, we consider
several quark (and lepton) FCNC processes, including the $B\to
K^{(*)}\nu\bar{\nu}$, $B\to K^{(*)}e\mu$, $B\to\pi e\mu$, $B_{s}\to e\mu$,
$B_{s}\to\mu^{+}\mu^{-}$ as well as $K\to\pi\nu\bar{\nu}$ decays. Besides the
input parameters listed in table 1, the lepton number breaking scale
$\Lambda_{\rm LN}=10^{14}{\,\rm GeV}$ and the lightest neutrino mass
$m_{\nu_{1}}=0.2\,\text{eV}$ are used in the leptonic MFV. Numerically, we
find that the dominant constraints come from the $B^{+}\to K^{+}\nu\bar{\nu}$,
$B^{0}\to K^{*0}\nu\bar{\nu}$, $B^{0}\to K^{*0}e^{-}\mu^{+}$ and
$B_{s}\to\mu^{+}\mu^{-}$ decays. The allowed parameter regions of
$\epsilon_{1}^{(i)}\kappa_{0}^{(i)}$ and $\epsilon_{1}^{(i)}\kappa_{1}^{(i)}$
are shown in figure 1. We can see that there is a small gap between the
allowed regions of the $B^{+}\to K^{+}\nu\bar{\nu}$ and $B^{0}\to
K^{*0}\nu\bar{\nu}$ decays. This reflects the discrepancy between the MFV
prediction in eq. (67) and the current Belle bound on ${\mathcal{B}}(B^{0}\to
K^{*0}\nu\bar{\nu})$.

Figure 1: Constraints on the MFV parameter products
$(\epsilon_{1}^{(i)}\kappa_{0}^{(i)},\epsilon_{1}^{(i)}\kappa_{1}^{(i)})$ for
${\mathcal{Q}}_{i}={\mathcal{Q}}_{lq}^{(1)}$ (left) and
${\mathcal{Q}}_{i}={\mathcal{Q}}_{lq}^{(3)}$ (right) in the cases of NO
(upper) and IO (lower) of neutrino masses. The allowed regions by the
measurement of ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ are at $2\sigma$
level, while the regions by ${\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})$,
${\mathcal{B}}(B^{0}\to K^{*0}\mu^{+}e^{-})$ and
${\mathcal{B}}(B_{s}\to\mu^{+}\mu^{-})$ at $90\%$ CL. The allowed regions by
all the last three observables are shown in blue. The black points $(0,0)$
denote the SM case.

Since the uncertainty of the recent Belle II measurement is still quite large,
we use instead all the processes except the $B^{+}\to K^{+}\nu\bar{\nu}$ decay
to derive the allowed parameter space. As shown in figure 1, the upper bounds
on the magnitudes of the products $\epsilon_{1}^{(i)}\kappa_{0}^{(i)}$ and
$\epsilon_{1}^{(i)}\kappa_{1}^{(i)}$ for the two operators
${\mathcal{Q}}_{lq}^{(1,3)}$ are all around $0.05$ in the cases of NO and IO
of neutrino masses, which are much smaller than the parameter regions required
to explain the $B^{+}\to K^{+}\nu\bar{\nu}$ excess. Using the allowed
parameter regions for ${\mathcal{Q}}_{lq}^{(i)}$, we can derive the following
numerical predictions:

$\displaystyle 3.4\,(3.5)\times 10^{-6}<{\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})_{\rm MFV}$ $\displaystyle<5.0\,(4.9)\times 10^{-6},$
$\displaystyle 7.4\,(7.5)\times 10^{-6}<{\mathcal{B}}(B^{0}\to
K^{*0}\nu\bar{\nu})_{\rm MFV}$ $\displaystyle<10.8\,(10.6)\times 10^{-6},$
$\displaystyle 1.1\,(1.2)\times
10^{-7}<\,{\mathcal{B}}(B^{+}\to\pi^{+}\nu\bar{\nu})_{\rm MFV}\,$
$\displaystyle<1.7\,(1.7)\times 10^{-7},$ (70)

which correspond the NO (IO) of neutrino masses. The numerical results show
less than $0.5\%$ difference between ${\mathcal{Q}}_{lq}^{(1)}$ and
${\mathcal{Q}}_{lq}^{(3)}$, implying that both of these two operators coincide
with each other under the above bounds. From these bounds, we can see that the
NP contributions to the $B\to K\nu\bar{\nu}$ and $B\to\pi\nu\bar{\nu}$ decays
should be less than $\sim 20\%$ of their SM predictions. It is noted that the
leptonic MFV has been used to derive the allowed parameter space and the above
bounds. Therefore, the bound on ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ is
more stringent than the one in eq. (66) derived without the leptonic MFV.

###  5.2 $\boldsymbol{M_{1}\to M_{2}+\text{DM}}$

The large excess of ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ reported by
Belle II can also be explained by light NP particles that contribute to the
signal as $B^{+}\to K^{+}+\text{missing energy}$. In this case, the
theoretical predictions of the $M_{1}\to M_{2}+\text{inv}$ decay can be
written as

$\displaystyle{\mathcal{B}}_{\rm th}(M_{1}\to
M_{2}+\text{inv})={\mathcal{B}}_{\rm SM}(M_{1}\to
M_{2}+\nu\bar{\nu})+{\mathcal{B}}_{\rm NP}(M_{1}\to M_{2}+\text{DM}).$ (71)

As discussed in section 4, the effects of light NP particles in the second
term can be generally described in the EFT framework. In the following, we
investigate this explanation in the DSMEFT framework, in which the Wilson
coefficients are taken to be the most general form or satisfy the MFV
hypothesis.

####  5.2.1 General flavour structure

In the DSMEFT with general flavour structure, only the $b\to s+\text{inv}$
processes are directly related to the $B^{+}\to K^{+}+\text{inv}$ decays,
which include $B^{0}\to K^{0}+\text{inv}$, $B^{0,+}\to K^{*\,0,+}+\text{inv}$,
$B_{s}\to\phi+\text{inv}$ and $B_{s}\to\text{inv}$ decays. In the following,
we investigate constraints on the DSMEFT operators from the experimental data
on these decays as summarized in table 2, and discuss future prospects at
Belle II, CEPC and FCC-ee. As discussed in section 4.1, seven DSMEFT operators
are not hermitian. In order to compare with the MFV analysis in the next
subsection, contributions from the hermitian conjugation of these operators
are neglected in the numerical analysis in this subsection, i.e.,
${\mathcal{C}}_{ji}=0$ is taken for $d_{i}\to d_{j}$ processes.555Let us
consider ${\mathcal{Q}}_{d\phi^{2}}$ as an example. For the $b\to s\phi\phi$
process, $[{\mathcal{C}}_{d\phi^{2}}]_{23}\neq 0$ and
$[{\mathcal{C}}_{d\phi^{2}}]_{32}=0$ are taken in the general DSMEFT analysis.
In the MFV hypothesis, the coupling matrix eq. (18) implies that
$[{\mathcal{C}}_{d\phi^{2}}]_{23}\gg[{\mathcal{C}}_{d\phi^{2}}]_{32}$, which
results in the approximation $[{\mathcal{C}}_{d\phi^{2}}]_{32}\approx 0$.  For
example, contributions to the $b\to s\phi\phi$ processes from the operator
${\mathcal{Q}}_{d\phi^{2}}$ are calculated by taking
$[{\mathcal{C}}_{d\phi^{2}}]_{23}\neq 0$ and
$[{\mathcal{C}}_{d\phi^{2}}]_{32}=0$. Two exceptions are the operators
${\mathcal{Q}}_{DqX^{2}}$ and ${\mathcal{Q}}_{DdX^{2}}$. Consider the former
as an example, we take
$[{\mathcal{C}}_{DqX^{2}}]_{ji}=[{\mathcal{C}}_{DqX^{2}}]_{ij}^{*}$ for the
$d_{i}\to d_{j}$ processes to compare with the MFV analysis, where its MFV
coupling matrix in eq. (17) also implies a hermitian
${\mathcal{C}}_{DqX^{2}}$. In this case, the hadronic matrix elements of the
weak currents can be parameterized in terms of the transition form factors of
the vector current as discussed in appendix A.2.

Figure 2: Constraints on the Wilson coefficients of the operators involving
scalar DM as a function of the scalar mass $m_{\phi}$. The blue band denotes
the $2\sigma$ allowed region by the recent Belle II measurement of the
$B^{+}\to K^{+}+\text{inv}$ decay. The upper bounds from other measurements
are shown by the solid lines, while the expected sensitivities at Belle II and
CEPC by the dotted lines. The sensitivities at FCC-ee are very similar to that
at CEPC and hence are not shown. For the operators ${\mathcal{Q}}_{\phi q}$
and ${\mathcal{Q}}_{\phi d}$, we have considered four different cases for the
DM masses, as indicated in the last four plots. Figure 3: Same as in figure 2,
but for the operators involving ALP.

For the scalar DM, there are three operators (${\mathcal{Q}}_{d\phi^{2}}$,
${\mathcal{Q}}_{\phi q}$ and ${\mathcal{Q}}_{\phi d}$, cf. eq. (4.1))
contributing to the three-body $d_{i}\to d_{j}\phi\phi$ and one operator
(${\mathcal{Q}}_{d\phi}$, cf. eq. (4.1)) to the two-body $d_{i}\to d_{j}\phi$
decay. Allowed regions of the Wilson coefficients for these four operators are
shown in figure 2. From this figure, we make the following observations:

  * •

The $B^{0}\to K^{*0}+\text{inv}$ decay provides the strongest constraints on
the parameter regions allowed by the Belle II measurement of the $B^{+}\to
K^{+}+\text{inv}$ decay. Especially, for the operators ${\mathcal{Q}}_{\phi
q}$ and ${\mathcal{Q}}_{\phi d}$, the regions of $m_{\phi_{2}}\lesssim
1.5{\,\rm GeV}$ ($m_{\phi_{1}}<m_{\phi_{2}}$ is assumed without loss of
generality) are almost excluded by the decay. For the operator
${\mathcal{Q}}_{d\phi}$, on the other hand, large part of the parameter region
with $m_{\phi}\lesssim 3.8{\,\rm GeV}$ is excluded.

  * •

Another stronger constraint arises from the $B_{s}\to\text{inv}$ decay, which
benefits from the large phase space. As a result, it excludes the parameter
region for large DM mass, e.g., the region with $m_{\phi}\gtrsim 2.3{\,\rm
GeV}$ is excluded for the operator ${\mathcal{Q}}_{d\phi^{2}}$. For the
operators ${\mathcal{Q}}_{\phi q}$ and ${\mathcal{Q}}_{\phi d}$, however, the
$B_{s}\to\phi_{1}\phi_{2}$ amplitude vanishes in the case of
$m_{\phi_{1}}=m_{\phi_{2}}$ due to the derivative term
$\phi_{1}\overleftrightarrow{\partial^{\mu}}\phi_{2}$ in the operators.
Therefore, only the parameter regions with $m_{\phi_{1}}\neq m_{\phi_{2}}$
receive constraint from the $B_{s}\to\text{inv}$ decay.

  * •

Compared to the $B^{0}\to K^{*0}+\text{inv}$ and $B_{s}\to\text{inv}$ decays,
The $B^{0}\to K^{0}+\text{inv}$ and $B^{+}\to K^{*+}+\text{inv}$ decays
provide almost no further constraints. The upper bounds from these two decays
are comparable to or weaker than the upper limits of the $2\sigma$ allowed
region required to explain the $B^{+}\to K^{+}+\text{inv}$ excess.

  * •

All the processes put very weak constraints on the parameter regions with
$m_{\phi}\lesssim 2.3{\,\rm GeV}$ for the operator
${\mathcal{Q}}_{d\phi^{2}}$.

In addition, constraints on the operators involving ALP are shown in figure 3,
which are very similar to the ones for the operator ${\mathcal{Q}}_{d\phi}$
but with much larger allowed Wilson coefficients.

For the fermionic DM, there are only two relevant operators
${\mathcal{Q}}_{q\chi}$ and ${\mathcal{Q}}_{d\chi}$ (cf. eq. (45)) in the
DSMEFT. The experimental constraints are shown in figure 4. In the cases of
$m_{\chi}\lesssim 0.6{\,\rm GeV}$ and $m_{\chi}\gtrsim 2.3{\,\rm GeV}$, the
parameter spaces required to account for the Belle II measurement are already
excluded by the $B^{0}\to K^{*0}+\text{inv}$ and $B_{s}\to\text{inv}$ decays,
respectively. In the range of $0.6\lesssim m_{\chi}\lesssim 2.3{\,\rm GeV}$,
however, large part of the parameter space survives all the constraints. This
is quite different from the case of the SMEFT operator
${\mathcal{Q}}_{lq}^{(1)}$ discussed in section 5.1, due to the massive
fermionic DM fields in the operators ${\mathcal{Q}}_{q\chi}$ and
${\mathcal{Q}}_{d\chi}$. In the massless limit, their effects on the decay
width are the same as that of the operator
${\mathcal{Q}}_{lq}^{(1)}$.666Although the leptonic and DM currents in these
operators have different chiralities, the total decay width does not depend on
the chirality of the currents.

Figure 4: Same as in figure 2, but for the operators involving fermionic DM.

Figure 5: Same as in figure 2, but for the operators involving one vector DM.

Figure 6: Same as in figure 2, but for the operators involving two vector DM
fields. For the operator ${\mathcal{Q}}_{dHX^{2}}$, we have considered four
different cases for the masses of the two DM particles, as indicated in the
last four plots.

For the vector DM, there are 8 operators (${\mathcal{Q}}_{dX^{2}}$,
${\mathcal{Q}}_{qXX}$, ${\mathcal{Q}}_{dXX}$, ${\mathcal{Q}}_{q\tilde{X}X}$,
${\mathcal{Q}}_{d\tilde{X}X}$, ${\mathcal{Q}}_{DqX^{2}}$,
${\mathcal{Q}}_{DdX^{2}}$ and ${\mathcal{Q}}_{dHX^{2}}$, cf. eq. (4.1))
contributing to the three-body $d_{i}\to d_{j}XX$ and 6 operators
(${\mathcal{Q}}_{dHX}$, ${\mathcal{Q}}_{dX}$, ${\mathcal{Q}}_{qX}$,
${\mathcal{Q}}_{HdX}$ and ${\mathcal{Q}}_{HqX}^{(1,3)}$, cf. eqs. (46) and
(4.1)) to the two-body $d_{i}\to d_{j}X$ decay. For the Wilson coefficient of
each operator, the parameter region required to account for the Belle II
measurement of the $B^{+}\to K^{+}+\text{inv}$ decay as well as the
experimental bounds from other processes are shown in figures 5 and 6. The
following observations are made:

  * •

For the operators involving one DM field, the high-mass region of $4.3\lesssim
m_{X}\lesssim 4.6{\,\rm GeV}$ does not suffer from any constraint, while large
part of the parameter region with $m_{X}\lesssim 4.3{\,\rm GeV}$ is excluded
by the $B^{0}\to K^{*0}+\text{inv}$ decay. Specifically, for the operator
${\mathcal{Q}}_{dHX}$, the unique one built from the field strength tensor
$X_{\mu\nu}$, the parameter region with $m_{X}\lesssim 4.3{\,\rm GeV}$ is
excluded; for the other operators, parts of the parameter region still survive
for $m_{X}\lesssim 2.5{\,\rm GeV}$ or $4.3\lesssim m_{X}<4.6{\,\rm GeV}$.

  * •

For the operators involving two DM fields, except for
${\mathcal{Q}}_{dHX^{2}}$, the allowed parameter regions with $m_{X}\lesssim
1.7{\,\rm GeV}$ suffer almost no constraints from any other processes. In the
case of a heavier DM, the $B_{s}\to\text{inv}$ decay excludes a large part of
the parameter space. Especially, the parameter spaces with $m_{X}\gtrsim
2{\,\rm GeV}$ are excluded by the $B_{s}\to\text{inv}$ decay for all these
operators.

  * •

Constraints on the operator ${\mathcal{Q}}_{dHX^{2}}$ are quite different from
the other ones. As the tensor hadronic matrix element $\langle
0|\bar{b}\sigma_{\mu\nu}s|B_{s}\rangle$ vanishes identically, this operator
does not contribute to the $B_{s}\to\text{inv}$ decay. Therefore, the high-
mass regions do not suffer from any constraint. However, most of the low-mass
regions are excluded by the $B^{0}\to K^{*0}+\text{inv}$ decay, e.g., the one
with $m_{X_{1}}\lesssim 1.8{\,\rm GeV}$ is excluded in the case of
$m_{X_{2}}-m_{X_{1}}=200{\,\rm MeV}$.

All the data on $b\to s$ decays discussed above will be significantly improved
in the future Belle-II:2022cgf ; Li:2022tov ; Belle-II:2018jsg ; Amhis:2023mpj
. Considering the expected sensitivities at Belle II ($5{\,\rm ab}^{-1}$),
CEPC and FCC-ee (Tera-$Z$ phase) listed in table 2, we also show in figures
2-6 the projected $90\%$ CL bounds on the DSMEFT operators. It is observed
that

  * •

At the Belle II with $5{\,\rm ab}^{-1}$, the relative uncertainty of the
measured ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ is expected to be $19\%$.
Using the SM prediction to estimate the future experimental central value, the
projected bound covers all the parameter space explaining the recent Belle II
excess. If the central value of the current Belle II data remains unchanged in
the future, the future experimental measurement of ${\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})$ is expected to be $(23\pm 4)\times 10^{-6}$, which would
deviate from the current SM prediction by about $4\sigma$.

  * •

For the $B^{0}\to K^{0}+\text{inv}$ decay, the projected bound can cover most
of the parameter regions required to explain the $B^{+}\to K^{+}+\text{inv}$
excess. Therefore, this decay can provide an independent probe of the proposed
mechanism behind the Belle II excess.

  * •

The $B_{s}\to\phi+\text{inv}$ and $B_{s}\to\text{inv}$ decays can strongly
constrain the low- and high-mass regions of the parameter space used to
explain the Belle II excess, respectively. Their combination can even exclude
almost all the parameter regions.

  * •

All the operators involving one DM (such as ${\mathcal{Q}}_{d\phi}$ and
${\mathcal{Q}}_{da}$), as well as the operator ${\mathcal{Q}}_{dHX^{2}}$, do
not contribute to the $B_{s}\to\text{inv}$ decay. Therefore, the high-mass
regions do not suffer from any constraint, even at the Belle II with $5{\,\rm
ab}^{-1}$ or at the CEPC and FCC-ee during the Tera-$Z$ phase.

####  5.2.2 MFV

As demonstrated in eq. (50), once the MFV hypothesis is assumed, 8 out of the
22 DSMEFT operators do not induce down-type FCNC interactions. Therefore,
these operators cannot account for the recent Belle II measurement of the
$B^{+}\to K^{+}+\text{inv}$ decay. For the remaining 14 operators, their MFV
couplings make the $b\to d$ and $s\to d$ processes also relevant. We consider
the constraints from the experimental measurements summarized in table 2,
which include the branching ratios of the $B\to K^{(*)},\pi,\rho+\text{inv}$,
$K^{+}\to\pi^{+}+\text{inv}$, $K_{L}\to\pi^{0}+\text{inv}$ and
$B_{s,d}\to\text{inv}$ decays. It is noted that, as discussed in the last
subsection, the $b\to s$ transition in the MFV hypothesis approximately
depends only on one Wilson coefficient, which is similar to the treatment in
the general DSMEFT discussed in the last section; see footnote 5 for an
example. Consequently, the relative strength of the constraints from the
various $b\to s$ processes are the same between the general case and the MFV
hypothesis. Therefore, in the following, we focus on whether the $s\to d$ and
$b\to d$ processes can provide further constraint with respect to that from
the $b\to s$ processes.

Figure 7: Constraints on the MFV parameter $\epsilon_{1}$ for the operators
involving scalar DM as a function of the scalar mass $m_{\phi}$. The blue band
denotes the $2\sigma$ allowed region by the recent Belle II measurement of the
$B^{+}\to K^{+}+\text{inv}$ decay, while the upper bounds obtained from all
the other processes are shown by the solid and dashed lines. The red points
denote the benchmark points used to investigate the $q^{2}$ distributions in
figure 12; see text for more details.

For the three operators with scalar DM (cf. eq. (4.1)), the allowed regions of
the MFV parameter $\epsilon_{1}^{i}$ are shown in figure 7. In the low-mass
region of $m_{\phi}\lesssim 200{\,\rm MeV}$, the $K^{+}\to\pi^{+}+\text{inv}$
decay puts the strongest constraint for the operators ${\mathcal{Q}}_{d\phi}$
and ${\mathcal{Q}}_{d\phi^{2}}$. For all the operators, the
$B^{+}\to\pi^{+}+\text{inv}$ decay provides the strongest constraint among all
the $b\to d$ processes. However, this decay provides no further constraint
after considering the constraints from the $b\to s$ processes. The only
exception is the constraint for the operator ${\mathcal{Q}}_{\phi q}$ in the
case of $m_{\phi_{1}}=m_{\phi_{2}}\equiv m_{\phi}$. In the region of
$m_{\phi}\gtrsim 2.3{\,\rm GeV}$, the strongest constraint is provided by the
$B^{+}\to\pi^{+}+\text{inv}$ decay (as well as the $\Upsilon(1S)\to\text{inv}$
decay; see the later discussion), because $m_{\phi_{1}}=m_{\phi_{2}}$ makes
the NP contribution to the $B_{s}\to\phi_{1}\phi_{2}$ decay vanishing. In
addition, ${\mathcal{Q}}_{qa}$ is the unique operator involving ALP in the MFV
framework. The constraints on this operator are shown in figure 8, which are,
as in the general case, very similar to the ones for ${\mathcal{Q}}_{d\phi}$
but with a much larger allowed Wilson coefficient.

Figure 8: Same as in figure 7, but for the operator involving ALP.

Among the operators with fermionic DM, there is only one relevant operator
(${\mathcal{Q}}_{q\chi}$, cf. eq. (45)) in the MFV hypothesis. The
experimental constraints are shown in figure 9. We can see that, after
considering the constraints from the $b\to s$ processes, all the $s\to d$ and
$b\to d$ decays provide no further constraints, except in a narrow window near
$m_{\chi}\approx 2.3{\,\rm GeV}$, where the $B^{+}\to\pi^{+}+\text{inv}$ decay
gives the strongest bound.

Figure 9: Same as in figure 7, but for the operator involving fermionic DM.

For the vector DM, in the MFV framework, there are five operators
(${\mathcal{Q}}_{qXX}$, ${\mathcal{Q}}_{q\tilde{X}X}$,
${\mathcal{Q}}_{dX^{2}}$, ${\mathcal{Q}}_{DqX^{2}}$ and
${\mathcal{Q}}_{dHX^{2}}$, cf. eq. (4.1)) contributing to the three-body
$d_{i}\to d_{j}XX$ and four operators (${\mathcal{Q}}_{dHX}$,
${\mathcal{Q}}_{qX}$ and ${\mathcal{Q}}_{HqX}^{(1,3)}$, cf. eqs. (46) and
(4.1)) to the two-body $d_{i}\to d_{j}X$ decay. The relevant experimental
constraints are shown in figures 10 and 11. For the operators
${\mathcal{Q}}_{qX}$ and ${\mathcal{Q}}_{HqX}^{(1,3)}$, the
$K^{+}\to\pi^{+}+\text{inv}$ decay provides the strongest constraint for
$m_{X}\lesssim 0.35{\,\rm GeV}$ and exclude the parameter region required to
account for the Belle II excess. For the four operators involving one DM, the
$B^{+}\to\pi^{+}+\text{inv}$ decay excludes the mass window around
$m_{X}\approx 4.8{\,\rm GeV}$. For the four operators involving two DM,
neither the $s\to d$ nor the $b\to d$ processes can put further constraints
after considering those from the $b\to s$ decays. One exception is the
operator ${\mathcal{Q}}_{dHX^{2}}$, where the decay
$B^{+}\to\pi^{+}+\text{inv}$ excludes the tail in the high-mass region.

Figure 10: Same as in figure 7, but for the operators involving one vector DM.

Figure 11: Same as in figure 7, but for the operators involving two vector DM
fields. For the operator ${\mathcal{Q}}_{dHX^{2}}$, we have considered four
different cases for the DM masses, as indicated in the last four plots.

In the MFV hypothesis, as indicated by eqs. (17) and (18), the flavour-
conserving couplings to the third generation are much larger than all the
other FCNC couplings.777As can be seen from eq. (50), the flavour-conserving
couplings can also be altered by the MFV parameter $\epsilon_{0}$. However,
this parameter does not affect the FCNC processes, e.g., the $B^{+}\to
K^{+}+\text{inv}$ decay. For simplicity, we take $\epsilon_{0}=0$ to derive
the bounds on $\epsilon_{1}$ from the invisible decays of charmoniums and
bottomoniums. For a non-zero $\epsilon_{0}$, it is straightforward to obtain
the corresponding bounds on $\epsilon_{1}$ by using eq. (50). Therefore, the
$\Upsilon\to{\rm DM}+{\rm DM}$ decay could be highly enhanced. However, the
current experimental bounds on the invisible $\Upsilon$ decays are still quite
weak. For the $\Upsilon(1S)\to\text{inv}$ decay, the $90\%$ CL upper limit on
its branching ratio is $3.0\times 10^{-4}$ BaBar:2009gco , which is about one
order of magnitude higher than the SM prediction $9.9\times 10^{-6}$
Chang:1997tq ; Fernandez:2014eja ; Bertuzzo:2017lwt . Constraints from this
process are also shown in figures 7-11. For the operator ${\mathcal{Q}}_{\phi
q}$ with $m_{\phi_{2}}-m_{\phi_{1}}\leq 200{\,\rm MeV}$, we can see that the
constraints from the $\Upsilon(1S)\to\text{inv}$ decay are much stronger than
from all the $b\to d$ decays in the high mass region of $m_{\phi_{2}}\gtrsim
2{\,\rm GeV}$ and, especially, exclude the regions of $m_{\phi_{2}}\gtrsim
2.3{\,\rm GeV}$ required to explain the Belle II excess. However, the
constraints for all the other operators are much weaker than from the
$B$-meson and kaon decays in the range of $m_{\rm DM}\lesssim m_{B}/2$. On the
other hand, the constraints from invisible charmonium decays are totally
negligible, because the NP contributions are suppressed by a factor of
${\mathcal{O}}(\lambda_{c}^{4})$ in the MFV hypothesis.

The $q^{2}$ distributions of the FCNC decays can also provide useful
information about the NP effects. Considering the large uncertainty of the
Belle II measurement, we take a conservative estimation of the NP
contributions to the $q^{2}$ distributions: For each DSMEFT operator, the
minimal value of the MFV parameter allowed by all the experimental constraints
is chosen as the benchmark point, which is also shown in figures 7-11.
Furthermore, we choose two typical values of the DM masses, $700{\,\rm MeV}$
and $1500{\,\rm MeV}$. By using these benchmark points, we show in figure 12
our predictions for the $q^{2}$ distributions of $B^{+}\to K^{+}+\text{inv}$
and $B^{0}\to K^{*0}+\text{inv}$ decays.888The $q^{2}$ distributions of other
$P\to P+\text{inv}$ and $P\to V+\text{inv}$ decays are similar to that of
$B^{+}\to K^{+}+\text{inv}$ and $B^{0}\to K^{*0}+\text{inv}$ respectively, and
hence are not shown here anymore. In addition, for comparison, we show the
$q^{2}$ distributions of these decays by considering the benchmark values of
the Wilson coefficients of the operators ${\mathcal{Q}}_{dHX}$ and
${\mathcal{Q}}_{dHX^{2}}$, although they are already excluded by the $B^{0}\to
K^{*0}+\text{inv}$ bound. In the case of $m_{\rm DM}=700{\,\rm MeV}$, we can
see that, for the operators ${\mathcal{Q}}_{\phi q}$ and
${\mathcal{Q}}_{d\phi^{2}}$, the distributions are close to each other in the
$B^{+}\to K^{+}+\text{inv}$ decay, but distinguishable in the $B^{0}\to
K^{*0}+\text{inv}$. However, the distributions for the operators
${\mathcal{Q}}_{qXX}$ and ${\mathcal{Q}}_{q\tilde{X}X}$ are quite close in
both of these two decays. In the case of $m_{\rm DM}=1500{\,\rm MeV}$, the
distributions are quite different for all the operators, and hence
distinguishable from each other.

Figure 12: Differential branching ratios of the $B^{+}\to K^{+}+\text{inv}$
and $B^{0}\to K^{*0}+\text{inv}$ decays both within the SM and in the DSMEFT
with MFV hypothesis. For each effective operator, the DM mass is chosen either
as $700{\,\rm MeV}$ or $1500{\,\rm MeV}$, and the MFV parameters are fixed at
the corresponding benchmark points shown in figures 7–11. For the operators
involving one DM field, $m_{\rm DM}=1.7$, $2.0$, $2.2$ and $2.4{\,\rm GeV}$
are chosen and the corresponding averaged branching ratios in the $[2,\,3]$,
$[3,\,4]$, $[4,\,5]$ and $[5,\,6]{\,\rm GeV}^{2}$ bins are shown,
respectively.

As can be seen from ref. Belle-II:2023esi , excesses of the $B^{+}\to
K^{+}+\text{inv}$ events mostly appear in the bins $2\leq q^{2}\leq 7{\,\rm
GeV}^{2}$, although the uncertainties are still quite large. From the
predictions shown in figure 12, one can see that the $q^{2}$ distribution
resulting from the operator ${\mathcal{Q}}_{q\chi}$ with $m_{\chi}=700{\,\rm
MeV}$ can closely match that of the Belle II excess. As another possibility,
the two-body decay with a DM mass around $2{\,\rm GeV}$ could also contribute
to the observed distribution. In this case, the DSMEFT operators involving one
DM field become relevant, and there are 6 such operators in the MFV
hypothesis, i.e., ${\mathcal{Q}}_{d\phi}$, ${\mathcal{Q}}_{qa}$,
${\mathcal{Q}}_{dHX}$, ${\mathcal{Q}}_{qX}$ and ${\mathcal{Q}}_{HqX}^{(1,3)}$.
In order to study the possible excesses in the $[2,\,3]$, $[3,\,4]$, $[4,\,5]$
and $[5,\,6]{\,\rm GeV}^{2}$ bins, we choose $m_{\rm DM}=1.7$, $2.0$, $2.2$
and $2.4{\,\rm GeV}$ as benchmarks, respectively. For the MFV couplings of
these 6 operators, we choose the minimal values allowed by all the
experimental constraints, as done for the operators involving two DM fields.
By using these benchmark points, we also show in figure 12 the bin-averaged
branching ratios of the $B^{+}\to K^{+}+\text{inv}$ and $B^{0}\to
K^{*0}+\text{inv}$ decays. Since the minimal values allowed by the Belle II
measurement are chosen, the distributions resulting from these operators are
the same for the $B^{+}\to K^{+}+\text{inv}$ decay spectrum. The resulting
$B^{0}\to K^{*0}+\text{inv}$ decay spectra are, however, distinguishable from
each other for the operators ${\mathcal{Q}}_{dHX}$, ${\mathcal{Q}}_{qa}$ and
${\mathcal{Q}}_{d\phi}$. They can also be distinguished from that resulting
from the operators ${\mathcal{Q}}_{qX}$ and ${\mathcal{Q}}_{HqX}^{(1,3)}$,
which have the same distribution. As the two-body decay kinematics is quite
different from that of the three-body decay, a dedicated statistical analysis
can provide more insightful information, for which we refer to ref.
Altmannshofer:2023hkn  for a recent study.

From figures 7-11, we can see that the benchmark points actually correspond to
the minimal allowed parameter values required to explain the $B^{+}\to
K^{+}+\text{inv}$ excess. They are irrelevant to the $b\to d+\text{inv}$ and
$s\to d+\text{inv}$ decays. Therefore, the predicted distributions of the
$B^{+}\to K^{+}+\text{inv}$ and $B^{0}\to K^{*0}+\text{inv}$ decays in figure
12 can also be applied to the general DSMEFT case discussed in the last
subsection.

##  6 Conclusion

The recent Belle II measurement of ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$
is about $2.7\sigma$ higher than the SM prediction. In this work, we have
deciphered the data with two different NP scenarios: the underlying quark-
level $b\to s\nu\bar{\nu}$ transition is, besides the SM contribution, further
affected by some heavy new mediators that are much heavier than the
electroweak scale, or amended by an additional decay channel with light new
final states that are sufficiently long-lived to escape the detector and
appear as missing-energy signals. To make our analyses model-independent as
much as possible, we have studied these two scenarios in the SMEFT and the
DSMEFT framework, respectively. Furthermore, the flavour structures of the
resulting effective operators are taken to be either of the most generic form
or satisfy the MFV hypothesis, both for the quark and lepton sectors.

In the first scenario, we have focused on the implications of the MFV
hypothesis. It is found that, once the MFV hypothesis is assumed for the quark
sector while without any assumption on the lepton sector, only one SM-like
low-energy effective operator ${\mathcal{O}}_{L}^{\nu_{i}\nu_{j}}$ can be
induced by the SMEFT dim-6 operators. As a consequence, the ratio of the
branching fractions of any two $b\to s\nu\bar{\nu}$ decays remains a constant,
which can be used to test the MFV hypothesis. In particular, the resulting
ratio ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})_{\rm
MFV}/{\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})_{\rm MFV}=0.46\pm 0.07$ is
inconsistent with the current $90\%$ CL lower limit $0.70$ derived by combing
the recent Belle II measurement of ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$
and the Belle upper bound on ${\mathcal{B}}(B^{0}\to K^{*0}\nu\bar{\nu})$.
Thus, the Belle II excess cannot be explained in the SMEFT with MFV hypothesis
for the quark sector. Such a conclusion is true irrespective of the flavour
structure of the lepton sector. As the uncertainty of the Belle II measurement
of ${\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})$ is still quite large, we have
also used all the relevant processes but without the $B^{+}\to
K^{+}\nu\bar{\nu}$ decay to constrain the parameter space of the MFV
parameters, and then predicted that ${\mathcal{B}}(B^{+}\to
K^{+}\nu\bar{\nu})_{\rm MFV}<10.5\times 10^{-6}$ and
$3.5<{\mathcal{B}}(B^{+}\to K^{+}\nu\bar{\nu})_{\rm MFV}<5.0\times 10^{-6}$,
which correspond to the cases without and with the leptonic MFV, respectively.
These predictions can be further tested at the Belle II with more statistics.

In the second scenario, it is found that the Belle II excess can be
accommodated by 22 DSMEFT operators with dim $\leq 6$, which include 4
operators involving scalar, 2 operators involving fermionic, 14 operators
involving vector DM, as well as 2 operators involving ALP. For each operator,
the current experimental constraints on the corresponding Wilson coefficient
and DM (ALP) mass are derived. For most of these operators, large parts of the
parameter spaces required to account for the $B^{+}\to K^{+}+\text{inv}$
excess are excluded by the $B^{0}\to K^{*0}+\text{inv}$ and
$B_{s}\to\text{inv}$ decays. Nevertheless, none of the operators is fully
excluded by the current data. The future Belle II ($5{\,\rm ab}^{-1}$), CEPC
and FCC-ee (Tera-$Z$ phase) data on $B_{s}\to\phi+\text{inv}$ and
$B_{s}\to\text{inv}$ decays are, however, expected to cover almost all the
parameter spaces needed for explaining the Belle II excess. Furthermore, the
$B^{0}\to K^{0}+\text{inv}$ decay can provide an independent probe of the NP
mechanism proposed for the Belle II excess. Once the MFV hypothesis is
assumed, 8 out of the 22 viable DSMEFT operators do not induce the down-type
FCNC interactions, and hence cannot account for the Belle II excess. For the
remaining 14 operators, their MFV couplings make the $b\to d$ and $s\to d$
processes also relevant. For ${\mathcal{Q}}_{d\phi}$, ${\mathcal{Q}}_{qa}$,
${\mathcal{Q}}_{qX}$ and ${\mathcal{Q}}_{HqX}^{(1,3)}$, which all involve one
DM field, the $K^{+}\to\pi^{+}+\text{inv}$ decay provides the strongest
constraint in the low-mass region $m_{\rm DM}\lesssim m_{K}-m_{\pi}$. For
${\mathcal{Q}}_{q\chi}$, ${\mathcal{Q}}_{dHX}$, ${\mathcal{Q}}_{qX}$,
${\mathcal{Q}}_{HqX}^{(1,3)}$ and ${\mathcal{Q}}_{dHX^{2}}$, especially for
the ones (e.g., ${\mathcal{Q}}_{dHX^{2}}$) that make no contribution to the
$B_{s}\to\text{inv}$ decay, the $B^{+}\to\pi^{+}+\text{inv}$ decay puts the
most stringent bound in the high-mass region. In addition, the
$\Upsilon(1S)\to\text{inv}$ decay provides the strongest constraint in the
high-mass region for the operator ${\mathcal{Q}}_{\phi q}$ with
$m_{\phi_{2}}-m_{\phi_{1}}<200{\,\rm MeV}$. Within the parameter space allowed
by all the current experimental data, the $q^{2}$ distributions of the $B\to
K^{(*)}+\text{inv}$ decays are then studied for each viable operator. We find
that the resulting prediction of the operator
$\mathcal{Q}_{q\chi}=\big{(}\bar{q}_{p}\gamma_{\mu}q_{r}\big{)}\big{(}\bar{\chi}\gamma^{\mu}\chi\big{)}$
with a fermionic dark matter mass $m_{\chi}\approx 700{\,\rm MeV}$ can closely
match the Belle II event distribution in the bins $2\leq q^{2}\leq 7{\,\rm
GeV}^{2}$.

At the Belle II with $5{\,\rm ab}^{-1}$ dataset, the branching ratio
${\mathcal{B}}(B^{+}\to K^{+}+\text{inv})$ is expected to be measured with an
uncertainty of $19\%$. If the central value of the current Belle II
measurement remains unchanged but the uncertainty is reduced accordingly in
the future, the discrepancy from the SM prediction will increase to $4\sigma$
level. At the same time, in the future, the $B^{0}\to K^{0}+\text{inv}$,
$B^{0}\to K^{*0}+\text{inv}$, $B^{+}\to\pi^{+}+\text{inv}$ and
$B_{s}\to\text{inv}$ decays at Belle II Belle-II:2018jsg ; Belle-II:2022cgf ,
the $B_{s}\to\phi+\text{inv}$ at CEPC Li:2022tov  and FCC-ee Amhis:2023mpj  as
well as the $K\to\pi+\text{inv}$ at NA62 NA62:2017rwk  and KOTO
Yamanaka:2012yma , are all expected to provide complementary information about
the (D)SMEFT operators as well as the validity of the MFV hypothesis studied
here.

The new singlets introduced in the DSMEFT could be the DM candidates in the
universe. In this case, the parameter space derived in our work could be
further constrained by cosmological observations, such as the DM relic density
and the Big Bang nucleosynthesis, as well as the DM direct detection
Nollett:2013pwa ; Cheung:2012gi ; Balazs:2014rsa ; Liem:2016xpm ;
GAMBIT:2021rlp . For the operators involving one new particle, their couplings
should also suffer from additional constraints to ensure that the particle is
sufficiently long-lived to escape the detector. In addition, the DM FCNC
couplings can affect the neutral-meson mixings, such as $B_{s}-\bar{B}_{s}$
mixing MartinCamalich:2020dfe ; Bauer:2021mvw . At the LHC, the large
couplings of the new singlets to the third generation in MFV can lead to the
DM production associated with top and bottom quarks Lin:2013sca . However, a
complete and dedicated study of all these constraints should go beyond the EFT
framework and specific UV completions are needed. Therefore, we leave these
explorations for our future work.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
under Grant Nos. 12135006, 12075097 and 11805077, as well as by the
Fundamental Research Funds for the Central Universities under Grant No.
CCNU22LJ004. XY is also supported in part by the Startup Research Funding from
CCNU.

##  Appendix A Hadronic matrix elements

For the $B$-meson and kaon decays discussed in sections 3 and 4, as well as
the invisible decays of charmoniums and bottomoniums, we need the hadronic
matrix elements between the vacuum and the initial hadronic states or between
the initial and the final hadronic state, which can be parameterized by the
decay constants or by the transition form factors. For convenience, they are
recapitulated in appendices A.1 and A.2, respectively. The matrix elements of
other bilinear quark currents not mentioned explicitly are zero due to parity
invariance of strong interaction.

###  A.1 Decay constants

For the $B_{q}$ meson ($q=u,d$), the non-vanishing annihilation matrix
elements are defined as

$\displaystyle\langle
0|\bar{q}\gamma_{\mu}\gamma_{5}b|\bar{B}_{q}(p)\rangle=if_{B_{q}}\,p_{\mu},\quad\quad\langle
0|\bar{q}\gamma_{5}b|\bar{B}_{q}(p)\rangle=-i\frac{m_{B_{q}}^{2}}{\overline{m}_{q}(\mu)+\overline{m}_{b}(\mu)}f_{B_{q}},$
(72)

with $\overline{m}_{b}(\mu)$ and $\overline{m}_{q}(\mu)$ denoting the $b$\-
and $q$-quark running masses at the scale $\mu$ defined in the
$\overline{\text{MS}}$ scheme. Here $f_{B_{q}}$ denote the $B_{q}$-meson decay
constants, and their up-to-date Lattice QCD values can be found in ref.
FlavourLatticeAveragingGroupFLAG:2021npn  and references therein. For the
vector quarkonium $V$ (such as $\Upsilon$ and $J/\psi$), the relevant hadronic
matrix elements are given by Aloni:2017eny ; Hatton:2021dvg

$\displaystyle\langle 0|\bar{q}\gamma^{\mu}q|V(p,\varepsilon)\rangle$
$\displaystyle=f_{V}m_{V}\varepsilon^{\mu},$ $\displaystyle\langle
0|\bar{q}\sigma^{\mu\nu}q|V(p,\varepsilon)\rangle$
$\displaystyle=if_{V}^{T}(\mu)\big{(}\varepsilon^{\mu}p^{\nu}-\varepsilon^{\nu}p^{\mu}\big{)},$
(73)

where $\varepsilon^{\mu}$ denotes the polarization vector of the vector
quarkonium, and $f_{V}^{(T)}$ is the decay constant. The matrix element of the
axial-tensor current can be related to that of the tensor current, while the
ones of all the other currents are zero. In the numerical analysis, we use the
values of $f_{\Upsilon}$ in ref. Hatton:2021dvg , $f_{J/\psi}$ in ref.
Hatton:2020qhk , and $f_{V}^{T}(2{\,\rm GeV})/f_{V}$ in ref. Konig:2015qat .

###  A.2 Transition form factors

The transition form factors for a $B$ meson decaying into a pseudoscalar meson
$P$ (e.g., $B\to K$) are defined by Beneke:2000wa ; Ball:2004ye

$\displaystyle c_{P}\langle
P(p^{\prime})|\bar{q}\gamma^{\mu}b|\bar{B}(p)\rangle$
$\displaystyle=f_{+}(q^{2})\left(P^{\mu}-\frac{m_{B}^{2}-m_{P}^{2}}{q^{2}}q^{\mu}\right)+f_{0}(q^{2})\frac{m_{B}^{2}-m_{P}^{2}}{q^{2}}q^{\mu},$
$\displaystyle c_{P}\langle P(p^{\prime})|\bar{q}b|\bar{B}(p)\rangle$
$\displaystyle=\frac{m_{B}^{2}-m_{P}^{2}}{\overline{m}_{b}(\mu)-\overline{m}_{q}(\mu)}f_{0}(q^{2}),$
$\displaystyle c_{P}\langle
P(p^{\prime})|\bar{q}\sigma^{\mu\nu}b(p)|\bar{B}(p)\rangle$
$\displaystyle=i\frac{P^{\mu}q^{\nu}-P^{\nu}q^{\mu}}{m_{B}+m_{P}}f_{T}(q^{2}),$
(74)

where $q=p-p^{\prime}$, $P=p+p^{\prime}$, and $c_{P}=-\sqrt{2}$ for $\pi^{0}$
and $c_{P}=1$ otherwise. For a $B$ meson decaying into a vector meson $V$, the
relevant transition form factors are introduced by the following Lorentz
decompositions of the bilinear quark current matrix elements Beneke:2000wa ;
Ball:2004rg :

$\displaystyle c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\gamma_{\mu}b|\bar{B}(p)\rangle$
$\displaystyle=-\epsilon_{\mu\nu\rho\sigma}\varepsilon^{*\nu}P^{\rho}q^{\sigma}\frac{V(q^{2})}{m_{B}+m_{V}},$
$\displaystyle c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\gamma^{\mu}\gamma^{5}b|\bar{B}(p)\rangle$
$\displaystyle=i\varepsilon^{*\mu}(m_{B}+m_{V})A_{1}(q^{2})-iP^{\mu}\frac{\varepsilon^{*}\cdot
q}{m_{B}+m_{V}}A_{2}(q^{2})$ $\displaystyle\qquad-
iq^{\mu}(\varepsilon^{*}\cdot
q)\frac{2m_{V}}{q^{2}}\left[A_{3}(q^{2})-A_{0}(q^{2})\right],$ $\displaystyle
c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\gamma^{5}b|\bar{B}(p)\rangle$
$\displaystyle=-2i\frac{m_{V}}{\overline{m}_{b}(\mu)+\overline{m}_{q}(\mu)}(\varepsilon^{*}\cdot
q)A_{0}(q^{2}),$ $\displaystyle c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\sigma_{\mu\nu}q^{\nu}b|\bar{B}(p)\rangle$
$\displaystyle=-i\epsilon_{\mu\nu\rho\sigma}\varepsilon^{*\nu}P^{\rho}q^{\sigma}T_{1}(q^{2}),$
$\displaystyle c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\sigma_{\mu\nu}q^{\nu}\gamma^{5}b|\bar{B}(p)\rangle$
$\displaystyle=\left[(m_{B}^{2}-m_{V}^{2})\varepsilon^{*}_{\mu}-(\varepsilon^{*}\cdot
q)P_{\mu}\right]T_{2}(q^{2})$ $\displaystyle\qquad+(\varepsilon^{*}\cdot
q)\left(q_{\mu}-\frac{q^{2}}{m_{B}^{2}-m_{V}^{2}}P_{\mu}\right)T_{3}(q^{2}),$
(75)

where $\varepsilon$ denotes the polarization vector of the vector meson, and
$c_{V}=-\sqrt{2}$ for $\rho^{0}$ and $c_{V}=1$ otherwise. We adopt the
convention $\epsilon^{0123}=+1$ and
$\sigma_{\mu\nu}=i/2\left[\gamma_{\mu},\gamma_{\nu}\right]$. With the above
parameterizations, the absence of dynamical singularities at $q^{2}=0$ implies
that $f_{+}(0)=f_{0}(0)$ and $A_{0}(0)=A_{3}(0)$. In addition, the algebraic
relation between $\sigma_{\mu\nu}$ and $\sigma_{\mu\nu}\gamma_{5}$,
$\sigma_{\mu\nu}\gamma_{5}=i/2\epsilon_{\mu\nu\alpha\beta}\sigma^{\alpha\beta}$,
gives rise to the identity $T_{1}(0)=T_{2}(0)$. It is also noted that the form
factor $A_{3}$ is redundant and can be written as a linear combination of
$A_{1}$ and $A_{2}$ Ball:2004rg

$\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}).$
(76)

The hadronic matrix elements of the tensor and axial-tensor currents can also
be decomposed as Khodjamirian:2020btr

$\displaystyle c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\sigma^{\mu\nu}b|\bar{B}(p)\rangle=$
$\displaystyle
i\epsilon_{\mu\nu\rho\sigma}\left[P^{\rho}\varepsilon^{*\sigma}-\frac{m_{B}^{2}-m_{V}^{2}}{q^{2}}q^{\rho}\varepsilon^{*\sigma}+\frac{\varepsilon^{*}\cdot
q}{q^{2}}q^{\rho}P^{\sigma}\right]T_{1}(q^{2})$
$\displaystyle+i\epsilon_{\mu\nu\rho\sigma}\left[\frac{m_{B}^{2}-m_{V}^{2}}{q^{2}}q^{\rho}\varepsilon^{*\sigma}+\frac{\varepsilon^{*}\cdot
q}{q^{2}}P^{\rho}q^{\sigma}\right]T_{2}(q^{2})$
$\displaystyle+i\epsilon_{\mu\nu\rho\sigma}\frac{\varepsilon^{*}\cdot
q}{m_{B}^{2}-m_{V}^{2}}P^{\rho}q^{\sigma}T_{3}(q^{2}),$ $\displaystyle
c_{V}\langle
V(p^{\prime},\varepsilon^{*})|\bar{q}\sigma^{\mu\nu}\gamma^{5}b|\bar{B}(p)\rangle=$
$\displaystyle\bigg{[}(\varepsilon^{*\mu}P^{\nu}-P^{\mu}\varepsilon^{*\nu})+\frac{m_{B}^{2}-m_{V}^{2}}{q^{2}}(q^{\mu}\varepsilon^{*\nu}-\varepsilon^{*\mu}q^{\nu})$
$\displaystyle+\frac{\varepsilon^{*}\cdot
q}{q^{2}}(P^{\mu}q^{\nu}-q^{\mu}P^{\nu})\bigg{]}T_{1}(q^{2})$
$\displaystyle+\left[\frac{m_{B}^{2}-m_{V}^{2}}{q^{2}}(\varepsilon^{*\mu}q^{\nu}-q^{\mu}\varepsilon^{*\nu})+\frac{\varepsilon^{*}\cdot
q}{q^{2}}(q^{\mu}P^{\nu}-P^{\mu}q^{\nu})\right]T_{2}(q^{2})$
$\displaystyle+\frac{\varepsilon^{*}\cdot
q}{m_{B}^{2}-m_{V}^{2}}\left(q^{\mu}P^{\nu}-P^{\mu}q^{\nu}\right)T_{3}(q^{2}),$
(77)

from which the forms given in eq. (A.2) can be obtained by multiplying both
sides by $q_{\nu}$. All the $B_{s}\to P\,(V)$ and $K\to\pi$ form factors can
be defined similarly.

Numerically, we use the values obtained in ref. Ball:2004ye  for the $B\to P$
and in ref. Ball:2004rg  for the $B\to V$ transition form factors. As the
Lattice QCD results are still not available for all of these form factors,
especially for the tensor ones (see ref.
FlavourLatticeAveragingGroupFLAG:2021npn  and references therein), we use
instead the calculations based on the light-cone sum rule approach Ball:2004ye
; Ball:2004rg . For the $K\to\pi$ vector and scalar form factors, the results
obtained in refs. Bernard:2006gy ; Bernard:2009zm  are adopted, while the
value in ref. Baum:2011rm  is used for the $K\to\pi$ tensor form factor.

For the operator ${\mathcal{Q}}_{DdX^{2}}$ defined in eq. (4.1), the relevant
form factors are not available in the literature so far. However, when the
Wilson coefficient ${\mathcal{C}}_{DdX^{2}}$ is assumed to be hermitian, they
can be derived from the above form factors. Taking the $\bar{B}\to\bar{K}^{0}$
transition as an example, we have the relation
$\langle\bar{K}^{0}(p^{\prime})\lvert(\partial^{\mu}\bar{s})\gamma^{\nu}P_{L}b+\bar{s}\gamma^{\nu}P_{L}(\partial^{\mu}b)\rvert\bar{B}(p)\rangle=-iq^{\mu}\langle\bar{K}^{0}(p^{\prime})\lvert\bar{s}\gamma^{\nu}P_{L}b\rvert\bar{B}(p)\rangle$.

## References

  * (1) S. L. Glashow, J. Iliopoulos, and L. Maiani, Weak Interactions with Lepton-Hadron Symmetry, Phys. Rev. D 2 (1970) 1285–1292.
  * (2) W. Altmannshofer, A. J. Buras, D. M. Straub, and M. Wick, New strategies for New Physics search in $B\to K^{*}\nu\bar{\nu}$, $B\to K\nu\bar{\nu}$ and $B\to X_{s}\nu\bar{\nu}$ decays, JHEP 04 (2009) 022, [arXiv:0902.0160].
  * (3) A. J. Buras, J. Girrbach-Noe, C. Niehoff, and D. M. Straub, $B\to{K}^{\left(\ast\right)}\nu\overline{\nu}$ decays in the Standard Model and beyond, JHEP 02 (2015) 184, [arXiv:1409.4557].
  * (4) D. Bečirević, G. Piazza, and O. Sumensari, Revisiting $B\rightarrow K^{(*)}\nu{\bar{\nu}}$ decays in the Standard Model and beyond, Eur. Phys. J. C 83 (2023), no. 3 252, [arXiv:2301.06990].
  * (5) T. Felkl, S. L. Li, and M. A. Schmidt, A tale of invisibility: constraints on new physics in b → s$\nu$$\nu$, JHEP 12 (2021) 118, [arXiv:2111.04327].
  * (6) R. Bause, H. Gisbert, M. Golz, and G. Hiller, Interplay of dineutrino modes with semileptonic rare B-decays, JHEP 12 (2021) 061, [arXiv:2109.01675].
  * (7) X. G. He and G. Valencia, $R_{K(\ast)}^{\nu}$ and non-standard neutrino interactions, Phys. Lett. B 821 (2021) 136607, [arXiv:2108.05033].
  * (8) T. E. Browder, N. G. Deshpande, R. Mandal, and R. Sinha, Impact of $B\to K\nu\bar{\nu}$ measurements on beyond the Standard Model theories, Phys. Rev. D 104 (2021), no. 5 053007, [arXiv:2107.01080].
  * (9) Belle-II Collaboration, I. Adachi et al., Evidence for $B^{+}\to K^{+}\nu\bar{\nu}$ Decays, arXiv:2311.14647.
  * (10) Belle-II Collaboration, W. Altmannshofer et al., The Belle II Physics Book, PTEP 2019 (2019), no. 12 123C01, [arXiv:1808.10567]. [Erratum: PTEP 2020, 029201 (2020)].
  * (11) Belle-II Collaboration, L. Aggarwal et al., Snowmass White Paper: Belle II physics reach and plans for the next decade and beyond, arXiv:2207.06307.
  * (12) P. Athron, R. Martinez, and C. Sierra, $B$ meson anomalies and large $B^{+}\to K^{+}\nu\bar{\nu}$ in non-universal $U(1)^{\prime}$ models, arXiv:2308.13426.
  * (13) R. Bause, H. Gisbert, and G. Hiller, Implications of an enhanced $B\to K\nu\bar{\nu}$ branching ratio, Phys. Rev. D 109 (2024), no. 1 015006, [arXiv:2309.00075].
  * (14) L. Allwicher, D. Becirevic, G. Piazza, S. Rosauro-Alcaraz, and O. Sumensari, Understanding the first measurement of $\mathcal{B}(B\to K\nu\bar{\nu})$, Phys. Lett. B 848 (2024) 138411, [arXiv:2309.02246].
  * (15) T. Felkl, A. Giri, R. Mohanta, and M. A. Schmidt, When energy goes missing: new physics in $b\rightarrow s\nu\nu$ with sterile neutrinos, Eur. Phys. J. C 83 (2023), no. 12 1135, [arXiv:2309.02940].
  * (16) M. Abdughani and Y. Reyimuaji, Constraining light dark matter and mediator with $B^{+}\rightarrow K^{+}\nu\bar{\nu}$ data, arXiv:2309.03706.
  * (17) H. K. Dreiner, J. Y. Günther, and Z. S. Wang, The Decay $B\to K\nu\bar{\nu}$ at Belle II and a Massless Bino in R-parity-violating Supersymmetry, arXiv:2309.03727.
  * (18) X.-G. He, X.-D. Ma, and G. Valencia, Revisiting models that enhance $B^{+}\to K^{+}\nu\bar{\nu}$ in light of the new Belle II measurement, arXiv:2309.12741.
  * (19) A. Berezhnoy and D. Melikhov, $B\to K^{*}M_{X}$ vs $B\to KM_{X}$ as a probe of a scalar-mediator dark matter scenario, EPL 145 (2024), no. 1 14001, [arXiv:2309.17191].
  * (20) A. Datta, D. Marfatia, and L. Mukherjee, $B\to K\nu\bar{\nu}$, MiniBooNE and muon g-2 anomalies from a dark sector, Phys. Rev. D 109 (2024), no. 3 L031701, [arXiv:2310.15136].
  * (21) W. Altmannshofer, A. Crivellin, H. Haigh, G. Inguglia, and J. Martin Camalich, Light New Physics in $B\to K^{(*)}\nu\bar{\nu}$?, arXiv:2311.14629.
  * (22) D. McKeen, J. N. Ng, and D. Tuckler, Higgs Portal Interpretation of the Belle II $B^{+}\to K^{+}\nu\nu$ Measurement, arXiv:2312.00982.
  * (23) K. Fridell, M. Ghosh, T. Okui, and K. Tobioka, Decoding the $B\to K\nu\nu$ excess at Belle II: kinematics, operators, and masses, arXiv:2312.12507.
  * (24) S.-Y. Ho, J. Kim, and P. Ko, Recent $B^{+}\\!\to K^{+}\nu\bar{\nu}$ Excess and Muon $g-2$ Illuminating Light Dark Sector with Higgs Portal, arXiv:2401.10112.
  * (25) F.-Z. Chen, Q. Wen, and F. Xu, Correlating $B\to K^{(\ast)}\nu\bar{\nu}$ and flavor anomalies in SMEFT, arXiv:2401.11552.
  * (26) E. Gabrielli, L. Marzola, K. Müürsepp, and M. Raidal, Explaining the $B^{+}\to K^{+}\nu\bar{\nu}$ excess via a massless dark photon, arXiv:2402.05901.
  * (27) T. Li, Z. Qian, M. A. Schmidt, and M. Yuan, The quark flavor-violating ALPs in light of B mesons and hadron colliders, arXiv:2402.14232.
  * (28) Particle Data Group Collaboration, R. L. Workman et al., Review of Particle Physics, PTEP 2022 (2022) 083C01. and 2023 update.
  * (29) G. Alonso-Álvarez and M. Escudero, The first limit on invisible decays of $B_{s}$ mesons comes from LEP, arXiv:2310.13043.
  * (30) I. Brivio and M. Trott, The Standard Model as an Effective Field Theory, Phys. Rept. 793 (2019) 1–98, [arXiv:1706.08945].
  * (31) G. Isidori, F. Wilsch, and D. Wyler, The Standard Model effective field theory at work, arXiv:2303.16922.
  * (32) G. Buchalla, A. J. Buras, and M. E. Lautenbacher, Weak decays beyond leading logarithms, Rev. Mod. Phys. 68 (1996) 1125–1144, [hep-ph/9512380].
  * (33) E. E. Jenkins, A. V. Manohar, and P. Stoffer, Low-Energy Effective Field Theory below the Electroweak Scale: Operators and Matching, JHEP 03 (2018) 016, [arXiv:1709.04486].
  * (34) W. Dekens and P. Stoffer, Low-energy effective field theory below the electroweak scale: matching at one loop, JHEP 10 (2019) 197, [arXiv:1908.05295]. [Erratum: JHEP 11, 148 (2022)].
  * (35) J. Aebischer, A. Crivellin, M. Fael, and C. Greub, Matching of gauge invariant dimension-six operators for $b\to s$ and $b\to c$ transitions, JHEP 05 (2016) 037, [arXiv:1512.02830].
  * (36) E. Del Nobile and F. Sannino, Dark Matter Effective Theory, Int. J. Mod. Phys. A 27 (2012) 1250065, [arXiv:1102.3116].
  * (37) M. Baumgart et al., Snowmass White Paper: Effective Field Theories for Dark Matter Phenomenology, arXiv:2203.08204.
  * (38) M. Duch, B. Grzadkowski, and J. Wudka, Classification of effective operators for interactions between the Standard Model and dark matter, JHEP 05 (2015) 116, [arXiv:1412.0520].
  * (39) J. F. Kamenik and C. Smith, FCNC portals to the dark sector, JHEP 03 (2012) 090, [arXiv:1111.6402].
  * (40) J. Brod, A. Gootjes-Dreesbach, M. Tammaro, and J. Zupan, Effective Field Theory for Dark Matter Direct Detection up to Dimension Seven, JHEP 10 (2018) 065, [arXiv:1710.10218]. [Erratum: JHEP 07, 012 (2023)].
  * (41) J. C. Criado, A. Djouadi, M. Perez-Victoria, and J. Santiago, A complete effective field theory for dark matter, JHEP 07 (2021) 081, [arXiv:2104.14443].
  * (42) J. Aebischer, W. Altmannshofer, E. E. Jenkins, and A. V. Manohar, Dark matter effective field theory and an application to vector dark matter, JHEP 06 (2022) 086, [arXiv:2202.06968].
  * (43) H. Song, H. Sun, and J.-H. Yu, Complete EFT Operator Bases for Dark Matter and Weakly-Interacting Light Particle, arXiv:2306.05999.
  * (44) X.-G. He, X.-D. Ma, and G. Valencia, FCNC B and K meson decays with light bosonic Dark Matter, JHEP 03 (2023) 037, [arXiv:2209.05223].
  * (45) J.-H. Liang, Y. Liao, X.-D. Ma, and H.-L. Wang, Dark sector effective field theory, JHEP 12 (2023) 172, [arXiv:2309.12166].
  * (46) R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38 (1977) 1440–1443.
  * (47) R. D. Peccei and H. R. Quinn, Constraints Imposed by CP Conservation in the Presence of Instantons, Phys. Rev. D 16 (1977) 1791–1797.
  * (48) S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40 (1978) 223–226.
  * (49) F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279–282.
  * (50) H. Georgi, D. B. Kaplan, and L. Randall, Manifesting the Invisible Axion at Low-energies, Phys. Lett. B 169 (1986) 73–78.
  * (51) I. Brivio, M. B. Gavela, L. Merlo, K. Mimasu, J. M. No, R. del Rey, and V. Sanz, ALPs Effective Field Theory and Collider Signatures, Eur. Phys. J. C 77 (2017), no. 8 572, [arXiv:1701.05379].
  * (52) M. Chala, G. Guedes, M. Ramos, and J. Santiago, Running in the ALPs, Eur. Phys. J. C 81 (2021), no. 2 181, [arXiv:2012.09017].