#  FCNC $B$ and $K$ Meson Decays with Light Bosonic Dark Matter

Xiao-Gang He<EMAIL_ADDRESS>Tsung-Dao Lee Institute (TDLI) & School of
Physics and Astronomy (SPA), Shanghai Jiao Tong University (SJTU), Shanghai
200240, China  Department of Physics, National Taiwan University, Taipei 10617
Xiao-Dong Ma<EMAIL_ADDRESS>Tsung-Dao Lee Institute (TDLI) & School of
Physics and Astronomy (SPA), Shanghai Jiao Tong University (SJTU), Shanghai
200240, China  German Valencia<EMAIL_ADDRESS>School of Physics
and Astronomy, Monash University, Wellington Road, Clayton, VIC-3800,
Australia

###### Abstract

We consider decays of $B$ and $K$ mesons into a pseudo-scalar or vector meson
plus missing energy. Within the SM, these modes originate from flavor changing
neutral current (FCNC) processes with two neutrinos in the final state. In
this paper we consider the experimental upper bounds on these modes and
interpret the difference between these bounds and the SM prediction as a
window into new light invisible particles. In particular we consider the case
where some new symmetry requires the new particles to be produced in pairs. We
first construct the general low energy effective Lagrangian coupling an FCNC
with two dark sector particles of spin zero, one-half and one. We then present
numerical estimates for the constraints that can be placed on these
interactions, finding that an effective new physics scale from ${\cal
O}(10)$-${\cal O}(10^{11})$ GeV can be probed, with the exact value strongly
depending on the interaction structure as well as the mass of the invisible
particle. For $K^{+}\to\pi^{+}\not{E}$ we incorporate into our constraints the
effect of using only the signal regions of NA62, and for $B^{+}\to
K^{+}\not{E}$ the $q^{2}$-dependent efficiency of Belle II.

###### Contents

  1. I Introduction
  2. II Quark-DM interaction in LEFT
  3. III $B\to(K,K^{*},\pi,\rho)$+DM+DM
    1. III.1 Form factors
    2. III.2  $B\to M+\phi\phi$ with scalar DM $\phi$
    3. III.3 $B\to M+XX$ with vector DM $X$: scenario A
    4. III.4 $B\to M+XX$ with vector DM $X$: scenario B
  4. IV $K\to\pi$+DM+DM
  5. V Summary and Conclusions
  6. A Differential decay rate
  7. B Detailed analysis for the operators with a vector DM
  8. C Lepton-DM interaction in LEFT
  9. D Form factors for the $B\to P(V)$ transitions
  10. E Specific renormalizable models to illustrate a possible origin of the LEFT operators

##  I Introduction

It is well known that flavor changing neutral current (FCNC) processes are
severely suppressed within the standard model (SM) due to both their loop
origin and the GIM mechanism Glashow:1970gm . For this reason, these processes
are very sensitive to new physics (NP) beyond the SM. We consider $B$ and $K$
decay modes with one meson and a neutrino pair in the final state because
experimental upper bounds are available for them. As the neutrinos are not
detected, these modes become decays with missing energy ($\not{E}$) in the
final state from the experimental perspective. At the quark level, any
transition of the form $b\to s(d)+\not{E}$ and $s\to d+\not{E}$, would
contribute to these modes and $\not{E}$ can originate from any sufficiently
light invisible particle. In addition to the SM neutrinos, the invisible
particles can be any hypothetical neutral particle that escapes the current
experimental detection. One well-motivated choice is to relate these light
invisible particles to dark matter (DM) or other dark sectors. These rare
meson decays can then be used to constrain light DM and are particularly
important in the face of current stringent experimental constraints on heavy
DM from direct detection experiments Roszkowski:2017nbc ; Bottaro:2021snn .

Topology Particle | Fermion $(\chi)$ | Scalar $(\phi)$ | Vector $(X)$
---|---|---|---
2-body decay |  $b\to s(d)+\chi$ Dib:2022ppx  |  $b\to s(d)+\phi$ Kamenik:2011vy ; Li:2021sqe  |  $b\to s(d)+X$ Kamenik:2011vy ; Li:2021sqe
— |  $s\to d+\phi$ Kamenik:2011vy ; He:2020jly  |  $s\to d+X$ Kamenik:2011vy
3-body decay |  $b\to s(d)+\chi\chi$ Kamenik:2011vy ; Li:2020dpc ; Felkl:2021uxi  |  $b\to s(d)+\phi\phi$ (✓) |  $b\to s(d)+XX$ (✓)
$s\to d+\chi\chi$ Kamenik:2011vy ; Li:2019fhz ; Deppisch:2020oyx ; He:2021yoz  |  $s\to d+\phi\phi$ (✓) |  $s\to d+XX$ (✓)

Table 1: List of possible $b\to s(d)+\not{E}$ and $s\to d+\not{E}$ modes: a
checkmark that they have not been studied before and are the subject of this
paper.

Without knowledge of the fundamental interactions in dark sectors, it is
suitable to study these FCNC processes in a model independent manner by
incorporating new light, neutral, degrees of freedom with an effective field
theory (EFT) approach. The new particles can be either scalar, fermion, or
vector in nature if we limit our study to spin less than or equal to one. Well
known examples for these invisible particles are the axion Weinberg:1977ma ;
Wilczek:1977pj , the sterile neutrino Dasgupta:2021ies , and the dark photon
Fabbrichesi:2020wbt . Model independent studies of the FCNC $b\to
s(d)+\not{E}$ and $s\to d+\not{E}$ transitions with the EFT framework have
been carried out in the past for these three types of invisible particles. A
summary of what has been done is shown in Tab. 1. These processes are further
divided into 2-body decays with a single new particle and 3-body decays with a
pair of new particles. All the 2-body channels plus the 3-body decay with a
pair of new fermions have been extensively studied before Dib:2022ppx ;
Kamenik:2011vy ; Li:2021sqe ; He:2020jly ; Li:2020dpc ; Felkl:2021uxi ;
Li:2019fhz ; Deppisch:2020oyx ; He:2021yoz . 111Note that the 2-body modes
with a single new fermion are only possible when the meson decays into a
baryon and $\chi$ Dib:2022ppx . For the 3-body decay with a pair of fermions,
$\chi$ can be either the neutrino or some DM particle, the former is studied
in Felkl:2021uxi ; Li:2019fhz ; Deppisch:2020oyx ; He:2021yoz , while the
latter in Li:2020dpc .  The table indicates with a checkmark processes that
have not received much attention, namely three body modes with a pair of
scalar or vector particles. We are only aware of a classification of relevant
SMEFT operators with some basic estimates in Kamenik:2011vy . In this paper,
we systematically investigate these channels using a general low energy
effective theory (LEFT), and find constraints on all the relevant effective
operators with the help of the most recent experimental results.

The observables we use to set the bounds are summarized in Tab. 2. In the
second column, we list the SM background with a pair of neutrinos. For most
modes in the table, we obtain the SM prediction using the flavio package
Straub:2018kue . In some cases, the SM calculation of the form factors is not
trivial and we use instead the predictions for $B\to\pi\bar{\nu}\nu$ from
Hambrock:2015wka . The kaon decay modes are very clean theoretically, but
their parametric uncertainty due to CKM angles can be large, with central
values changing by up to 20% Buras:2015qea , we quote the values from the PDG
in this case Workman:2022ynf .

The new physics we discuss, will always add incoherently to the SM di-neutrino
background in $B(K)\to M\not{E}$ processes (with $M$ representing the final
state mesons shown in Tab. 2). In view of this we define the “room for new
physics” in these modes as the difference between the experimental upper bound
and the standard model prediction. We adopt the simple prescription of
subtracting from the 90 % experimental upper bound the lower limit of the
$90\,\%$ C.L. SM range and show this number in the last column. The
constraints can be easily adapted for more sophisticated subtractions if
desired. The mode $K^{+}\to\pi^{+}\bar{\nu}\nu$ has been measured by both
BNL787/949 E949:2008btt ; BNL-E949:2009dza  and NA62 NA62:2020fhy ;
NA62:2021zjw , so in this case we use the upper limit of the $90\,\%$ range
quoted by PDG Workman:2022ynf  as the experimental upper limit.

Observable |  SM prediction ($\not{E}=\nu\bar{\nu}$) |  90 % C.L. upper bound 222Except for $K^{+}\to\pi^{+}\not{E}$ channel, for which the $1\,\sigma$ measured value quoted by PDG Workman:2022ynf  is used. | New physics bound ${\cal B}^{\rm UL}$
---|---|---|---
${\cal B}(B^{+}\to K^{+}\not{E})$ | $(4.4\pm 0.6)\times 10^{-6}$ |  $1.6\times 10^{-5}$ Workman:2022ynf  | $1.3\times 10^{-5}$
${\cal B}(B^{0}\to K^{0}\not{E})$ | $(4.1\pm 0.6)\times 10^{-6}$ |  $2.6\times 10^{-5}$ Belle:2017oht  | $2.3\times 10^{-5}$
${\cal B}(B^{+}\to K^{*+}\not{E})$ | $(1.0\pm 0.1)\times 10^{-5}$ |  $4.0\times 10^{-5}$ Belle:2013tnz  | $3.1\times 10^{-5}$
${\cal B}(B^{0}\to K^{*0}\not{E})$ | $(9.5\pm 1.0)\times 10^{-6}$ |  $1.8\times 10^{-5}$ Belle:2017oht  | $1.0\times 10^{-5}$
${\cal B}(B^{+}\to\pi^{+}\not{E})$ |  $(2.39_{-0.28}^{+0.30})\times 10^{-7}$ Hambrock:2015wka  |  $1.4\times 10^{-5}$ Belle:2017oht  | $1.4\times 10^{-5}$
${\cal B}(B^{0}\to\pi^{0}\not{E})$ |  $(1.2_{-0.14}^{+0.15})\times 10^{-7}$ Hambrock:2015wka  |  $9.0\times 10^{-6}$ Belle:2017oht  | $8.9\times 10^{-6}$
${\cal B}(B^{+}\to\rho^{+}\not{E})$ | $(4.5\pm 1.0)\times 10^{-7}$ |  $3.0\times 10^{-5}$ Belle:2017oht  | $3.0\times 10^{-5}$
${\cal B}(B^{0}\to\rho^{0}\not{E})$ | $(2.0\pm 0.4)\times 10^{-7}$ |  $4.0\times 10^{-5}$ Belle:2017oht  | $4.0\times 10^{-5}$
${\cal B}(K^{+}\to\pi^{+}\not{E})$ |  $(8.1\pm 0.4)\times 10^{-11}$ Workman:2022ynf  |  $(1.14_{-0.33}^{+0.40})\times 10^{-10}$ Workman:2022ynf  | $1.1\times 10^{-10}$
${\cal B}(K_{L}\to\pi^{0}\not{E})$ |  $(2.8\pm 0.2)\times 10^{-11}$ Workman:2022ynf  |  $4.9\times 10^{-9}$ Workman:2022ynf  | $4.9\times 10^{-9}$

Table 2: Summary on the status of FCNC $B\to(K,K^{*},\pi,\rho)\not{E}$ and
$K\to\pi\not{E}$ decays with missing energy.

The paper is organized as follows. In section II we classify all the relevant
effective interactions in the framework of low energy effective field theory.
In section III we consider the FCNC $B$ meson decay
$B\to(K,K^{*},\pi,\rho)$+DM+DM with a pair of scalar or vector DM, and use the
current experimental bounds to constrain the relevant effective new physics
scale involving $(bs)$ and $(bd)$ quark flavors. In section IV, we use chiral
perturbation theory to analyze the FCNC kaon decay $K\to\pi$+DM+DM to
constrain the effective scale involving $(ds)$ quark flavors. In all cases we
will refer to the invisible light particles as “DM” regardless of their
origin. In section V, we draw our conclusions. Supplementary material
presented in the Appendix includes: the phase space integration in A; the
reduction of operators with vector DM fields in B; operators for lepton-DM
interactions in C; a collection of form factors involving $B$ meson decays in
D, and specific renormalizable model realizations for some scalar and vector
DM operators as an illustration in appendix E.

##  II Quark-DM interaction in LEFT

In this section we list the most general local quark-DM interactions in the
framework of low energy effective field theory, with the FCNC interactions
being a subset of these operators with appropriate flavor indices. In LEFT,
only the unbroken $SU(3)_{c}\times U(1)_{\rm em}$ symmetry of the standard
model is imposed. Operators at the weak scale in SMEFT have been listed before
Grzadkowski:2010es ; Lehman:2014jma ; Liao:2016hru  along with their
contributions to the LEFT operators Jenkins:2017jig ; Liao:2020zyx . Since we
are considering cases with light new particles it is more general to start
from LEFT Jenkins:2017jig . In this way we cover scenarios that contain both
light DM as well as weak scale new mediators that can be integrated out to
reach the LEFT description containing only the relevant degrees of freedom
Lehmann:2020lcv .

The scenario in which modes with two DM particles becomes relevant is that in
which they must be pair produced due to underlying symmetries in the dark
sector which can also guarantee that the DM is stable. In this case the
relevant quark-DM interactions have to involve a pair of quark fields and a
pair of DM fields. These interactions can be classified in terms of the spin
of the DM particle, and we consider the cases of spin zero, one-half and one
with corresponding scalar, fermion and vector DM fields. We denote the light
SM quarks as $q\in\\{d,s,b,u,c\\}$, the fermionic DM particle as $\chi$ (Dirac
or Majorana fermion), the scalar DM as $\phi$ (complex or real scalar), and
the vector DM as $X$ (complex or real vector), respectively. For each type of
DM particle, the relevant quark-DM interactions are enumerated below.

Fermion case: The leading dimension-6 (dim-6) operators for this case have
been considered before Kumar:2013iva ; Badin:2010uh  and we list them (with
slightly different convention) here for completeness. They are,

$\displaystyle{\cal O}_{q\chi 1}^{S}$
$\displaystyle=(\overline{q}q)(\overline{\chi}\chi),$ $\displaystyle{\cal
O}_{q\chi 2}^{S}$
$\displaystyle=(\overline{q}q)(\overline{\chi}i\gamma_{5}\chi),$ (1a)
$\displaystyle{\cal O}_{q\chi 1}^{P}$
$\displaystyle=(\overline{q}i\gamma_{5}q)(\overline{\chi}\chi),$
$\displaystyle{\cal O}_{q\chi 2}^{P}$
$\displaystyle=(\overline{q}\gamma_{5}q)(\overline{\chi}\gamma_{5}\chi),$ (1b)
$\displaystyle{\cal O}_{q\chi 1}^{V}$
$\displaystyle=(\overline{q}\gamma_{\mu}q)(\overline{\chi}\gamma_{\mu}\chi),\,(\times)$
$\displaystyle{\cal O}_{q\chi 2}^{V}$
$\displaystyle=(\overline{q}\gamma_{\mu}q)(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi),$
(1c) $\displaystyle{\cal O}_{q\chi 1}^{A}$
$\displaystyle=(\overline{q}\gamma_{\mu}\gamma_{5}q)(\overline{\chi}\gamma_{\mu}\chi),\,(\times)$
$\displaystyle{\cal O}_{q\chi 2}^{A}$
$\displaystyle=(\overline{q}\gamma_{\mu}\gamma_{5}q)(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi),$
(1d) $\displaystyle{\cal O}_{q\chi 1}^{T}$
$\displaystyle=(\overline{q}\sigma_{\mu\nu}q)(\overline{\chi}\sigma_{\mu\nu}\chi),\,(\times)$
$\displaystyle{\cal O}_{q\chi 2}^{T}$
$\displaystyle=(\overline{q}\sigma_{\mu\nu}q)(\overline{\chi}\sigma_{\mu\nu}\gamma_{5}\chi),\,(\times)$
(1e)

where the quark flavor indices have been omitted for notational simplicity,
but are understood to be those required for the specific FCNC process
throughout the paper. The “$(\times)$” indicates the accompanying operator
vanishes for the Majorana DM case due to the fermion bilinear identity,
$\overline{\chi}\Gamma\chi=-\overline{\chi^{\tt C}}\Gamma\chi^{\tt C}$ for
$\Gamma\in\\{\gamma_{\mu},\sigma_{\mu\nu},\sigma_{\mu\nu}\gamma_{5}\\}$, where
$\chi^{\tt C}\equiv-i\gamma^{2}\chi^{*}$ is the charge conjugation of $\chi$.

Scalar case: The leading order operators for the scalar DM appear at dimension
5 and 6,333They would all arise at dimension 6 in the SMEFT framework
Kamenik:2011vy . A specific renormalizable model realization for the operators
${\cal O}_{q\phi}^{S}$ and ${\cal O}_{q\phi}^{P}$ is given in appendix E. The
realization for the operators ${\cal O}_{qX}^{S}$ and ${\cal O}_{qX}^{P}$
involving vector DM is also given there. they can be parametrized in the
following manner,

$\displaystyle{\cal O}_{q\phi}^{S}$ $\displaystyle=$
$\displaystyle(\overline{q}q)(\phi^{\dagger}\phi),$ (2a) $\displaystyle{\cal
O}_{q\phi}^{P}$ $\displaystyle=$
$\displaystyle(\overline{q}i\gamma_{5}q)(\phi^{\dagger}\phi),$ (2b)
$\displaystyle{\cal O}_{q\phi}^{V}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma^{\mu}q)(\phi^{\dagger}i\overleftrightarrow{\partial_{\mu}}\phi),\,(\times)$
(2c) $\displaystyle{\cal O}_{q\phi}^{A}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma^{\mu}\gamma_{5}q)(\phi^{\dagger}i\overleftrightarrow{\partial_{\mu}}\phi),\,(\times).$
(2d)

Once again the implicit quark flavor indices should be understood. The symbol
“$(\times)$” indicates the related operator vanishes for real scalar DM, and
the double arrow derivative is defined as
$A\overleftrightarrow{\partial_{\mu}}B\equiv
A(\partial_{\mu}B)-(\partial_{\mu}A)B$.

Vector case A: For the vector DM, we consider separately two cases: when the
DM field is represented by the four-vector potential $X_{\mu}$ (scenario A) or
by the field strength tensor $X_{\mu\nu}$ (scenario B). For scenario A, we
find there are 4 independent dim-5 operators with the field content
$\bar{q}qX^{\dagger}X$ and 12 dim-6 operators with the content
$\bar{q}qX^{\dagger}XD_{\mu}$. By requiring the flavor diagonal operators to
be self-conjugate, we can parametrize those operators in the following way,

$\displaystyle{\cal O}_{qX}^{S}$ $\displaystyle=$
$\displaystyle(\overline{q}q)(X_{\mu}^{\dagger}X^{\mu}),$ (3a)
$\displaystyle{\cal O}_{qX}^{P}$ $\displaystyle=$
$\displaystyle(\overline{q}i\gamma_{5}q)(X_{\mu}^{\dagger}X^{\mu}),$ (3b)
$\displaystyle{\cal O}_{qX1}^{T}$ $\displaystyle=$ $\displaystyle{i\over
2}(\overline{q}\sigma^{\mu\nu}q)(X_{\mu}^{\dagger}X_{\nu}-X_{\nu}^{\dagger}X_{\mu}),\,(\times)$
(3c) $\displaystyle{\cal O}_{qX2}^{T}$ $\displaystyle=$ $\displaystyle{1\over
2}(\overline{q}\sigma^{\mu\nu}\gamma_{5}q)(X_{\mu}^{\dagger}X_{\nu}-X_{\nu}^{\dagger}X_{\mu}),\,(\times)$
(3d) $\displaystyle{\cal O}_{qX1}^{V}$ $\displaystyle=$ $\displaystyle{1\over
2}[\overline{q}\gamma_{(\mu}i\overleftrightarrow{D_{\nu)}}q](X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(3e) $\displaystyle{\cal O}_{qX2}^{V}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}q)\partial_{\nu}(X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(3f) $\displaystyle{\cal O}_{qX3}^{V}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}q)(X_{\rho}^{\dagger}\overleftrightarrow{\partial_{\nu}}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},$
(3g) $\displaystyle{\cal O}_{qX4}^{V}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma^{\mu}q)(X_{\nu}^{\dagger}i\overleftrightarrow{\partial_{\mu}}X^{\nu}),\,(\times)$
(3h) $\displaystyle{\cal O}_{qX5}^{V}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}q)i\partial_{\nu}(X^{\mu\dagger}X^{\nu}-X^{\nu\dagger}X^{\mu}),\,(\times)$
(3i) $\displaystyle{\cal O}_{qX6}^{V}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}q)i\partial_{\nu}(X^{\dagger}_{\rho}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},\,(\times)$
(3j) $\displaystyle{\cal O}_{qX1}^{A}$ $\displaystyle=$ $\displaystyle{1\over
2}[\overline{q}\gamma_{(\mu}\gamma_{5}i\overleftrightarrow{D_{\nu)}}q](X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(3k) $\displaystyle{\cal O}_{qX2}^{A}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}\gamma_{5}q)\partial_{\nu}(X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(3l) $\displaystyle{\cal O}_{qX3}^{A}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}\gamma_{5}q)(X_{\rho}^{\dagger}\overleftrightarrow{\partial_{\nu}}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},$
(3m) $\displaystyle{\cal O}_{qX4}^{A}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma^{\mu}\gamma_{5}q)(X_{\nu}^{\dagger}i\overleftrightarrow{\partial_{\mu}}X^{\nu}),\,(\times)$
(3n) $\displaystyle{\cal O}_{qX5}^{A}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}\gamma_{5}q)i\partial_{\nu}(X^{\mu\dagger}X^{\nu}-X^{\nu\dagger}X^{\mu}),\,(\times)$
(3o) $\displaystyle{\cal O}_{qX6}^{A}$ $\displaystyle=$
$\displaystyle(\overline{q}\gamma_{\mu}\gamma_{5}q)i\partial_{\nu}(X^{\dagger}_{\rho}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},\,(\times)$
(3p)

where the current
$\overline{q}\gamma_{(\mu}i\overleftrightarrow{D_{\nu)}}q\equiv\overline{q}\gamma_{\mu}i\overleftrightarrow{D_{\nu}}q+\mu\leftrightarrow\nu$,
and similarly for the current
$\overline{q}\gamma_{(\mu}\gamma_{5}i\overleftrightarrow{D_{\nu)}}q$. The
symbol “$(\times)$” indicates the corresponding operator vanishes for real
vector DM. Other operators with different Lorentz contractions or derivatives
are not independent and can always be reduced to those given above by using
the Dirac gamma identities (DI), integration by parts (IBP), and equation of
motions (EoM). The construction of the above dim-6 operators with a derivative
is expanded upon in appendix B. We have also checked our result by using the
Hilbert series method with a modified conformal representation for the vector
DM field Henning:2017fpj . 444 However, when using the Basisgen package
Criado:2019ugp , we only find 4 dim-6 operators containing a derivative and
the other 8 operators are missing.

Using the above operators to calculate the amplitudes for physical processes
results in rates that are divergent in the limit of massless DM particles.
This divergence originates from the longitudinal part in the polarization sum
and is a well known problem. One way to deal with this problem is to assume
that the vectors are gauge bosons of some dark symmetry and gauge invariance
under that symmetry forbids the direct appearance of the field $X_{\mu}$.
These would appear instead in covariant derivatives acting on other dark
matter fields which are not present in our effective Lagrangian. The net
effect of such a scenario is that $X_{\mu}$ adquires mass by some Higgs
mechanism Kamenik:2011vy ; Williams:2011qb  and the effective operators
inherit a coefficient that vanishes in the limit of massless $X_{\mu}$. An
example of how this could work in a specific model is given in appendix E.
Operationally, for our numerical analysis, we require the Wilson coefficients
for the above operators to contain an explicit factor of the DM mass to the
minimum power necessary to cancel potential divergences as $m\to 0$.

Vector case B: With the vector DM entering through field strength tensors the
minimal dimensionality of the operators is 7.555If working with SMEFT
assumption, the resulting operators would be at dimension 8. In this case two
explicit DM fields are required again for processes with two DM particles
because a single field strength tensor referring to a dark non-abelian
symmetry could not couple to a quark current that is not charged under the
dark group. The operators in this scenario produce amplitudes that are well
behaved as $m\to 0$ and we find there are 6 operators as follows,

$\displaystyle\tilde{\cal O}_{qX1}^{S}$ $\displaystyle=$
$\displaystyle(\overline{q}q)X_{\mu\nu}^{\dagger}X^{\mu\nu},$ (4a)
$\displaystyle\tilde{\cal O}_{qX2}^{S}$ $\displaystyle=$
$\displaystyle(\overline{q}q)X_{\mu\nu}^{\dagger}\tilde{X}^{\mu\nu},$ (4b)
$\displaystyle\tilde{\cal O}_{qX1}^{P}$ $\displaystyle=$
$\displaystyle(\overline{q}i\gamma_{5}q)X_{\mu\nu}^{\dagger}X^{\mu\nu},$ (4c)
$\displaystyle\tilde{\cal O}_{qX2}^{P}$ $\displaystyle=$
$\displaystyle(\overline{q}i\gamma_{5}q)X_{\mu\nu}^{\dagger}\tilde{X}^{\mu\nu},$
(4d) $\displaystyle\tilde{\cal O}_{qX1}^{T}$ $\displaystyle=$
$\displaystyle{i\over
2}(\overline{q}\sigma^{\mu\nu}q)(X^{\dagger}_{\mu\rho}X^{\rho}_{\,\nu}-X^{\dagger}_{\nu\rho}X^{\rho}_{\,\mu}),\,(\times)$
(4e) $\displaystyle\tilde{\cal O}_{qX2}^{T}$ $\displaystyle=$
$\displaystyle{1\over
2}(\overline{q}\sigma^{\mu\nu}\gamma_{5}q)(X^{\dagger}_{\mu\rho}X^{\rho}_{\,\nu}-X^{\dagger}_{\nu\rho}X^{\rho}_{\,\mu}).\,(\times)$
(4f)

The dual field strength is defined as $\tilde{X}^{\mu\nu}={1\over
2}\epsilon^{\mu\nu\rho\sigma}X_{\alpha\beta}$, and the symbol “$(\times)$”
denotes an operator that vanishes for real vector fields.

The method used to construct the LEFT quark-DM interactions given above can
also be used to obtain the corresponding lepton-DM interactions. We list those
in appendix C for completeness.

In the literature, Ref. Kumar:2013iva  provides a list of LEFT operators for
the fermion, scalar, and vector DM (scenario A) cases coupled to flavor
diagonal currents. The fermion and scalar cases agree with our list, but the
vector case does not.666The list of dim-6 operators with a derivative in
Kumar:2013iva  is not complete. For example, it does not contain operators
corresponding to ${\cal O}_{qX3}^{V(A)}$. Here we discuss the more general
case with flavor non-diagonal operators relevant to the FCNC processes we
study.

As mentioned above, we only impose the unbroken $SU(3)_{\rm c}\times U(1)_{\rm
em}$ symmetry to obtain the LEFT operators. If we instead assume that all
mediators are far beyond the weak scale, we can start from a SMEFT with the
full SM gauge symmetry $SU(3)_{\rm c}\times SU(2)_{L}\times U(1)_{Y}$ and
obtain the corresponding LEFT in the Higgs phase. The relevant DM EFT
operators in this picture can be found in Kamenik:2011vy ; Brod:2017bsw ;
Criado:2021trs ; Aebischer:2022wnl . 777 For the vector DM case, we find that
the DM-quark (lepton) operators given in Tab. 8 in Criado:2021trs  are neither
independent nor complete. For example, for vector DM coupled to the down-type
right-handed quark current, that paper provides 2 operators:
$(\rho^{\dagger}_{\mu}D_{\nu}\rho^{\mu})(\bar{d}\gamma^{\nu}d)$ and
$(\rho^{\dagger}_{\mu}D_{\nu}\rho^{\nu})(\bar{d}\gamma^{\mu}d)$. The first one
is equivalent to ${\cal O}_{qX1}^{V+A}$ in our list, while the second one
vanishes from the on-shell condition $\partial_{\mu}\rho^{\mu}=0$.  The FCNC
$B$ and $K$ meson decays into fermion DM (or similar invisible particles like
sterile neutrinos) have been studied in Kamenik:2011vy ; Li:2020dpc ;
Felkl:2021uxi ; Li:2019fhz ; Deppisch:2020oyx ; He:2021yoz . For this reason,
in the following we restrict ourselves to the scalar and vector DM cases and
investigate the experimental sensitivity to the interactions in Eqs. (2-4) in
detail. We will first consider $B$ meson decay in the next section, followed
by the $K$ meson decay after that.

##  III $B\to(K,K^{*},\pi,\rho)$+DM+DM

To calculate the decay rate for the $B\to M$ transition (where $M$ denotes
either a pseudo-scalar meson $P=K,\pi$ or a vector meson $V=K^{*},\rho$) from
the effective interactions in Eqs. (2-4), we first need to know the hadronic
transition matrix elements $\langle M|\bar{q}\Gamma b|B\rangle$. These are
usually parametrized in terms of scalar form factors associated with each
possible allowed Lorentz structure. The Lorentz structures can be organized
according to parity and charge conjugation. While some of the form factors can
be determined from experimental data, others require theoretical models for
the non-perturbative aspects of QCD. In the following subsections, we first
collect the relevant form factors and their determination using light-cone sum
rules. We then consider the decay rates for both scalar and vector DM
scenarios and the implications for the parameter space.

###  III.1 Form factors

We follow the parametrization of $B$ meson form factors in Gubernari:2018wyi .
For the $B\to P(J^{P}=0^{-})$ transition with a final state pseudo-scalar
$P=\pi,\,K$, the non-vanishing hadronic matrix elements from the scalar,
vector, and tensor quark currents are parametrized by the form factors
$f_{0},\ f_{+}$, and $f_{T}$,

$\displaystyle\langle P(k)|{\bar{q}}b|B(p)\rangle$ $\displaystyle=$
$\displaystyle{m_{B}^{2}-m_{P}^{2}\over m_{b}-m_{q}}f_{0}(q^{2}),$ (5a)
$\displaystyle\langle P(k)|{\bar{q}}\gamma^{\mu}b|B(p)\rangle$
$\displaystyle=$
$\displaystyle\left[\left(p+k\right)^{\mu}-\frac{m_{B}^{2}-m_{P}^{2}}{q^{2}}q^{\mu}\right]f_{+}(q^{2})+\frac{m_{B}^{2}-m_{P}^{2}}{q^{2}}q^{\mu}f_{0}(q^{2}),$
(5b) $\displaystyle\langle P(k)|{\bar{q}}\sigma^{\mu\nu}b|B(p)\rangle$
$\displaystyle=$
$\displaystyle\frac{2i}{m_{B}+m_{P}}(p^{\mu}q^{\nu}-p^{\nu}q^{\mu})f_{T}(q^{2}),$
(5c)

where $k$ and $p$ are the 4-momenta of $P$ and $B$ respectively,
$q^{\mu}=p^{\mu}-k^{\mu}$, $m_{B}$ and $m_{P}$ are the masses of the initial
$B$ meson and final state $P$ meson, and $m_{b}$ and $m_{q}$ are the masses of
the quarks appearing in the currents. In the $q^{2}\to 0$ limit,
$f_{+}(0)=f_{0}(0)$, and we follow the light-cone sum rule (LCSR) methods
Ball:2004ye  to parametrize the dependence on the momentum transfer $s\equiv
q^{2}$ as,

$\displaystyle f_{0}(s)={r_{2}\over 1-s/m_{\rm
fit}^{2}},\,f_{+(T)}^{\pi}(s)={r_{1}\over 1-s/m_{R}^{2}}+{r_{2}\over
1-s/m_{\rm fit}^{2}},\,f_{+(T)}^{K}(s)={r_{1}\over
1-s/m_{R}^{2}}+{r_{2}\over\left(1-s/m_{R}^{2}\right)^{2}}.$ (6)

Above, $r_{1,2}$, $m_{R}^{2}$, and $m_{\rm fit}^{2}$, are parameters with the
preferred values given in Ball:2004ye  and collected in appendix D for
reference.

For the transition into a vector meson $V$, $B\to V(J^{P}=1^{-})$ with
$V=\rho,\,K^{*}$, the non-vanishing form factors $V_{0},\ A_{0,1,2,3},\
T_{1,2,3}$ are defined as

$\displaystyle\langle V(k)|{\bar{q}}\gamma_{5}b|B(p)\rangle$ $\displaystyle=$
$\displaystyle-i\epsilon_{V,\nu}^{*}q^{\nu}{2m_{V}\over m_{b}+m_{q}}A_{0},$
(7a) $\displaystyle\langle V(k)|{\bar{q}}\gamma^{\mu}b|B(p)\rangle$
$\displaystyle=$
$\displaystyle\epsilon^{\mu\nu\rho\sigma}\epsilon_{V,\nu}^{*}p_{\rho}k_{\sigma}\frac{2}{m_{B}+m_{V}}V_{0},$
(7b) $\displaystyle\langle V(k)|{\bar{q}}\gamma^{\mu}\gamma_{5}b|B(p)\rangle$
$\displaystyle=$ $\displaystyle
i\epsilon_{V,\nu}^{*}\left[g^{\mu\nu}(m_{B}+m_{V})A_{1}-\frac{(p+k)^{\mu}q^{\nu}}{m_{B}+m_{V}}A_{2}-q^{\mu}q^{\nu}\frac{2m_{V}}{q^{2}}(A_{3}-A_{0})\right],$
(7c) $\displaystyle\langle V(k)|{\bar{q}}\sigma_{\mu\nu}b|B(p)\rangle$
$\displaystyle=$ $\displaystyle
i\epsilon_{\mu\nu\rho\sigma}\epsilon_{V,\alpha}^{*}\left\\{g^{\alpha\rho}(p+k)^{\sigma}T_{1}-g^{\alpha\rho}q^{\sigma}{m_{B}^{2}-m_{V}^{2}\over
q^{2}}(T_{1}-T_{2})\right.$ (7d) $\displaystyle+$
$\displaystyle\left.2q^{\alpha}p^{\rho}k^{\sigma}\left[{1\over
m_{B}^{2}-m_{V}^{2}}T_{3}-{1\over q^{2}}(T_{1}-T_{2})\right]\right\\},$

where $\epsilon_{V}$ is the polarization vector of the spin-one meson and
$m_{V}$ its mass.888In literature, the tensor current is usually parametrized
by multiplying by the four-momentum $q_{\nu}$. From the tensor current in Eq.
(7), they can be directly calculated to take the form,  $\displaystyle\langle
V(k)|{\bar{q}}i\sigma^{\mu\nu}q_{\nu}b|B(p)\rangle$ $\displaystyle=$
$\displaystyle
2\epsilon^{\mu\nu\rho\sigma}\epsilon_{V,\nu}^{*}p_{\rho}k_{\sigma}T_{1},$ (8a)
$\displaystyle\langle
V(k)|{\bar{q}}i\sigma^{\mu\nu}\gamma_{5}q_{\nu}b|B(p)\rangle$ $\displaystyle=$
$\displaystyle
i\epsilon_{V,\nu}^{*}\left\\{\left[g^{\mu\nu}(m_{B}^{2}-m_{V}^{2})-(p+k)^{\mu}q^{\nu}\right]T_{2}+q^{\nu}\left[q^{\mu}-\frac{q^{2}(p+k)^{\mu}}{m_{B}^{2}-m_{V}^{2}}\right]T_{3}\right\\}.$
(8b) The axial-vector and pseudo-scalar current matrix elements can be related
by EoM via the relation $i\partial_{\mu}\langle
P|\bar{q}\gamma^{\mu}\gamma_{5}b|B\rangle=\langle
P|\bar{q}i(\not{D}+\overleftarrow{\not{D}})\gamma_{5}b|B\rangle$. Equivalently
in momentum space, $q_{\mu}\langle
P(k)|\bar{q}\gamma^{\mu}\gamma_{5}b|B(p)\rangle=-(m_{b}+m_{q})\langle
P(k)|\bar{q}\gamma_{5}b|B(p)\rangle$, which implies that the form factor
$A_{3}$ is a redundant and can be expressed in terms of $A_{1}$ and $A_{2}$
as,

$\displaystyle
A_{3}\equiv\frac{m_{B}+m_{V}}{2m_{V}}A_{1}-\frac{m_{B}-m_{V}}{2m_{V}}A_{2}.$
(9)

It is common practice to replace $A_{2}$ and $T_{3}$ by

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

where the K$\ddot{a}$llen function $\lambda(x,y,x)$ is the usual,

$\displaystyle\lambda(x,y,z)\equiv x^{2}+y^{2}+z^{2}-2(xy+yz+zx).$ (11)

In the $q^{2}\to 0$ limit, some of form factors are related with as follows,

$\displaystyle A_{0}(0)=A_{3}(0),\quad T_{1}(0)=T_{2}(0),\quad
A_{12}(0)=\frac{m_{B}^{2}-m_{V}^{2}}{8m_{B}m_{V}}A_{0}(0).$ (12)

Relabeling the form factors
$\\{A_{0},A_{1},A_{12},V_{0},T_{1},T_{2},T_{23}\\}$ as $F_{1,2,3,4,5,6,7}$,
their momentum transfer dependence can be parametrized as Bharucha:2015bzk ,

$\displaystyle F_{i}(s)={1\over
1-s/m_{R,i}^{2}}\sum_{k}\alpha_{k}^{i}[z(s)-z(0)]^{k},\quad
z(s)\equiv{\sqrt{s_{+}-s}-\sqrt{s_{+}-s_{0}}\over\sqrt{s_{+}-s}+\sqrt{s_{+}-s_{0}}},$
(13)

where $s_{\pm}\equiv(m_{B}\pm m_{V})^{2}$ and $s_{0}\equiv
s_{+}(1-\sqrt{1-s_{-}/s_{+}})$. $M_{R,i}$ are the resonance masses associated
with the transition modes and are taken from Tab. 3 in Bharucha:2015bzk . The
parameters $\alpha_{k}^{i}$ are truncated at quadratic order in $z$, $k_{\rm
max}=2$, so that three fit parameters
$\alpha_{0}^{i},\alpha_{1}^{i},\alpha_{3}^{i}$ are needed for each form factor
$i$. They are given in Tab. 14 of Bharucha:2015bzk  and we collect them in
appendix D for convenience.

###  III.2  $B\to M+\phi\phi$ with scalar DM $\phi$

For the quark-scalar DM interactions in Eq. (2), with the hadronic matrix
elements given in Eq. (5) and Eq. (7), the non-vanishing amplitudes for the
processes $B(p)\to P(k)\phi(k_{1})\phi^{*}(k_{2})$ and $B(p)\to
V(k)\phi(k_{1})\phi^{*}(k_{2})$ take the following general form,

$\displaystyle i{\cal M}_{B\to P\phi\phi}$ $\displaystyle=$ $\displaystyle
C_{q\phi}^{S,xb}\langle
P(k)|\overline{q}_{x}b|B(p)\rangle+C_{q\phi}^{V,xb}(k_{1}^{\mu}-k_{2}^{\mu})\langle
P(k)|\overline{q}_{x}\gamma_{\mu}b|B(p)\rangle,$ (14a) $\displaystyle i{\cal
M}_{B\to V\phi\phi}$ $\displaystyle=$ $\displaystyle C_{q\phi}^{P,xb}\langle
V(k)|\overline{q}_{x}i\gamma_{5}b|B(p)\rangle+C_{q\phi}^{V,xb}(k_{1}^{\mu}-k_{2}^{\mu})\langle
V(k)|\overline{q}_{x}\gamma_{\mu}b|B(p)\rangle$ (14b) $\displaystyle+$
$\displaystyle C_{q\phi}^{A,xb}(k_{1}^{\mu}-k_{2}^{\mu})\langle
V(k)|\overline{q}_{x}\gamma_{\mu}\gamma_{5}b|B(p)\rangle,$

where $x=d,s$ is a quark flavor label characterizing the final state meson
$P(V)=\pi,K(\rho,K^{*})$. Using the hadronic matrix elements in Eq. (5) and
Eq. (7), with the help of Feyncalc Shtabovenko:2016sxi , the differential
decay widths take the following compact form

$\displaystyle{d\Gamma_{B\to P\phi\phi}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{(m_{B}^{2}-m_{P}^{2})^{2}\over
256\pi^{3}m_{B}^{3}(m_{b}-m_{q_{x}})^{2}}\lambda^{1\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{1\over
2}(m^{2},s)f_{0}^{2}\left|C_{q\phi}^{S,xb}\right|^{2}$ (15) $\displaystyle+$
$\displaystyle{1\over 768\pi^{3}m_{B}^{3}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over
2}(m^{2},s)f_{+}^{2}\left|C_{q\phi}^{V,xb}\right|^{2},$
$\displaystyle{d\Gamma_{B\to V\phi\phi}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{1\over 256\pi^{3}m_{B}^{3}(m_{b}+m_{q_{x}})^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)A_{0}^{2}\left|C_{q\phi}^{P,xb}\right|^{2}$ (16) $\displaystyle+$
$\displaystyle{s\over 384\pi^{3}m_{B}^{3}(m_{B}+m_{V})^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)V_{0}^{2}\left|C_{q\phi}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{1\over 384\pi^{3}m_{B}^{3}}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)\left[(m_{B}+m_{V})^{2}sA_{1}^{2}+32m_{B}^{2}m_{V}^{2}A_{12}^{2}\right]\left|C_{q\phi}^{A,xb}\right|^{2},$

where $\kappa(m^{2},s)$ is a kinetic “endpoint” function defined as

$\displaystyle\kappa(m^{2},s)\equiv 1-{4m^{2}\over s},$ (17)

and the $q^{2}$ dependence of the form factors is left implicit for notational
simplicity. In the above results, there are no interference effects between
any pair of operators because the relevant hadronic matrix elements have
different parity and/or charge conjugation properties and thus cannot mix. The
dependence on the DM mass enters only through kinematics, unlike the vector DM
case discussed below. For the case of real scalar DM, only the scalar and
pseudo-scalar quark currents appear (due to some operators vanishing as noted
in Eq. (2)); there is also an additional factor of two in the decay width.

Figure 1: Normalized differential decay width for $B\to K^{(*)+}\phi\phi$
decay from different types of operators. Left panel: $m=100\,\rm MeV$; Right
panel: $m=1\,\rm GeV$.

Figure 2: Constraints on the effective new physics scale for each operator
involving a $b$ quark as a function of the DM mass $m$ from all possible
relevant $B$ meson decay channels.

The different operators result in different $q^{2}$ distributions in the $B\to
M\not{E}$ decay. For example, we illustrate the normalized distributions in
the $b\to s$ transition for different operators in Fig. 1. 999In the figures
we add a subscript to differentiate the type of DM in question, but in the
text we refer generically to the mass of any DM particle as $m$. In the left
(right) panel, we consider scalar DM mass $m=100\,\rm MeV(1\,GeV)$ and in both
cases a solid (dashed) line is used for the $B\to K^{+}\phi\phi$ ($K\to
K^{*+}\phi\phi$) channels. One can clearly see that the distribution varies
significantly between operators and as a function of the DM mass. This feature
could be exploited to differentiate the various cases in the upcoming
experimental search from Belle II.

To quantify the constraints set on the parameter by the current experimental
bounds listed in Tab. 2, we use an effective scale $\Lambda_{\rm eff}$
associated with each operator from dimensional analysis:
$C_{q\phi}^{S(P)}\equiv\Lambda_{\rm eff}^{-1}$ and
$C_{q\phi}^{V(A)}\equiv\Lambda_{\rm eff}^{-2}$. Fig. 2 shows the current
experimental sensitivity in the $m$-$\Lambda$ plane for each operator with
$(bs)$ and $(bd)$ flavor changing quark combinations. The solid (dashed) lines
correspond to the constraints from charged (neutral) decay modes for both
pseudo-scalar and vector final state mesons. The largest possible DM mass is
restricted by the kinematic relation $m\leq(m_{B}-m_{M})/2$ as reflected in
each panel. Generally, for the scalar and vector current operators on the left
two panels, the charged mode $B^{+}\to K^{+}$ sets the stronger constraints
for $(bs)$ transitions and the neutral mode $B^{0}\to\pi^{0}$ for the $(bd)$
transitions. For the pseudo-scalar and axial-vector quark current operators
shown on the right two panels, the stronger bounds for the $(bs)$ and $(bd)$
transitions are set respectively by the neutral mode $B^{0}\to K^{*0}$ and the
charged mode $B^{+}\to\rho^{+}$. This feature is just a reflection of the
current experimental bounds as can be seen in Tab. 2.

The situation depicted in Fig. 2 will, of course, be modified by experimental
considerations. For example, Belle II has reported with its current
measurement of $B^{+}\to K^{+}\not{E}$ a signal efficiency that varies with
$q^{2}$, peaking at low values and becoming very small for $q^{2}\gtrsim
12~{}{\rm GeV}^{2}$ Belle-II:2021rof . Fig. 1 then suggests that searches for
very light DM will be more sensitive than searches for heavier DM. Inclusion
of this experimental sensitivity changes the corresponding constraint, the
solid red line in Fig. 2, into the solid gray line (where E.E. stands for
“experimental efficiency”). As can be seen the sensitivity in this case is
weaker by a factor of 2-5 for $m=0$, and drops rapidly as $m$ increases.

###  III.3 $B\to M+XX$ with vector DM $X$: scenario A

For spin one DM, we first consider scenario A in which the operators were
constructed using the vector field as given in Eq. (3). Two of these operators
involve a covariant derivative acting on the quark current, ${\cal
O}_{qX1}^{V}$ and ${\cal O}_{qX1}^{A}$. To evaluate their contribution one
needs form factors that have not been studied before. Within specific quark
models for mesons one could estimate these form factors by replacing the
derivatives with the corresponding quark momentum. As this is beyond the scope
of this paper, we simply ignore these two operators in the following numerical
study. The non-vanishing amplitudes for the processes $B(p)\to
P(k)X(k_{1})X^{*}(k_{2})$ and $B(p)\to V(k)X(k_{1})X^{*}(k_{2})$ from the
remaining operators in Eq. (3) take the following form,

$\displaystyle i{\cal M}_{B\to PXX}^{A}$ $\displaystyle=$
$\displaystyle\epsilon^{*}_{\rho}(k_{1})\epsilon^{*}_{\sigma}(k_{2})\left\\{g^{\rho\sigma}C_{qX}^{S,xb}\langle
P(k)|\bar{q}_{x}b|B(p)\rangle\right.$ (18a) $\displaystyle+$
$\displaystyle{i\over
2}\left(2g^{\mu\rho}g^{\nu\sigma}C_{qX1}^{T,xb}+\epsilon^{\mu\nu\rho\sigma}C_{qX2}^{T,xb}\right)\langle
P(k)|\bar{q}_{x}\sigma_{\mu\nu}b|B(p)\rangle$ $\displaystyle+$
$\displaystyle\left.\left[i(g^{\mu\rho}k_{1\sigma}+g^{\mu\sigma}k_{2\rho})C_{qX2}^{V,xb}-i\epsilon^{\mu\nu\rho\sigma}(k_{1}-k_{2})_{\nu}C_{qX3}^{V,xb}+g^{\rho\sigma}(k_{1}-k_{2})^{\mu}C_{qX4}^{V,xb}\right.\right.$
$\displaystyle-$
$\displaystyle\left.\left.(g^{\mu\rho}k_{1\sigma}-g^{\mu\sigma}k_{2\rho})C_{qX5}^{V,xb}-\epsilon^{\mu\nu\rho\sigma}(k_{1}+k_{2})_{\nu}C_{qX6}^{V,xb}\right]\langle
P(k)|\bar{q}_{x}\gamma_{\mu}b|B(p)\rangle\right\\},$ $\displaystyle i{\cal
M}_{B\to VXX}^{A}$ $\displaystyle=$
$\displaystyle\epsilon^{*}_{\rho}(k_{1})\epsilon^{*}_{\sigma}(k_{2})\left\\{g^{\rho\sigma}C_{qX}^{P,xb}\langle
V(k)|\bar{q}_{x}i\gamma_{5}b|B(p)\rangle\right.$ (18b) $\displaystyle+$
$\displaystyle{i\over
2}\left(2g^{\mu\rho}g^{\nu\sigma}C_{qX1}^{T,xb}+\epsilon^{\mu\nu\rho\sigma}C_{qX2}^{T,xb}\right)\langle
V(k)|\bar{q}_{x}\sigma_{\mu\nu}b|B(p)\rangle$ $\displaystyle+$
$\displaystyle\left.\left[i(g^{\mu\rho}k_{1\sigma}+g^{\mu\sigma}k_{2\rho})C_{qX2}^{V,xb}-i\epsilon^{\mu\nu\rho\sigma}(k_{1}-k_{2})_{\nu}C_{qX3}^{V,xb}+g^{\rho\sigma}(k_{1}-k_{2})^{\mu}C_{qX4}^{V,xb}\right.\right.$
$\displaystyle-$
$\displaystyle\left.(g^{\mu\rho}k_{1\sigma}-g^{\mu\sigma}k_{2\rho})C_{qX5}^{V,xb}-\epsilon^{\mu\nu\rho\sigma}(k_{1}+k_{2})_{\nu}C_{qX6}^{V,xb}\right]\langle
V(k)|\bar{q}_{x}\gamma_{\mu}b|B(p)\rangle$ $\displaystyle+$
$\displaystyle\left.\left[i(g^{\mu\rho}k_{1\sigma}+g^{\mu\sigma}k_{2\rho})C_{qX2}^{A,xb}-i\epsilon^{\mu\nu\rho\sigma}(k_{1}-k_{2})_{\nu}C_{qX3}^{A,xb}+g^{\rho\sigma}(k_{1}-k_{2})^{\mu}C_{qX4}^{A,xb}\right.\right.$
$\displaystyle-$
$\displaystyle\left.\left.(g^{\mu\rho}k_{1\sigma}-g^{\mu\sigma}k_{2\rho})C_{qX5}^{A,xb}-\epsilon^{\mu\nu\rho\sigma}(k_{1}+k_{2})_{\nu}C_{qX6}^{A,xb}\right]\langle
V(k)|\bar{q}_{x}\gamma_{\mu}\gamma_{5}b|B(p)\rangle\right\\},$

For the final state with a pseudo-scalar meson $P$, using the hadronic matrix
elements in Eq. (5), leads to the differential decay width

$\displaystyle{d\Gamma_{B\to PXX}^{A}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{(m_{B}^{2}-m_{P}^{2})^{2}(s^{2}-4m^{2}s+12m^{4})\over
1024\pi^{3}m_{B}^{3}(m_{b}-m_{q_{x}})^{2}m^{4}}\lambda^{1\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{1\over
2}(m^{2},s)f_{0}^{2}\left|C_{qX}^{S,xb}\right|^{2}$ (19) $\displaystyle+$
$\displaystyle{s(s+4m^{2})\over
3072\pi^{3}m_{B}^{3}(m_{B}+m_{P})^{2}m^{4}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over
2}(m^{2},s)f_{T}^{2}\left|C_{qX1}^{T,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s+2m^{2}\over
768\pi^{3}m_{B}^{3}(m_{B}+m_{P})^{2}m^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{1\over
2}(m^{2},s)f_{T}^{2}\left|C_{qX2}^{T,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s\over 3072\pi^{3}m_{B}^{3}m^{4}}\lambda^{1\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left[3(s-4m^{2})(m_{B}^{2}-m_{P}^{2})^{2}f_{0}^{2}+4m^{2}\lambda(m_{B}^{2},m_{P}^{2},s)f_{+}^{2}\right]\left|C_{qX2}^{V,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{1\over
768\pi^{3}m_{B}^{3}m^{2}}\lambda^{1\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left[6m^{2}(m_{B}^{2}-m_{P}^{2})^{2}f_{0}^{2}+(s-4m^{2})\lambda(m_{B}^{2},m_{P}^{2},s)f_{+}^{2}\right]\left|C_{qX3}^{V,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{s^{2}-4m^{2}s+12m^{4}\over
3072\pi^{3}m_{B}^{3}m^{4}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over
2}(m^{2},s)f_{+}^{2}\left|C_{qX4}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s(s+4m^{2})\over 3072\pi^{3}m_{B}^{3}m^{4}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over
2}(m^{2},s)f_{+}^{2}\left|C_{qX5}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s+2m^{2}\over 768\pi^{3}m_{B}^{3}m^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{1\over
2}(m^{2},s)f_{+}^{2}\left|C_{qX6}^{V,xb}\right|^{2}+\cdots$

We have dropped interference between different operators (represented by
“$\cdots$” above) as our numerical study will deal only with one operator at a
time.

For the final state with a vector meson $V$, using the form factors in Eq.
(7), we obtain

$\displaystyle{d\Gamma_{B\to VXX}^{A}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{s^{2}-4m^{2}s+12m^{4}\over
1024\pi^{3}m_{B}^{3}(m_{b}+m_{q_{x}})^{2}m^{4}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)A_{0}^{2}\left|C_{qX}^{P,xb}\right|^{2}$ (20) $\displaystyle+$
$\displaystyle{1\over 1536\pi^{3}m_{B}^{3}m^{4}s}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left\\{(s^{2}-16m^{4})\lambda(m_{B}^{2},m_{V}^{2},s)T_{1}^{2}+4m^{2}(s+2m^{2})\left[(m_{B}^{2}-m_{V}^{2})^{2}T_{2}^{2}+{8m_{B}^{2}m_{V}^{2}s\over(m_{B}+m_{V})^{2}}T_{23}^{2}\right]\right\\}\left|C_{qX1}^{T,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{1\over
1536\pi^{3}m_{B}^{3}m^{4}s}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left\\{4m^{2}(s+2m^{2})\lambda(m_{B}^{2},m_{V}^{2},s)T_{1}^{2}+(s^{2}-16m^{4})\left[(m_{B}^{2}-m_{V}^{2})^{2}T_{2}^{2}+{8m_{B}^{2}m_{V}^{2}s\over(m_{B}+m_{V})^{2}}T_{23}^{2}\right]\right\\}\left|C_{qX2}^{T,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{s^{2}\over
384\pi^{3}m_{B}^{3}(m_{B}+m_{V})^{2}m^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)V_{0}^{2}\left|C_{qX2}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s^{2}\over
384\pi^{3}m_{B}^{3}(m_{B}+m_{V})^{2}m^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{5\over
2}(m^{2},s)V_{0}^{2}\left|C_{qX3}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s(s^{2}-4m^{2}s+12m^{4})\over
1536\pi^{3}m_{B}^{3}(m_{B}+m_{V})^{2}m^{4}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)V_{0}^{2}\left|C_{qX4}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s^{2}(s+4m^{2})\over
1536\pi^{3}m_{B}^{3}(m_{B}+m_{V})^{2}m^{4}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)V_{0}^{2}\left|C_{qX5}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s(s+2m^{2})\over
384\pi^{3}m_{B}^{3}(m_{B}+m_{V})^{2}m^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)V_{0}^{2}\left|C_{qX6}^{V,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s\over 3072\pi^{3}m_{B}^{3}m^{4}}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left\\{3(s-4m^{2})\lambda(m_{B}^{2},m_{V}^{2},s)A_{0}^{2}+8m^{2}\left[(m_{B}+m_{V})^{2}sA_{1}^{2}+32m_{B}^{2}m_{V}^{2}A_{12}^{2}\right]\right\\}\left|C_{qX2}^{A,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{1\over
384\pi^{3}m_{B}^{3}m^{2}}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left\\{3m^{2}\lambda(m_{B}^{2},m_{V}^{2},s)A_{0}^{2}+(s-4m^{2})\left[(m_{B}+m_{V})^{2}sA_{1}^{2}+32m_{B}^{2}m_{V}^{2}A_{12}^{2}\right]\right\\}\left|C_{qX3}^{A,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{s^{2}-4m^{2}s+12m^{4}\over
1536\pi^{3}m_{B}^{3}m^{4}}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)\left[(m_{B}+m_{V})^{2}sA_{1}^{2}+32m_{B}^{2}m_{V}^{2}A_{12}^{2}\right]\left|C_{qX4}^{A,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{s(s+4m^{2})\over
1536\pi^{3}m_{B}^{3}m^{4}}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)\left[(m_{B}+m_{V})^{2}sA_{1}^{2}+32m_{B}^{2}m_{V}^{2}A_{12}^{2}\right]\left|C_{qX5}^{A,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{s+2m^{2}\over
384\pi^{3}m_{B}^{3}m^{2}}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)\left[(m_{B}+m_{V})^{2}sA_{1}^{2}+32m_{B}^{2}m_{V}^{2}A_{12}^{2}\right]\left|C_{qX6}^{A,xb}\right|^{2}+\cdots.$

where again we have dropped interference between different operators.

Figure 3: Constraints on the effective new physics scale for the 4 dim-5
operators ${\cal O}_{qX}^{S,P}$ and ${\cal O}_{qX1,2}^{T}$ as a function of
the DM mass $m$ from $B\to K(\pi)\not{E}$ channels.

It is evident that the above differential decay widths (and also decay widths)
diverge in the limit of vanishing DM mass, $m\to 0$. As discussed above, we
assume that each relevant Wilson coefficient also depends on the DM mass to
some power determined by the number of independent vector four-potentials that
cannot be reduced to field strength tensors. Operationally we use effective
scales $\Lambda_{\rm eff}$ defined as follows,

$\displaystyle C_{qX}^{S,P}\equiv{m^{2}\over\Lambda_{\rm eff}^{3}},\quad
C_{qX1,2}^{T}\equiv{m^{2}\over\Lambda_{\rm eff}^{3}},\quad
C_{qX2,4,5}^{V,A}\equiv{m^{2}\over\Lambda_{\rm eff}^{4}},\quad
C_{qX3,6}^{V,A}\equiv{m\over\Lambda_{\rm eff}^{3}}.$ (21)

In Eq. (21), $C_{qX3,6}^{V,A}$ depends linearly on $m$ because one of two
vector fields in the operators ${\cal O}_{qX3,6}^{V,A}$ can be rewritten as a
field strength tensor.

The possible origin of these mass factors is illustrated with an example in
appendix E. The divergence as $m\to 0$ also affects the kaon decay mode,
$K\to\pi XX$, where we use a parametrization similar to Eq. (21). As already
mentioned, the problem with the $m\to 0$ limit can be avoided by assuming
these vector particles are gauge bosons of a dark gauge symmetry and requiring
them to enter the LEFT as field strength tensors as in Eq. (4). We elaborate
on this second scenario for vector DM in the next subsection.

Figure 4: Constraints on the effective new physics scale for the 4 dim-6
operators ${\cal O}_{qX3,6}^{V,A}$ as a function of the DM mass $m$ from $B\to
K(\pi)\not{E}$ channels.

Figure 5: Constraints on the effective new physics scale for the 6 dim-6
operators ${\cal O}_{qX2,4,5}^{V,A}$ as a function of the DM mass $m$ from
$B\to K(\pi)\not{E}$ channels.

Fig. 3 shows the current experimental sensitivity in the $m$-$\Lambda_{\rm
eff}$ plane for the 4 dim-5 operators ${\cal O}_{qX}^{S,P}$ and ${\cal
O}_{qX1,2}^{T}$ following Eq. (21). For the quark scalar current operator
${\cal O}_{qX}^{S}$ (left upper panel), the charged mode ($B^{+}\to K^{+}XX$)
gives the strongest constraint for $(sb)$ flavor indices in the whole DM mass
range. Similarly, for $(db)$ flavor indices, the strongest constraint arises
from the neutral mode $B^{0}\to\pi^{0}XX$. For the quark pseudo-scalar current
operator ${\cal O}_{qX}^{P}$ (right upper panel), the constraints also exhibit
a similar behavior to those for the pseudo-scalar DM operator ${\cal
O}_{q\phi}^{P}$. However, due to the different dimensionality of the quark
(pseudo-)scalar current operators in the two cases, the numerical results are
very different for the two cases, as clearly seen in Fig. 2 and Fig. 3.

For the tensor operator ${\cal O}_{qX1}^{T}({\cal O}_{qX2}^{T})$ (lower two
panels), the neutral mode $B^{0}\to K^{*0}XX$ gives a stronger bound for DM
mass $m\lesssim 2.2(1.9)\,\rm GeV$ and the charged mode $B^{+}\to K^{+}XX$ for
$m\gtrsim 2.2(1.9)\,\rm GeV$ for $(sb)$ flavor indices. For $(db)$ flavor
indices instead, the neutral mode $B^{0}\to\pi^{0}XX$ gives a stronger bound
in the full DM mass range for ${\cal O}_{qX1}^{T}$. For ${\cal O}_{qX2}^{T}$
with $(db)$ flavor indices, it is the charged mode $B^{+}\to\rho^{+}XX$ that
gives a stronger bound for $m\lesssim 1.5\,\rm GeV$ and the neutral mode
$B^{0}\to\pi^{0}XX$ for $m\gtrsim 1.5\,\rm GeV$. The different behavior of the
operator ${\cal O}_{qX2}^{T}$ from the $B\to K(\pi)XX$ modes is due to a
quadratic (rather than quartic) inverse dependence on $m$ in the decay widths.
In all the four cases, the effective scale is constrained to be above a few
hundreds of GeV, validating our use of an EFT framework for this discussion.

Fig. 4 shows the constraints for the 4 dim-6 operators with quark
(axial-)vector currents, ${\cal O}_{qX3,6}^{V,A}$, and following Eq. (21). It
can be seen that these constraints observe a similar behavior to those for the
quark (axial-)vector current operators ${\cal O}_{q\phi}^{V,A}$ in the scalar
DM case (Fig. 2). The limit on $\Lambda_{\rm eff}$ is weaker by roughly an
order of magnitude, due to the different dimensionality of the operators. Fig.
5 shows the results for the remaining dim-6 operators, ${\cal
O}_{qX2,4,5}^{V,A}$. The constraints for ${\cal O}_{qX2}^{V}$ from the $B\to
K^{*}(\rho)XX$ exhibit similar behaviour to those for ${\cal O}_{qX2}^{T}$ in
Fig. 3. They are much weaker due to the higher dimensionality and result in
limits on $\Lambda_{\rm eff}$ of order a few tens of GeV. These are still much
larger than the $B$ meson mass implying that the LEFT framework is valid. On
the other hand, it may be difficult to interpret them within SMEFT.

For the $B^{+}\to K^{+}XX$ mode, the solid gray line in Figs. 3, 4, and 5
shows the sensitivity curve with the experimental efficiency included. It can
be seen that the sensitivity for $m=0$ is weaker by a factor of about 2 when
the efficiency is included, and then gradually decreases as $m$ increases.
Futhermore, the sensitivity in this case is limited to $m\lesssim
2\,$\mathrm{G}\mathrm{e}\mathrm{V}$$ by $4m^{2}\lesssim q^{2}_{\rm max}$, with
the maximum $q^{2}_{\rm max}\approx 16\,$\mathrm{G}\mathrm{e}\mathrm{V}$^{2}$
corresponding to the non-vanishing signal efficiency region in Belle II.

###  III.4 $B\to M+XX$ with vector DM $X$: scenario B

For these operators there is no issue with the $m\to 0$ limit, and in addition
all the form factors needed have been estimated in the literature before. The
non-vanishing amplitudes for the two processes $B(p)\to
P(k)X(k_{1})X^{*}(k_{2})$ and $B(p)\to V(k)X(k_{1})X^{*}(k_{2})$ take the
following form,

$\displaystyle i{\cal M}_{B\to PXX}^{B}$ $\displaystyle=$
$\displaystyle\epsilon^{*}_{\rho}(k_{1})\epsilon^{*}_{\sigma}(k_{2})\left\\{\left[2(k_{1}^{\sigma}k_{2}^{\rho}-k_{1}\cdot
k_{2}g^{\rho\sigma})\tilde{C}_{qX1}^{S,xb}+2\epsilon^{\mu\nu\rho\sigma}k_{1\mu}k_{2\nu}\tilde{C}_{qX2}^{S,xb}\right]\langle
P(k)|\bar{q}_{x}b|B(p)\rangle\right.$ (22a) $\displaystyle+$
$\displaystyle{i\over
2}\left[(k_{1}^{\alpha}k_{2}^{\beta}g^{\rho\sigma}+k_{1}\cdot
k_{2}g^{\alpha\rho}g^{\beta\sigma}-k_{1}^{\alpha}k_{2}^{\rho}g^{\beta\sigma}-k_{1}^{\sigma}k_{2}^{\beta}g^{\alpha\rho})\left(2g_{\mu\alpha}g_{\nu\beta}\tilde{C}_{qX1}^{T,xb}+\epsilon_{\mu\nu\alpha\beta}\tilde{C}_{qX2}^{T,xb}\right)\right]$
$\displaystyle\times$ $\displaystyle\left.\langle
P(k)|\bar{q}_{x}\sigma^{\mu\nu}b|B(p)\rangle\right\\},$ $\displaystyle i{\cal
M}_{B\to VXX}^{B}$ $\displaystyle=$
$\displaystyle\epsilon^{*}_{\rho}(k_{1})\epsilon^{*}_{\sigma}(k_{2})\left\\{\left[2(k_{1}^{\sigma}k_{2}^{\rho}-k_{1}\cdot
k_{2}g^{\rho\sigma})\tilde{C}_{qX1}^{P,xb}+2\epsilon^{\mu\nu\rho\sigma}k_{1\mu}k_{2\nu}\tilde{C}_{qX2}^{P,xb}\right]\langle
V(k)|\bar{q}_{x}i\gamma_{5}b|B(p)\rangle\right.$ (22b) $\displaystyle+$
$\displaystyle{i\over
2}\left[(k_{1}^{\alpha}k_{2}^{\beta}g^{\rho\sigma}+k_{1}\cdot
k_{2}g^{\alpha\rho}g^{\beta\sigma}-k_{1}^{\alpha}k_{2}^{\rho}g^{\beta\sigma}-k_{1}^{\sigma}k_{2}^{\beta}g^{\alpha\rho})\left(2g_{\mu\alpha}g_{\nu\beta}\tilde{C}_{qX1}^{T,xb}+\epsilon_{\mu\nu\alpha\beta}\tilde{C}_{qX2}^{T,xb}\right)\right]$
$\displaystyle\times$ $\displaystyle\left.\langle
V(k)|\bar{q}_{x}\sigma^{\mu\nu}b|B(p)\rangle\right\\}.$

From these amplitudes, the differential decay widths result in the following
compact forms,

$\displaystyle{d\Gamma_{B\to PXX}^{B}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{(m_{B}^{2}-m_{P}^{2})^{2}(s^{2}-4m^{2}s+6m^{4})\over
128\pi^{3}m_{B}^{3}(m_{b}-m_{q_{x}})^{2}}\lambda^{1\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{1\over
2}(m^{2},s)f_{0}^{2}\left|\tilde{C}_{qX1}^{S,xb}\right|^{2}$ (23)
$\displaystyle+$ $\displaystyle{(m_{B}^{2}-m_{P}^{2})^{2}s^{2}\over
128\pi^{3}m_{B}^{3}(m_{b}-m_{q_{x}})^{2}}\lambda^{1\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over
2}(m^{2},s)f_{0}^{2}\left|\tilde{C}_{qX2}^{S,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s(s+2m^{2})\over
1536\pi^{3}m_{B}^{3}(m_{B}+m_{P})^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{3\over
2}(m^{2},s)f_{T}^{2}\left|\tilde{C}_{qX1}^{T,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{s^{2}-2m^{2}s+4m^{4}\over
1536\pi^{3}m_{B}^{3}(m_{B}+m_{P})^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{P}^{2},s)\kappa^{1\over
2}(m^{2},s)f_{T}^{2}\left|\tilde{C}_{qX2}^{T,xb}\right|^{2},$
$\displaystyle{d\Gamma_{B\to VXX}^{B}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{s^{2}-4m^{2}s+6m^{4}\over
128\pi^{3}m_{B}^{3}(m_{b}+m_{q_{x}})^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)A_{0}^{2}\left|\tilde{C}_{qX1}^{P,xb}\right|^{2}$ (24)
$\displaystyle+$ $\displaystyle{s^{2}\over
128\pi^{3}m_{B}^{3}(m_{b}+m_{q_{x}})^{2}}\lambda^{3\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{3\over
2}(m^{2},s)A_{0}^{2}\left|\tilde{C}_{qX1}^{P,xb}\right|^{2}$ $\displaystyle+$
$\displaystyle{1\over 768\pi^{3}m_{B}^{3}s}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)\left\\{(s+2m^{2})(s-4m^{2})\lambda(m_{B}^{2},m_{V}^{2},s)T_{1}^{2}\right.$
$\displaystyle+$
$\displaystyle\left.(s^{2}-2m^{2}s+4m^{4})\left[(m_{B}^{2}-m_{V}^{2})^{2}T_{2}^{2}+{8m_{B}^{2}m_{V}^{2}s\over(m_{B}+m_{V})^{2}}T_{23}^{2}\right]\right\\}\left|\tilde{C}_{qX1}^{T,xb}\right|^{2}$
$\displaystyle+$ $\displaystyle{1\over 768\pi^{3}m_{B}^{3}s}\lambda^{1\over
2}(m_{B}^{2},m_{V}^{2},s)\kappa^{1\over
2}(m^{2},s)\left\\{(s^{2}-2m^{2}s+4m^{4})\lambda(m_{B}^{2},m_{V}^{2},s)T_{1}^{2}\right.$
$\displaystyle+$
$\displaystyle\left.(s+2m^{2})(s-4m^{2})\left[(m_{B}^{2}-m_{V}^{2})^{2}T_{2}^{2}+{8m_{B}^{2}m_{V}^{2}s\over(m_{B}+m_{V})^{2}}T_{23}^{2}\right]\right\\}\left|\tilde{C}_{qX2}^{T,xb}\right|^{2}.$

Similarly to the case of scalar DM, there is no interference between the
different operators due to their different parity and/or charge conjugation.
Interestingly, in the above expressions, the contributions from each pair of
operators with similar quark Lorentz structure (characterized by the same
superscript “S/P/T”) 101010Note that the two tensor currents are “similar” to
each other because $\bar{q}\sigma^{\mu\nu}\gamma_{5}q={i\over
2}\epsilon^{\mu\nu\rho\sigma}\bar{q}\sigma_{\rho\sigma}q$. become the same in
the limit of $m\to 0$.

Figure 6: Normalized differential decay width for $B\to K^{(*)+}XX$ from
different operators in the second scenario for vector DM. Left panel: $m=0$;
Right panel: $m=1\,\rm GeV$.

Figure 7: Constraints on the effective new physics scale for each operator as
a function of the DM mass $m$ from $B\to K(\pi)\not{E}$ channels.

Fig. 6 shows the normalized distributions for $B\to K^{+}XX$ (solid lines) and
$B\to K^{*+}XX$ (dashed lines) for different operators with $(sb)$ flavor
indices. In the left (right) panel, we set the vector DM mass to $m=0$ and $1$
GeV respectively. One can see in the left panel that the distributions are the
same for the two operators with the same superscript “S/P/T”. Even for
$m=1\,\rm GeV$ as shown in the right panel, this degeneracy still holds except
for a rather narrow range of $q^{2}$ value, which implies the
$d\tilde{\Gamma}/dq^{2}$ observable alone cannot be used to distinguish these
interactions.

To set numerical constraints we include the new physics scale in the Wilson
coefficients as $\tilde{C}_{qX1,2}^{S,P,T}\equiv\Lambda_{\rm eff}^{-3}$. Fig.
7 shows the bounds on the $m$-$\Lambda_{\rm eff}$ plane for each operator.
These bounds are comparable to those shown in Fig. 3 for operators with the
same quark current as they have the same dependence on $\Lambda_{\rm eff}$.
The effective scale $\Lambda_{\rm eff}$ which can be probed depends on the
value of DM mass, and is above a few hundreds of GeV for $m\lesssim
2\,$\mathrm{G}\mathrm{e}\mathrm{V}$$. Finally, the solid gray line in Fig. 7
shows the sensitivity after including the experimental efficiency for
$B^{+}\to K^{+}XX$. Similarly to the other two DM cases, the constraint is
weaker by a factor of 2 or so for $m=0$ when the efficiency is included, with
a sensitivity limit for $m$ around 2 GeV.

##  IV $K\to\pi$+DM+DM

Since the kaon decay process $K\rightarrow\pi\not{E}$ involves only the light
quarks $q=u,~{}d,~{}s$, its transition matrix element due to effective
interactions in LEFT can be evaluated by matching onto chiral perturbation
theory ($\chi$PT), which is the low energy effective field theory of QCD.
$\chi$PT is based on the fact that the QCD Lagrangian has the approximate
chiral symmetry $SU(3)_{L}\times SU(3)_{R}$ for the three light quarks which
is spontaneously broken by the quark condensate $\langle\bar{q}q\rangle$ to
the diagonal flavor $SU(3)_{V}$. The symmetry breakdown results in eight
pseudo-Nambu-Goldstone bosons (pNGBs), which are identified with the octet of
the lowest-lying pseudoscalars
$\pi^{\pm},\,\pi^{0},\,K^{\pm},\,K^{0},\,\overline{K^{0}},\,\eta$. In the
$\chi$PT formalism they are represented by the element in the coset space
$SU(3)_{L}\times SU(3)_{R}/SU(3)_{V}$ and take the matrix form,

$\displaystyle
U(x)=\exp\left(\frac{i\sqrt{2}\Pi(x)}{F_{0}}\right),\quad\Pi=\begin{pmatrix}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&\pi^{+}&K^{+}\\\
\pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-\sqrt{\frac{2}{3}}\eta\end{pmatrix},$ (25)

where $F_{0}$ is the pion decay constant in the chiral limit. Corresponding to
the chiral transformations of quarks $q_{L}\to Lq_{L}$ and $q_{R}\to Rq_{R}$,
$U$ transforms as $U\rightarrow LUR^{\dagger}$ with $L\in SU(3)_{L}$ and $R\in
SU(3)_{R}$.

The interactions of pNGBs with DM due to the effective operators in Eqs. (2-4)
can be realized through the external source method in which the global chiral
symmetry is promoted to a local one Gasser:1983yg ; Gasser:1984gg ;
Cata:2007ns . At the quark-gluon level, the QCD Lagrangian with all possible
external sources associated with quark bilinear currents is parameterized as
follows,

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\mathcal{L}_{\textrm{QCD}}+\overline{q}\hat{v}_{\mu}\gamma^{\mu}q+\overline{q}\hat{a}_{\mu}\gamma^{\mu}\gamma_{5}q-\overline{q}\hat{s}q+\overline{q}\hat{p}i\gamma_{5}q+\overline{q}\hat{t}_{\mu\nu}\sigma^{\mu\nu}q$
(26a) $\displaystyle=$
$\displaystyle\mathcal{L}_{\textrm{QCD}}+\overline{q_{L}}l_{\mu}\gamma^{\mu}q_{L}+\overline{q_{R}}r_{\mu}\gamma^{\mu}q_{R}-\left[\overline{q_{R}}(s+ip)q_{L}-\overline{q_{R}}t_{l}^{\mu\nu}\sigma_{\mu\nu}q_{L}+{\rm
h.c.}\right],$ (26b)

where $\mathcal{L}_{\textrm{QCD}}$ is the QCD Lagrangian for massless
$u,\,d,\,s$ quarks. The external sources,
$\hat{v}_{\mu}=\hat{v}_{\mu}^{\dagger}$,
$\hat{a}_{\mu}=\hat{a}_{\mu}^{\dagger}$, $\hat{s}=\hat{s}^{\dagger}$,
$\hat{p}=\hat{p}^{\dagger}$, $\hat{t}_{\mu\nu}=(\hat{t}_{\mu\nu})^{\dagger}$,
are $3\times 3$ Hermitian traceless matrices in flavor space. In the second
line, we rewrite the Lagrangian in terms of chiral quark fields to make the
chiral symmetry manifest. The external sources in the second line are related
to those in the first line by the relations,

$\displaystyle l_{\mu}=\hat{v}_{\mu}-\hat{a}_{\mu},\quad
r_{\mu}=\hat{v}_{\mu}+\hat{a}_{\mu},\quad\chi=2B(\hat{s}-i\hat{p}),\quad
t_{l}^{\mu\nu}=P_{L}^{\mu\nu\alpha\beta}\hat{t}_{\alpha\beta},\quad
t_{r}^{\mu\nu}=P_{R}^{\mu\nu\alpha\beta}\hat{t}_{\alpha\beta}=t_{l}^{\mu\nu\dagger},$
(27)

where the chiral projection for the tensor currents is defined by
$P_{R,L}^{\mu\nu\alpha\beta}={1\over
4}(g^{\mu\alpha}g^{\nu\beta}-g^{\mu\beta}g^{\nu\alpha}\pm
i\epsilon^{\mu\nu\alpha\beta})$ Cata:2007ns , with the property
$\hat{t}^{\mu\nu}=t_{l}^{\mu\nu}+t_{r}^{\mu\nu}$. Under chiral
transformations, $l_{\mu}\to Ll_{\mu}L^{\dagger}+iL\partial_{\mu}L^{\dagger}$,
$r_{\mu}\to Rr_{\mu}R^{\dagger}+iR\partial_{\mu}R^{\dagger}$, $\chi\to L\chi
R^{\dagger}$, and $t_{l}^{\mu\nu}\to Rt_{l}^{\mu\nu}L^{\dagger}$,
respectively. One should note that the tensor external sources have mass
dimension one in our convention. 111111Our convention for $U$ and $\chi$ are
equivalent to $U^{\dagger}$ and $\chi^{\dagger}$ in Cata:2007ns . The constant
$B$ is related to the quark condensate and $F_{0}$ by the relation
$B=-\langle\bar{q}q\rangle/(3F_{0}^{2})$. For numerical estimates, we use
$F_{0}=87\,\rm MeV$ Colangelo:2003hf  and $B\approx 2.8\,\rm GeV$.

By comparing the external sources in Eq. (26b) with the effective interactions
in LEFT in Eqs. (2-4) we see that the correspondence needed to calculate the
$K\to\pi$ transition with scalar DM is,

$\displaystyle(\hat{v}^{\mu})_{ds}=C_{q\phi}^{V,ds}\phi^{\dagger}i\overleftrightarrow{\partial_{\mu}}\phi,\quad(\hat{a}^{\mu})_{ds}=C_{q\phi}^{A,ds}\phi^{\dagger}i\overleftrightarrow{\partial_{\mu}}\phi,\quad(\hat{s})_{ds}=-C_{q\phi}^{S,ds}\phi^{\dagger}\phi,\quad(\hat{p})_{ds}=C_{q\phi}^{P,ds}\phi^{\dagger}\phi,$
(28)

and the corresponding Hermitian conjugates interchanging labels $s$ and $d$.
For the vector DM case, it is also easily to identify the relevant external
sources from the operators in Eqs. (3-4). For the second scenario of Eq. (4),
the non-vanishing sources are

$\displaystyle(\hat{s})_{ds}$ $\displaystyle=$
$\displaystyle-\tilde{C}_{qX1}^{S,ds}X_{\mu\nu}^{\dagger}X^{\mu\nu}-\tilde{C}_{qX2}^{S,ds}X_{\mu\nu}^{\dagger}\tilde{X}^{\mu\nu},$
(29a) $\displaystyle(\hat{p})_{ds}$ $\displaystyle=$
$\displaystyle\tilde{C}_{qX1}^{P,ds}X_{\mu\nu}^{\dagger}X^{\mu\nu}+\tilde{C}_{qX2}^{P,ds}X_{\mu\nu}^{\dagger}\tilde{X}^{\mu\nu},$
(29b) $\displaystyle(\hat{t}^{\mu\nu})_{ds}$ $\displaystyle=$
$\displaystyle{i\over
2}\tilde{C}_{qX1}^{T,ds}(X^{\dagger\mu}_{\,\,\,\,\rho}X^{\rho\nu}-X^{\dagger\nu}_{\,\,\,\,\rho}X^{\rho\mu})+{i\over
2}\tilde{C}_{qX2}^{T,ds}\epsilon^{\mu\nu\alpha\beta}X^{\dagger}_{\alpha\rho}X^{\rho}_{\,\beta},$
(29c)

and the corresponding Hermitian conjugates exchanging the $s$ and $d$ labels.

In $\chi$PT the vector and scalar sources first appear at order ${\cal
O}(p^{2})$ in the chiral power counting scheme Gasser:1983yg ; Gasser:1984gg

$\displaystyle\mathcal{L}^{(2)}_{\chi\rm PT}=\frac{F_{0}^{2}}{4}{\rm
Tr}\left[D_{\mu}U(D^{\mu}U)^{\dagger}\right]+\frac{F_{0}^{2}}{4}{\rm
Tr}\left[\chi U^{\dagger}+U\chi^{\dagger}\right],\quad
D_{\mu}U\equiv\partial_{\mu}U-il_{\mu}U+iUr_{\mu},$ (30)

whereas the tensor sources first appear at ${\cal O}(p^{4})$ Cata:2007ns

$\displaystyle\mathcal{L}^{(4)}_{\chi\rm PT}\ni i\Lambda_{2}{\rm
Tr}\left[t_{r}^{\mu\nu}(D_{\mu}U)^{\dagger}U(D_{\nu}U)^{\dagger}+t_{l}^{\mu\nu}D_{\mu}UU^{\dagger}D_{\nu}U\right],$
(31)

with $\Lambda_{2}$ being a new low energy constant with mass dimension one.
For our numerical estimates we follow Jiang:2012ir  and use
$|\Lambda_{2}|\approx 0.018\,\rm GeV$, noting this is comparable with naive
dimensional analysis estimates which find $\Lambda_{2}\sim{\Lambda_{\chi}\over
16\pi^{2}}\sim 0.008\,\rm GeV$ with the chiral symmetry breaking scale
$\Lambda_{\chi}\sim 1.2\,\rm GeV$. To find the relevant local interactions
mediating $K\to\pi\not{E}$ in question from the $\chi$PT formalism, it
suffices to expand the above Lagrangian in Eq. (30) and Eq. (31) to linear
order in the kaon and pion fields as well as each external source. This leads
to the following local interactions mediating $K\to\pi\not{E}$,

$\displaystyle\mathcal{L}_{K\to\pi}$ $\displaystyle=$
$\displaystyle-B\hat{s}_{sd}\pi^{-}K^{+}+{B\over\sqrt{2}}\hat{s}_{sd}\pi^{0}K^{0}$
$\displaystyle+$ $\displaystyle
i(\hat{v}_{\mu})_{sd}(\partial^{\mu}\pi^{-}K^{+}-\pi^{-}\partial^{\mu}K^{+})-{i\over\sqrt{2}}(\hat{v}_{\mu})_{sd}(\partial^{\mu}\pi^{0}K^{0}-\pi^{0}\partial^{\mu}K^{0})$
$\displaystyle-$ $\displaystyle i{\Lambda_{2}\over
F_{0}^{2}}\hat{t}^{\mu\nu}_{sd}(\partial_{\mu}\pi^{-}\partial_{\nu}K^{+}-\partial_{\nu}\pi^{-}\partial_{\mu}K^{+})+i{\Lambda_{2}\over\sqrt{2}F_{0}^{2}}\hat{t}^{\mu\nu}_{sd}(\partial_{\mu}\pi^{0}\partial_{\nu}K^{0}-\partial_{\nu}\pi^{0}\partial_{\mu}K^{0})+{\rm
h.c.}$ $\displaystyle\ni$ $\displaystyle-B\hat{s}_{sd}\pi^{-}K^{+}+{B\over
2}(\hat{s}_{sd}+\hat{s}_{ds})\pi^{0}K_{L}$ $\displaystyle+$ $\displaystyle
i(\hat{v}_{\mu})_{sd}(\partial^{\mu}\pi^{-}K^{+}-\pi^{-}\partial^{\mu}K^{+})-{i\over
2}[(\hat{v}_{\mu})_{sd}-(\hat{v}_{\mu})_{ds}](\partial^{\mu}\pi^{0}K_{L}-\pi^{0}\partial^{\mu}K_{L})$
$\displaystyle-$ $\displaystyle i{\Lambda_{2}\over
F_{0}^{2}}\hat{t}^{\mu\nu}_{sd}(\partial_{\mu}\pi^{-}\partial_{\nu}K^{+}-\partial_{\nu}\pi^{-}\partial_{\mu}K^{+})+i{\Lambda_{2}\over
2F_{0}^{2}}(\hat{t}^{\mu\nu}_{sd}-\hat{t}^{\mu\nu}_{ds})(\partial_{\mu}\pi^{0}\partial_{\nu}K_{L}-\partial_{\nu}\pi^{0}\partial_{\mu}K_{L}).$
(32b)

In the second line we ignore CP violation in kaon mixing to write
$K^{0}(\bar{K}^{0})\approx{1\over\sqrt{2}}(K_{L}\pm K_{S})$. It can be seen in
the above effective Lagrangian that the neutral mode $K_{L}\to\pi^{0}\not{E}$
can be obtained from the charged mode by replacing the relevant Wilson
coefficients $C_{i}^{sd}$ with $(C_{i}^{sd}\pm
C_{i}^{ds})/2\sim\Re[C_{i}^{ds}](\Im[C_{i}^{ds}])$. The plus sign (real)
applies to the scalar current while the minus sign (imaginary) applies to the
vector and tensor currents. The above Lagrangian leads to the following form
factors

$\displaystyle\langle\pi^{+}(k)|{\bar{s}}d|K^{+}(p)\rangle$
$\displaystyle\simeq$ $\displaystyle B,$ (33a)
$\displaystyle\langle\pi^{+}(k)|{\bar{s}}\gamma^{\mu}d|K^{+}(p)\rangle$
$\displaystyle\simeq$ $\displaystyle\left(p+k\right)^{\mu},$ (33b)
$\displaystyle\langle\pi^{+}(k)|{\bar{s}}\sigma^{\mu\nu}d|K^{+}(p)\rangle$
$\displaystyle\simeq$ $\displaystyle i{\Lambda_{2}\over
F_{0}^{2}}(p^{\mu}k^{\nu}-p^{\nu}k^{\mu}),$ (33c)

These matrix elements have exactly the same Lorentz structure as those for
$B\to P$ transitions given in Eq. (5) in the limit of $q^{2}\to 0$ and
noticing that $f_{+}(0)=f_{0}(0)$. Thus, the differential decay width for
$K\to\pi\not{E}$ can be obtained directly from the result for $B$ decay with
suitable replacements of variables. The differential decay width for the
charged mode and scalar DM is then

$\displaystyle{d\Gamma_{K^{+}\to\pi^{+}\phi\phi}\over dq^{2}}$
$\displaystyle=$ $\displaystyle{B^{2}\over 256\pi^{3}m_{K}^{3}}\lambda^{1\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{1\over
2}(m^{2},s)\left|C_{q\phi}^{S,ds}\right|^{2}$ (34) $\displaystyle+$
$\displaystyle{1\over 768\pi^{3}m_{K}^{3}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over
2}(m^{2},s)\left|C_{q\phi}^{V,ds}\right|^{2}.$

Figure 8: Constraints on the effective new physics scale for the two operators
${\cal O}_{a\phi}^{S,ds}$ and ${\cal O}_{a\phi}^{V,ds}$ involving $d,s$ quarks
as a function of the DM mass $m$ from $K\to\pi\not{E}$ channel.

Figure 9: Constraints on the effective new physics scale for the vector DM
operators with $d,s$ quarks as a function of the DM mass $m$ from
$K\to\pi\not{E}$ channel.

For the two vector DM scenarios, using Eqs. (2-4) and Eq. (33), the
differential decay widths are

$\displaystyle{d\Gamma_{K^{+}\to\pi^{+}XX}^{A}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{B^{2}(s^{2}-4m^{2}s+12m^{4})\over
1024\pi^{3}m_{K}^{3}m^{4}}\lambda^{1\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{1\over
2}(m^{2},s)\left|C_{qX}^{S,ds}\right|^{2}$ (35) $\displaystyle+$
$\displaystyle{\Lambda_{2}^{2}s(s+4m^{2})\over
12288\pi^{3}F_{0}^{4}m_{K}^{3}m^{4}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over
2}(m^{2},s)\left|C_{qX1}^{T,ds}\right|^{2}$ $\displaystyle+$
$\displaystyle{\Lambda_{2}^{2}(s+2m^{2})\over
3072\pi^{3}F_{0}^{4}m_{K}^{3}m^{2}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{1\over
2}(m^{2},s)\left|C_{qX2}^{T,ds}\right|^{2}$ $\displaystyle+$
$\displaystyle{s\over 3072\pi^{3}m_{K}^{3}m^{4}}\lambda^{1\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left[3(s-4m^{2})(m_{K}^{2}-m_{\pi}^{2})^{2}+4m^{2}\lambda(m_{K}^{2},m_{\pi}^{2},s)\right]\left|C_{qX2}^{V,ds}\right|^{2}$
$\displaystyle+$ $\displaystyle{1\over
768\pi^{3}m_{K}^{3}m^{2}}\lambda^{1\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over 2}(m^{2},s)$ $\displaystyle\times$
$\displaystyle\left[6m^{2}(m_{K}^{2}-m_{\pi}^{2})^{2}+(s-4m^{2})\lambda(m_{K}^{2},m_{\pi}^{2},s)\right]\left|C_{qX3}^{V,ds}\right|^{2}$
$\displaystyle+$ $\displaystyle{s^{2}-4m^{2}s+12m^{4}\over
3072\pi^{3}m_{K}^{3}m^{4}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over
2}(m^{2},s)\left|C_{qX4}^{V,ds}\right|^{2}$ $\displaystyle+$
$\displaystyle{s(s+4m^{2})\over 3072\pi^{3}m_{K}^{3}m^{4}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over
2}(m^{2},s)\left|C_{qX5}^{V,ds}\right|^{2}$ $\displaystyle+$
$\displaystyle{s+2m^{2}\over 768\pi^{3}m_{K}^{3}m^{2}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{1\over
2}(m^{2},s)\left|C_{qX6}^{V,ds}\right|^{2}+\cdots,$
$\displaystyle{d\Gamma_{K^{+}\to\pi^{+}XX}^{B}\over dq^{2}}$ $\displaystyle=$
$\displaystyle{B^{2}(s^{2}-4m^{2}s+6m^{4})\over
128\pi^{3}m_{K}^{3}}\lambda^{1\over 2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{1\over
2}(m^{2},s)\left|\tilde{C}_{qX1}^{S,ds}\right|^{2}$ (36) $\displaystyle+$
$\displaystyle{B^{2}s^{2}\over 128\pi^{3}m_{K}^{3}}\lambda^{1\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over
2}(m^{2},s)\left|\tilde{C}_{qX2}^{S,ds}\right|^{2}$ $\displaystyle+$
$\displaystyle{\Lambda_{2}^{2}s(s+2m^{2})\over
6144\pi^{3}F_{0}^{4}m_{K}^{3}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{3\over
2}(m^{2},s)\left|\tilde{C}_{qX1}^{T,ds}\right|^{2}$ $\displaystyle+$
$\displaystyle{\Lambda_{2}^{2}(s^{2}-2m^{2}s+4m^{4})\over
6144\pi^{3}F_{0}^{4}m_{K}^{3}}\lambda^{3\over
2}(m_{K}^{2},m_{\pi}^{2},s)\kappa^{1\over
2}(m^{2},s)\left|\tilde{C}_{qX2}^{T,ds}\right|^{2}.$

The corresponding results for the neutral kaon mode $K_{L}\to\pi^{0}\not{E}$
are obtained from these ones replacing $C_{i}^{S,ds}$ by their real parts
$\Re[C_{i}^{S,ds}]$ for scalar quark currents and $C_{i}^{V(T),ds}$ by their
imaginary parts $\Im[C_{i}^{V(T),ds}]$ for vector and tensor quark currents,
respectively.

Fig. 8 shows the bounds obtained on $\Lambda_{\rm eff}$ from $K\to\pi\phi\phi$
for scalar DM. Kaon decays only cover the low mass region, $m\lesssim 180$
MeV, but the constraints on $\Lambda_{\rm eff}$ in this region are much
stronger than the corresponding ones in $B$ meson decay shown in Fig. 2. The
difference, a factor of ${\cal O}(10^{4})$ (${\cal O}(10)$) for scalar
(vector) current operators, is due both to the much longer kaon lifetime and
the much stronger experimental bounds on kaon modes.

Fig. 9 shows the constraints for the vector DM cases: the upper four panels
for scenario A and the lower two panels for scenario B. Unlike the case of
scalar DM, the constraints for vector DM case from kaon and $B$ meson decays
are similar as seen in Figs. 3, 4, 5, and 7. This can be understood from
dimensional arguments: the decay width from vector DM operators scales
approximately as $\Gamma\sim{m_{K,B}^{7}/\Lambda_{\rm
eff}^{6}}({m_{K,B}^{9}/\Lambda_{\rm eff}^{8}})$ for operators ${\cal
O}^{S}_{qX},{\cal O}^{T}_{qX1,2},{\cal O}^{V}_{qX3,6},\tilde{\cal
O}^{S}_{qX1,2},\tilde{\cal O}^{T}_{qX1,2},$ (${\cal O}^{V}_{qX2,4,5}$). This
large enhancement of $(m_{B}/m_{K})^{7(9)}$ for $B$ mesons compensates for its
shorter lifetime and weaker experimental bounds, resulting in comparable
constraints on $\Lambda_{\rm eff}$ (of course, for different flavor indices).

Cases | Operators | Lower bound on $\Lambda_{\rm eff}\,\rm[GeV]$
---|---|---
|  | $(sb)$-flavor | $(db)$-flavor | $(ds)$-flavor
|  | $m=0$ | $m=2\,\rm GeV$ | $m=0$ | $m=2\,\rm GeV$ | $m=0$ | $m=150\,\rm MeV$
scalar | ${\cal O}_{q\phi}^{S}\,[\Lambda_{\rm eff}^{-1}]$ | $2.4\times 10^{7}$ | $7.8\times 10^{6}$ | $2.6\times 10^{7}$ | $1.1\times 10^{7}$ | $1.6\times 10^{11}$ | $3.5\times 10^{10}$
${\cal O}_{q\phi}^{P}\,[\Lambda_{\rm eff}^{-1}]$ | $1.9\times 10^{7}$ | $1.7\times 10^{6}$ | $1.2\times 10^{7}$ | $1.6\times 10^{6}$ | — | —
${\cal O}_{q\phi}^{V}\,[\Lambda_{\rm eff}^{-2}]$ | $7.0\times 10^{3}$ | $2.2\times 10^{3}$ | $7.0\times 10^{3}$ | $2.9\times 10^{3}$ | $7.4\times 10^{4}$ | $1.3\times 10^{4}$
${\cal O}_{q\phi}^{A}\,[\Lambda_{\rm eff}^{-2}]$ | $8.0\times 10^{3}$ | $1.9\times 10^{3}$ | $6.2\times 10^{3}$ | $1.7\times 10^{3}$ | — | —
vector A | ${\cal O}_{qX}^{S}\,[m^{2}\Lambda_{\rm eff}^{-3}]$ | $540$ | $400$ | $590$ | $470$ | $1.7\times 10^{3}$ | $1.2\times 10^{3}$
${\cal O}_{qX}^{P}\,[m^{2}\Lambda_{\rm eff}^{-3}]$ | $430$ | $230$ | $390$ | $230$ | — | —
${\cal O}_{qX1}^{T}\,[m^{2}\Lambda_{\rm eff}^{-3}]$ | $410$ | $310$ | $410$ | $300$ | $520$ | $250$
${\cal O}_{qX2}^{T}\,[m^{2}\Lambda_{\rm eff}^{-3}]$ | $470$ | $290$ | $400$ | $340$ | — | $330$
${\cal O}_{qX2}^{V}\,[m^{2}\Lambda_{\rm eff}^{-4}]$ | $160$ | $94$ | $170$ | $110$ | $140$ | $73$
${\cal O}_{qX3}^{V}\,[m\Lambda_{\rm eff}^{-3}]$ | $540$ | $340$ | $550$ | $410$ | $1.0\times 10^{3}$ | $610$
${\cal O}_{qX4}^{V}\,[m^{2}\Lambda_{\rm eff}^{-4}]$ | $130$ | $80$ | $130$ | $93$ | $100$ | $53$
${\cal O}_{qX5}^{V}\,[m^{2}\Lambda_{\rm eff}^{-4}]$ | $130$ | $90$ | $130$ | $100$ | $100$ | $60$
${\cal O}_{qX6}^{V}\,[m\Lambda_{\rm eff}^{-3}]$ | $540$ | $400$ | $550$ | $450$ | $1.0\times 10^{3}$ | $590$
${\cal O}_{qX2}^{A}\,[m^{2}\Lambda_{\rm eff}^{-4}]$ | $140$ | $74$ | $120$ | $72$ | — | —
${\cal O}_{qX3}^{A}\,[m\Lambda_{\rm eff}^{-3}]$ | $580$ | $200$ | $500$ | $200$ | — | —
${\cal O}_{qX4}^{A}\,[m^{2}\Lambda_{\rm eff}^{-4}]$ | $140$ | $72$ | $120$ | $70$ | — | —
${\cal O}_{qX5}^{A}\,[m^{2}\Lambda_{\rm eff}^{-4}]$ | $140$ | $81$ | $120$ | $79$ | — | —
${\cal O}_{qX6}^{A}\,[m\Lambda_{\rm eff}^{-3}]$ | $580$ | $390$ | $500$ | $360$ | — | —
vector B | $\tilde{\cal O}_{qX1}^{S}\,[\Lambda_{\rm eff}^{-3}]$ | $760$ | $520$ | $830$ | $630$ | $2.4\times 10^{3}$ | $1.5\times 10^{3}$
$\tilde{\cal O}_{qX2}^{S}\,[\Lambda_{\rm eff}^{-3}]$ | $760$ | $450$ | $830$ | $570$ | $2.4\times 10^{3}$ | $1.3\times 10^{3}$
$\tilde{\cal O}_{qX1}^{P}\,[\Lambda_{\rm eff}^{-3}]$ | $610$ | $300$ | $550$ | $300$ | — | —
$\tilde{\cal O}_{qX2}^{P}\,[\Lambda_{\rm eff}^{-3}]$ | $610$ | $230$ | $550$ | $240$ | — | —
$\tilde{\cal O}_{qX1}^{T}\,[\Lambda_{\rm eff}^{-3}]$ | $560$ | $310$ | $480$ | $320$ | $580$ | $270$
$\tilde{\cal O}_{qX2}^{T}\,[\Lambda_{\rm eff}^{-3}]$ | $560$ | $310$ | $480$ | $360$ | $580$ | $340$

Table 3: The strongest bounds on the effective scale $\Lambda_{\rm eff}$
associated with all FCNC operators with two representative DM masses: $m=0$
and $m=2\,\rm GeV$ (or $m=150\,\rm MeV$ for $(ds)$\- flavor). For the
$(ds)$-flavor, all the constraints in the table come from the charged channel
$K^{+}\to\pi^{+}\not{E}$.

In practice, the decay-in-flight search of $K^{+}\to\pi^{+}+\not{E}$ by NA62
experiment NA62:2020fhy ; NA62:2018ctf  has only two signal regions
corresponding to $0<q^{2}<0.01\,$\mathrm{G}\mathrm{e}\mathrm{V}$^{2}$ (region
1) and
$0.026\,$\mathrm{G}\mathrm{e}\mathrm{V}$^{2}<q^{2}<0.068\,$\mathrm{G}\mathrm{e}\mathrm{V}$^{2}$
(region 2), with both signal regions also being constrained in the pion
momentum by
$15\,$\mathrm{G}\mathrm{e}\mathrm{V}$<|\boldsymbol{p}_{\pi}|<35\,$\mathrm{G}\mathrm{e}\mathrm{V}$$
(equivalently, the pion energy by
$15\,$\mathrm{G}\mathrm{e}\mathrm{V}$<E_{\pi}<35\,$\mathrm{G}\mathrm{e}\mathrm{V}$$).
When restricting the phase space to these kinematic windows, the corresponding
sensitivity bounds on the parameter space shift from the solid color lines
into the solid gray lines in Figs. 8 and 9. The bound weakens by factors of a
few relative to what could be obtained from the entire phase space. The
sensitivity drops quickly for heavier DM masses and vanishes around $m\approx
130\,$\mathrm{M}\mathrm{e}\mathrm{V}$$, the cutoff value being determined by
the largest $q^{2}$, $q^{2}_{\rm max}\approx
0.068\,$\mathrm{G}\mathrm{e}\mathrm{V}$^{2}$, in the NA62 signal region.

Finally, in Tab. 3, we summarize the strongest constraints on $\Lambda_{\rm
eff}$ for all FCNC interactions with two representative DM masses: $m=0$ and
$m=2\,\rm GeV$ (or $m=150\,\rm MeV$ for kaon decays). In the second column, we
show the scaling behavior of the Wilson coefficients employed to obtain the
bounds. One should note that all results we present correspond to the case of
complex DM fields. For the case of real DM fields, the operators ${\cal
O}_{q\phi}^{V,A}$ in Eq. (2), ${\cal O}_{qX1,2}^{T},{\cal O}_{qX4,5,6}^{V,A}$
in Eq. (3), and $\tilde{\cal O}_{qX1,2}^{T}$ in Eq. (4) do not exist. For the
remaining ones, the bounds on $\Lambda_{\rm eff}$ will be enhanced by a factor
of $2^{1/2n}$ with $n$ being the power of $\Lambda_{\rm eff}^{-1}$ in the
corresponding Wilson coefficient, shown in the second column of Tab. 3.

##  V Summary and Conclusions

In this paper we have carried out a systematic study of possible flavor
changing neutral current $B$ and $K$ meson decays with a pair of light scalar
or vector invisible particles in the final state using the effective field
theory approach. This completes the existing studies of $B(K)\to M\not{E}$
transitions where the missing energy is attributed to a pair of new invisible
particles in the context of effective field theory. The case of two invisible
fermions was studied before, and we have now addressed the cases of two
invisible scalars or vectors. This study is particularly relevant when new
symmetries forbid the appearance of single DM particles.

We first constructed the local quark-DM interactions relevant to these
processes in the low energy effective field theory framework at leading order.
We tabulated results for invisible scalar, fermion, or vector particles
completing and correcting existing lists in the literature. We then used the
effective interactions to consider $B\to(K,\pi,K^{*},\rho)$+DM+DM transitions
using the form factor formalism, followed by $K\to\pi$+DM+DM transitions using
chiral perturbation theory. We describe the different characteristic features
of each operator in the differential decay rate, which could be exploited in
future detailed experimental searches to differentiate between various
possibilities of new physics. Finally, with the help of the most recent
experimental results on these modes, we set constraints on all the relevant
effective operators. The sensitivity to the new physics scale strongly depends
on the operator structure as well as the DM mass. In the massless limit, for
the $B$ meson decay with scalar (vector) DM, we find that the current
experimental data can probe the effective new physics scale up to ${\cal
O}(10^{7})({\cal O}(10^{3}))$ GeV for some operators, while for the kaon
decays ${\cal O}(10^{11})({\cal O}(10^{3}))$ GeV can be reached.

In two cases, $B^{+}\to K^{+}\not{E}$ and $K^{+}\to\pi^{+}\not{E}$, we
considered the effect that experimental efficiency affects the theoretical
constraints. For the former we relied on the current Belle II sensitivity to
bins of different $q^{2}$ and for the latter on the signal window in NA62.

###### Acknowledgements.

XDM would like to thank Shao-Zhou Jiang for clarifying their result for the
LEC $\Lambda_{2}$ associated with one $p^{4}$ chiral Lagrangian term for the
tensor external source in Jiang:2012ir . This work was supported in part by
the National Natural Science Foundation of China (Nos. 12090064, 11975149,
11735010), the Shanghai Pujiang Program (20PJ1407800), Chinese Academy of
Sciences Center for Excellence in Particle Physics (CCEPP), and Key Laboratory
for Particle Physics, Astrophysics and Cosmology, Ministry of Education, and
Shanghai Key Laboratory for Particle Physics and Cosmology. XGH was also
supported in part by the MOST (Grant No. MOST 106- 2112-M-002-003-MY3 ). GV
and XG also was supported in part by the Australian Government through the
Australian Research Council.

##  Appendix A Differential decay rate

For reference we detail here the kinematics relevant for the modes discussed
in this paper. For the three-body decay, $B(p)\to M(k)+{\rm DM}(k_{1})+{\rm
DM}(k_{2})$, there are two independent Mandelstam variables describing the
kinematics that are denoted as $s\equiv(p-k)^{2}=(k_{1}+k_{2})^{2}$ and
$t\equiv(p-k_{1})^{2}=(k_{2}+k)^{2}$. Then the scalar product of any pair of
four vectors can be expressed in terms of the two Mandelstam variables and
masses as follows,

$\displaystyle p\cdot k_{1}={m_{B}^{2}+m_{1}^{2}-t\over 2},\quad p\cdot
k_{2}={s+t-m_{1}^{2}-m_{M}^{2}\over 2},\quad p\cdot
k={m_{B}^{2}+m_{M}^{2}-s\over 2},$ (37) $\displaystyle k_{1}\cdot
k_{2}={s-m_{1}^{2}-m_{2}^{2}\over 2},\quad k_{1}\cdot
k={m_{B}^{2}+m_{2}^{2}-s-t\over 2},\quad k_{2}\cdot
k={t-m_{2}^{2}-m_{M}^{2}\over 2},$ (38)

where $m_{B}$ is the initial $B$ meson mass while $m_{M}$ for the mass of
final state meson; $m_{1}=m_{2}\equiv m$ is the mass of the DM particle. For
the process to happen, the largest possible DM mass is restricted by
kinematics to be $m\leq(M_{B}-m_{M})/2$.

The differential decay width can be expressed as

$\displaystyle{d\Gamma\over dq^{2}}={1\over S}{1\over
256\pi^{3}m_{B}^{3}}\int_{t_{-}}^{t_{+}}dt\overline{|{\cal M}|^{2}},$ (39)

where $S$ is a possible symmetry factor for identical DM particles at the
final state, and $\overline{|{\cal M}|^{2}}$ is the spin-averaged squared
matrix element for the relevant decay process. The integration domain for $t$
is

$\displaystyle
t_{\pm}=(E_{2}^{*}+E_{3}^{*})^{2}-\left(\sqrt{E_{2}^{*2}-m_{2}^{2}}\mp\sqrt{E_{3}^{*2}-m_{M}^{2}}\right)^{2},\,E_{2}^{*}={s-m_{1}^{2}+m_{2}^{2}\over
2\sqrt{s}},\,E_{3}^{*}={m_{B}^{2}-s-m_{M}^{2}\over 2\sqrt{s}}.$ (40)

To extract the constraints we take one operator at a time, therefore ignoring
interference between different operators. As the DM are not SM particles,
there is never interference with the SM. If these interactions exist, they
thus contribute additively to the SM and can be directly constrained by the
“room for new physics” of the last column in Tab. 2. The branching ratio can
be written as a sum of numerical coefficients times the squares of the
couplings of the new operators, schematically

$\displaystyle{\cal B}_{\rm NP}={1\over\Gamma_{B}^{\rm tot}}\int
dq^{2}{d\Gamma\over dq^{2}}\equiv\sum_{i}\hat{B}_{i}(m)|C_{i}|^{2},$ (41)

where the generic range of $q^{2}$ goes from $(m_{1}+m_{2})^{2}$ to
$(m_{B}-m_{M})^{2}$.

We allow for two modifications to the integration range to accommodate
reported details of existing experiments. For $B^{+}\to K^{+}\not{E}$ we
multiply the differential branching ratio with the efficiency as a function of
$q^{2}$ reported in Fig. 3 (supplementary material) of Belle-II:2021rof . For
$K^{+}\to\pi^{+}\not{E}$, we limit the integration region to the two signal
regions of NA62 NA62:2020fhy ; NA62:2018ctf

$\displaystyle 15\leq|\boldsymbol{p}_{\pi}|\leq 35\,{\rm GeV}\text{ for the
pion momentum in the NA62 rest frame,}$ $\displaystyle 0\leq q^{2}\leq
0.1\,{\rm GeV}^{2}~{}{\rm or}~{}0.026\leq q^{2}\leq 0.068\,{\rm GeV}^{2}\text{
for the two signal regions.}$ (42)

In this case we work on the NA62 lab frame, with the kaon momentum
$|\boldsymbol{p}_{K}|=75\,$\mathrm{G}\mathrm{e}\mathrm{V}$$.

Requiring that the NP contribution of each operator to the branching ratio
does not exceed the value ${\cal B}^{\rm UL}$ given in Tab. 2, we set the
constraints

$\displaystyle|C_{i}|^{2}\leq\left({{\cal B}^{\rm
UL}\over\hat{B}_{i}(m)}\right).$ (43)

To interpret this as a bound on new physics we then write
$C_{i}\equiv\Lambda_{\rm eff}^{-n}$ (with the power $n$ depending on the
dimension of the corresponding operator), leading to

$\displaystyle\Lambda_{\rm eff}(m)\geq\left({\hat{B}_{i}(m)\over{\cal B}^{\rm
UL}}\right)^{1\over 2n}.$ (44)

##  Appendix B Detailed analysis for the operators with a vector DM

Here we describe in detail how to obtain the operators with a vector DM in Eq.
(3) for scenario A. First, we can always choose the operators to be self-
conjugate for the flavor diagonal case as given in Eq. (3), and we denote the
Hermitian and anti-Hermitian combination of vector DM fields as,
$S^{\mu\nu}\equiv X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}$ and
$A^{\mu\nu}\equiv X^{\mu\dagger}X^{\nu}-X^{\nu\dagger}X^{\mu}$. It can be seen
that they are automatically symmetric and anti-symmetric in their two Lorentz
indices, respectively. For a combination of the vector quark current
$\overline{q}\gamma_{\mu}q$ with two DM fields $X^{\dagger}_{\alpha}X_{\beta}$
to form dim-6 operators by attaching an additional derivative, we have the
following possibilities

$\displaystyle(\overline{q}\gamma_{\mu}q)X^{\dagger}_{\nu}X^{\nu}$
$\displaystyle\overset{i\partial_{\mu}}{\Rightarrow}$
$\displaystyle(\overline{q}\gamma^{\mu}q)X^{\dagger}_{\nu}i\overleftrightarrow{\partial_{\mu}}X^{\nu},$
(45a) $\displaystyle(\overline{q}\gamma_{\mu}q)S^{\mu\nu}$
$\displaystyle\overset{i\partial_{\nu}}{\Rightarrow}$
$\displaystyle(\overline{q}\gamma_{\mu}q)\partial_{\nu}S^{\mu\nu},\,(\overline{q}\gamma_{\mu}i\overleftrightarrow{D_{\nu}}q)S^{\mu\nu},\,$
(45b) $\displaystyle(\overline{q}\gamma_{\mu}q)A^{\mu\nu}$
$\displaystyle\overset{i\partial_{\nu}}{\Rightarrow}$
$\displaystyle(\overline{q}\gamma_{\mu}q)i\partial_{\nu}A^{\mu\nu},\,\boxed{(\overline{q}\gamma_{\mu}\overleftrightarrow{D_{\nu}}q)A^{\mu\nu}},$
(45c)
$\displaystyle(\overline{q}\gamma_{\mu}q)X^{\dagger}_{\rho}X_{\sigma}\epsilon^{\mu\nu\rho\sigma}$
$\displaystyle\overset{\partial_{\nu}}{\Rightarrow}$
$\displaystyle(\overline{q}\gamma_{\mu}q)(X^{\dagger}_{\rho}\overleftrightarrow{\partial_{\nu}}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},\,(\overline{q}\gamma_{\mu}q)\partial_{\nu}A_{\rho\sigma}\epsilon^{\mu\nu\rho\sigma},\,\boxed{(\overline{q}\gamma_{\mu}\overleftrightarrow{D_{\nu}}q)A_{\rho\sigma}\epsilon^{\mu\nu\rho\sigma}},$
(45d)

where we have used IBP and the on-shell condition $\partial_{\mu}X^{\mu}=0$.
For the quark axial-vector current $\overline{q}\gamma_{\mu}\gamma_{5}q$, a
similar operator realization can be obtained by replacing the vector gamma
matrix $\gamma_{\mu}$ in the quark current by $\gamma_{\mu}\gamma_{5}$.

Above, we have written 8 possible operators. However, the two operators in a
“$\Box$” are redundant and they can be transformed into others appearing in
our basis given in Eq. (3). Using the Dirac gamma matrix identities (DIs),

$\displaystyle\gamma^{\mu}\gamma^{\nu}=g^{\mu\nu}-i\sigma^{\mu\nu},\quad\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}=g^{\mu\nu}\gamma^{\rho}+g^{\nu\rho}\gamma^{\mu}-g^{\mu\rho}\gamma^{\nu}+i\epsilon^{\mu\nu\rho\sigma}\gamma_{\sigma}\gamma_{5},\quad\sigma^{\mu\nu}\gamma_{5}={i\over
2}\epsilon^{\mu\nu\rho\sigma}\sigma_{\rho\sigma},$ (46)

the two operators can be manipulated as follows,

$\displaystyle
2(\overline{q_{a}}\gamma_{\mu}\overleftrightarrow{D_{\nu}}q_{b})A^{\mu\nu}$
$\displaystyle=$
$\displaystyle\left[\overline{q_{a}}(\gamma_{\mu}\gamma_{\nu}\not{D}-\overleftarrow{\not{D}}\gamma_{\nu}\gamma_{\mu})q_{b}+\overline{q_{a}}(\gamma_{\mu}\not{D}\gamma_{\nu}-\gamma_{\nu}\overleftarrow{\not{D}}\gamma_{\mu})q_{b}\right]A^{\mu\nu}$
(47a) $\displaystyle\overset{\rm DI}{=}$
$\displaystyle-\left[\overline{q_{a}}(\sigma_{\mu\nu}i\not{D}+i\overleftarrow{\not{D}}\sigma_{\mu\nu})q_{b}+\overline{q_{a}}(\gamma_{\mu}\not{D}\gamma_{\nu}+\gamma_{\mu}\overleftarrow{\not{D}}\gamma_{\nu})q_{b}\right]A^{\mu\nu}$
$\displaystyle\overset{\rm EoM}{=}$
$\displaystyle(m_{a}-m_{b})(\overline{q_{a}}\sigma_{\mu\nu}q_{b})A^{\mu\nu}-\partial^{\alpha}(\overline{q_{a}}i\epsilon_{\mu\alpha\nu\beta}\gamma^{\beta}\gamma_{5}q_{b})A^{\mu\nu}$
$\displaystyle\overset{\rm IBP}{=}$
$\displaystyle(m_{a}-m_{b})(\overline{q_{a}}\sigma_{\mu\nu}q_{b})A^{\mu\nu}+(\overline{q_{a}}\gamma^{\mu}\gamma_{5}q_{b})i\partial^{\nu}A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma},$
$\displaystyle
2(\overline{q_{a}}\gamma^{\mu}\overleftrightarrow{D^{\nu}}q_{b})A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma}$
$\displaystyle\overset{\rm DI}{=}$
$\displaystyle-\left[\overline{q_{a}}(\sigma^{\mu\nu}i\not{D}+i\overleftarrow{\not{D}}\sigma^{\mu\nu})q_{b}+\overline{q_{a}}(\gamma^{\mu}\not{D}\gamma^{\nu}+\gamma^{\mu}\overleftarrow{\not{D}}\gamma^{\nu})q_{b}\right]A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma}$
(47b) $\displaystyle\overset{\rm EoM}{=}$
$\displaystyle(m_{a}-m_{b})(\overline{q_{a}}\sigma^{\mu\nu}q_{b})A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma}-\partial_{\alpha}(\overline{q_{a}}i\epsilon^{\mu\alpha\nu\beta}\gamma_{\beta}\gamma_{5}q_{b})A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma}$
$\displaystyle\overset{\rm IBP}{=}$ $\displaystyle
2i(m_{b}-m_{a})(\overline{q_{a}}\sigma^{\mu\nu}\gamma_{5}q_{b})A_{\mu\nu}+(\overline{q_{a}}\gamma_{\beta}\gamma_{5}q_{b})i\partial_{\alpha}A^{\rho\sigma}2(\delta^{\alpha}_{\rho}\delta^{\beta}_{\sigma}-\delta^{\beta}_{\rho}\delta^{\alpha}_{\sigma})$
$\displaystyle=$ $\displaystyle
2i(m_{b}-m_{a})(\overline{q_{a}}\sigma_{\mu\nu}\gamma_{5}q_{b})A^{\mu\nu}-4(\overline{q_{a}}\gamma_{\mu}\gamma_{5}q_{b})i\partial_{\nu}A^{\mu\nu}.$

In the reduction of the second operator, we used the identity
$\epsilon^{\mu\nu\alpha\beta}\epsilon_{\mu\nu\rho\sigma}=-2(\delta^{\alpha}_{\rho}\delta^{\beta}_{\sigma}-\delta^{\beta}_{\rho}\delta^{\alpha}_{\sigma})$.
Since the final 4 operators are already in our basis, we conclude that the two
operators in a “$\Box$” are redundant. Similarly, for the axial-vector quark
current, the two corresponding operators are redundant and can be transformed
into those in our basis as follows,

$\displaystyle
2(\overline{q_{a}}\gamma_{\mu}\overleftrightarrow{D_{\nu}}\gamma_{5}q_{b})A^{\mu\nu}$
$\displaystyle=$
$\displaystyle(m_{a}+m_{b})(\overline{q_{a}}\sigma_{\mu\nu}\gamma_{5}q_{b})A^{\mu\nu}+(\overline{q_{a}}\gamma^{\mu}q_{b})i\partial^{\nu}A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma},$
(48a)
$\displaystyle(\overline{q_{a}}\gamma^{\mu}\overleftrightarrow{D^{\nu}}\gamma_{5}q_{b})A^{\rho\sigma}\epsilon_{\mu\nu\rho\sigma}$
$\displaystyle=$
$\displaystyle-i(m_{a}+m_{b})(\overline{q_{a}}\sigma_{\mu\nu}q_{b})A^{\mu\nu}-2(\overline{q_{a}}\gamma_{\mu}q_{b})i\partial_{\nu}A^{\mu\nu}.$
(48b)

In conclusion, there are 6 independent operators for each quark current and
they are listed in Eq. (3) after normalization and (anti-)symmetrization.

##  Appendix C Lepton-DM interaction in LEFT

For completeness we list here the independent operators involving a lepton
current and two dark sector particles. Denoting the charged leptons as
$\ell\in\\{e,\mu,\tau\\}$, the charged lepton-DM interactions can be directly
obtained from the quark-DM interactions given in section II by exchanging the
quark flavor label $q$ by the lepton label $\ell$. Following the conventions
for quark-DM interactions in section II, they are:

Fermion case:

$\displaystyle{\cal O}_{\ell\chi 1}^{S}$
$\displaystyle=(\overline{\ell}\ell)(\overline{\chi}\chi),$
$\displaystyle{\cal O}_{\ell\chi 2}^{S}$
$\displaystyle=(\overline{\ell}\ell)(\overline{\chi}i\gamma_{5}\chi),$ (49a)
$\displaystyle{\cal O}_{\ell\chi 1}^{P}$
$\displaystyle=(\overline{\ell}i\gamma_{5}\ell)(\overline{\chi}\chi),$
$\displaystyle{\cal O}_{\ell\chi 2}^{P}$
$\displaystyle=(\overline{\ell}\gamma_{5}\ell)(\overline{\chi}\gamma_{5}\chi),$
(49b) $\displaystyle{\cal O}_{\ell\chi 1}^{V}$
$\displaystyle=(\overline{\ell}\gamma_{\mu}\ell)(\overline{\chi}\gamma_{\mu}\chi),\,(\times)$
$\displaystyle{\cal O}_{\ell\chi 2}^{V}$
$\displaystyle=(\overline{\ell}\gamma_{\mu}\ell)(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi),$
(49c) $\displaystyle{\cal O}_{\ell\chi 1}^{A}$
$\displaystyle=(\overline{\ell}\gamma_{\mu}\gamma_{5}\ell)(\overline{\chi}\gamma_{\mu}\chi),\,(\times)$
$\displaystyle{\cal O}_{\ell\chi 2}^{A}$
$\displaystyle=(\overline{\ell}\gamma_{\mu}\gamma_{5}\ell)(\overline{\chi}\gamma_{\mu}\gamma_{5}\chi),$
(49d) $\displaystyle{\cal O}_{\ell\chi 1}^{T}$
$\displaystyle=(\overline{\ell}\sigma_{\mu\nu}\ell)(\overline{\chi}\sigma_{\mu\nu}\chi),\,(\times)$
$\displaystyle{\cal O}_{\ell\chi 2}^{T}$
$\displaystyle=(\overline{\ell}\sigma_{\mu\nu}\ell)(\overline{\chi}\sigma_{\mu\nu}\gamma_{5}\chi),\,(\times)$
(49e)

Scalar case:

$\displaystyle{\cal O}_{\ell\phi}^{S}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\ell)(\phi^{\dagger}\phi),$ (50a)
$\displaystyle{\cal O}_{\ell\phi}^{P}$ $\displaystyle=$
$\displaystyle(\overline{\ell}i\gamma_{5}\ell)(\phi^{\dagger}\phi),$ (50b)
$\displaystyle{\cal O}_{\ell\phi}^{V}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma^{\mu}\ell)(\phi^{\dagger}i\overleftrightarrow{\partial_{\mu}}\phi),\,(\times)$
(50c) $\displaystyle{\cal O}_{\ell\phi}^{A}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma^{\mu}\gamma_{5}\ell)(\phi^{\dagger}i\overleftrightarrow{\partial_{\mu}}\phi),\,(\times).$
(50d)

Vector case A:

$\displaystyle{\cal O}_{\ell X}^{S}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\ell)(X_{\mu}^{\dagger}X^{\mu}),$ (51a)
$\displaystyle{\cal O}_{\ell X}^{P}$ $\displaystyle=$
$\displaystyle(\overline{\ell}i\gamma_{5}\ell)(X_{\mu}^{\dagger}X^{\mu}),$
(51b) $\displaystyle{\cal O}_{\ell X1}^{T}$ $\displaystyle=$
$\displaystyle{i\over
2}(\overline{\ell}\sigma^{\mu\nu}\ell)(X_{\mu}^{\dagger}X_{\nu}-X_{\nu}^{\dagger}X_{\mu}),\,(\times)$
(51c) $\displaystyle{\cal O}_{\ell X2}^{T}$ $\displaystyle=$
$\displaystyle{1\over
2}(\overline{\ell}\sigma^{\mu\nu}\gamma_{5}\ell)(X_{\mu}^{\dagger}X_{\nu}-X_{\nu}^{\dagger}X_{\mu}),\,(\times)$
(51d) $\displaystyle{\cal O}_{\ell X1}^{V}$ $\displaystyle=$
$\displaystyle{1\over
2}[\overline{\ell}\gamma_{(\mu}i\overleftrightarrow{D_{\nu)}}\ell](X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(51e) $\displaystyle{\cal O}_{\ell X2}^{V}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\ell)\partial_{\nu}(X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(51f) $\displaystyle{\cal O}_{\ell X3}^{V}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\ell)(X_{\rho}^{\dagger}\overleftrightarrow{\partial_{\nu}}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},$
(51g) $\displaystyle{\cal O}_{\ell X4}^{V}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma^{\mu}\ell)(X_{\nu}^{\dagger}i\overleftrightarrow{\partial_{\mu}}X^{\nu}),\,(\times)$
(51h) $\displaystyle{\cal O}_{\ell X5}^{V}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\ell)i\partial_{\nu}(X^{\mu\dagger}X^{\nu}-X^{\nu\dagger}X^{\mu}),\,(\times)$
(51i) $\displaystyle{\cal O}_{\ell X6}^{V}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\ell)i\partial_{\nu}(X^{\dagger}_{\rho}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},\,(\times)$
(51j) $\displaystyle{\cal O}_{\ell X1}^{A}$ $\displaystyle=$
$\displaystyle{1\over
2}[\overline{\ell}\gamma_{(\mu}\gamma_{5}i\overleftrightarrow{D_{\nu)}}\ell](X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(51k) $\displaystyle{\cal O}_{\ell X2}^{A}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\gamma_{5}\ell)\partial_{\nu}(X^{\mu\dagger}X^{\nu}+X^{\nu\dagger}X^{\mu}),$
(51l) $\displaystyle{\cal O}_{\ell X3}^{A}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\gamma_{5}\ell)(X_{\rho}^{\dagger}\overleftrightarrow{\partial_{\nu}}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},$
(51m) $\displaystyle{\cal O}_{\ell X4}^{A}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma^{\mu}\gamma_{5}\ell)(X_{\nu}^{\dagger}i\overleftrightarrow{\partial_{\mu}}X^{\nu}),\,(\times)$
(51n) $\displaystyle{\cal O}_{\ell X5}^{A}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\gamma_{5}\ell)i\partial_{\nu}(X^{\mu\dagger}X^{\nu}-X^{\nu\dagger}X^{\mu}),\,(\times)$
(51o) $\displaystyle{\cal O}_{\ell X6}^{A}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\gamma_{\mu}\gamma_{5}\ell)i\partial_{\nu}(X^{\dagger}_{\rho}X_{\sigma})\epsilon^{\mu\nu\rho\sigma},\,(\times)$
(51p)

Vector case B:

$\displaystyle\tilde{\cal O}_{\ell X1}^{S}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\ell)X_{\mu\nu}^{\dagger}X^{\mu\nu},$ (52a)
$\displaystyle\tilde{\cal O}_{\ell X2}^{S}$ $\displaystyle=$
$\displaystyle(\overline{\ell}\ell)X_{\mu\nu}^{\dagger}\tilde{X}^{\mu\nu},$
(52b) $\displaystyle\tilde{\cal O}_{\ell X1}^{P}$ $\displaystyle=$
$\displaystyle(\overline{\ell}i\gamma_{5}\ell)X_{\mu\nu}^{\dagger}X^{\mu\nu},$
(52c) $\displaystyle\tilde{\cal O}_{\ell X2}^{P}$ $\displaystyle=$
$\displaystyle(\overline{\ell}i\gamma_{5}\ell)X_{\mu\nu}^{\dagger}\tilde{X}^{\mu\nu},$
(52d) $\displaystyle\tilde{\cal O}_{q\ell X1}^{T}$ $\displaystyle=$
$\displaystyle{i\over
2}(\overline{q}\sigma^{\mu\nu}\ell)(X^{\dagger}_{\mu\rho}X^{\rho}_{\,\nu}-X^{\dagger}_{\nu\rho}X^{\rho}_{\,\mu}),\,(\times)$
(52e) $\displaystyle\tilde{\cal O}_{\ell X2}^{T}$ $\displaystyle=$
$\displaystyle{1\over
2}(\overline{\ell}\sigma^{\mu\nu}\gamma_{5}\ell)(X^{\dagger}_{\mu\rho}X^{\rho}_{\,\nu}-X^{\dagger}_{\nu\rho}X^{\rho}_{\,\mu}),\,(\times)$
(52f)

In addition to the usually considered electron-DM scattering process for
direct DM detection, the above lepton-DM interactions can induce lepton flavor
violating transitions, $\ell_{j}\to\ell_{i}$+DM+DM. The charged leptons can be
replaced by SM neutrinos leading to exotic interactions between neutrinos and
DM particles. We defer a study of these possibilities to a future publication.

##  Appendix D Form factors for the $B\to P(V)$ transitions

In Tabs. 4 and 5 we collect the fitted parameters for the form factor
parameterizations with $q^{2}\neq 0$ as given in Eq. (6) and Eq. (13) relevant
to $B\to P$ and $B\to V$ transitions respectively.

Form factor | $r_{1}$ | $r_{2}$ | $m_{\rm fit}^{2}($\mathrm{G}\mathrm{e}\mathrm{V}^{2}$)$ | $m_{R}($\mathrm{G}\mathrm{e}\mathrm{V}$)$
---|---|---|---|---
$f_{+}^{\pi}$ | 0.744 | -0.486 | 40.73 | 5.32
$f_{0}^{\pi}$ | 0 | 0.258 | 33.81 | —
$f_{T}^{\pi}$ | 1.387 | -1.134 | 32.22 | 5.32
$f_{+}^{K}$ | 0.162 | 0.173 | — | 5.41
$f_{0}^{K}$ | 0 | 0.330 | 37.46 | —
$f_{T}^{K}$ | 0.161 | 0.198 | — | 5.41
Table 4: Parameters appearing in the form factors of the $B\to K(\pi)$ transitions Ball:2004ye . $f_{T}$ is a scale-dependent quantity and the value is given at $\mu=4.8\,$\mathrm{G}\mathrm{e}\mathrm{V}$$. $F_{i}$ | $B\to K^{*}$ |  $m_{R,i}^{b\to s}/$GeV | $B\to\rho$ |  $m_{R,i}^{b\to d}/$GeV
---|---|---|---|---
$\alpha_{0}^{A_{0}}$ | $0.36\pm 0.05$ |  | $0.36\pm 0.04$ |
$\alpha_{1}^{A_{0}}$ | $-1.04\pm 0.27$ | $5.366$ | $-0.83\pm 0.20$ | $5.279$
$\alpha_{2}^{A_{0}}$ | $1.12\pm 1.35$ |  | $1.33\pm 1.05$ |
$\alpha_{0}^{A_{1}}$ | $0.27\pm 0.03$ |  | $0.26\pm 0.03$ |
$\alpha_{1}^{A_{1}}$ | $0.30\pm 0.19$ | $5.829$ | $0.39\pm 0.14$ | $5.724$
$\alpha_{2}^{A_{1}}$ | $-0.11\pm 0.48$ |  | $0.16\pm 0.41$ |
$\alpha_{0}^{A_{12}}$ | $0.26\pm 0.03$ |  | $0.30\pm 0.03$ |
$\alpha_{1}^{A_{12}}$ | $0.60\pm 0.20$ | $5.829$ | $0.76\pm 0.20$ | $5.724$
$\alpha_{2}^{A_{12}}$ | $0.12\pm 0.84$ |  | $0.46\pm 0.76$ |
$\alpha_{0}^{V_{0}}$ | $0.34\pm 0.04$ |  | $0.33\pm 0.03$ |
$\alpha_{1}^{V_{0}}$ | $-1.05\pm 0.24$ | $5.415$ | $-0.86\pm 0.18$ | $5.325$
$\alpha_{2}^{V_{0}}$ | $2.37\pm 1.39$ |  | $1.80\pm 0.97$ |
$\alpha_{0}^{T_{1}}$ | $0.28\pm 0.03$ |  | $0.27\pm 0.03$ |
$\alpha_{1}^{T_{1}}$ | $-0.89\pm 0.19$ | $5.415$ | $-0.74\pm 0.14$ | $5.325$
$\alpha_{2}^{T_{1}}$ | $1.95\pm 1.10$ |  | $1.45\pm 0.77$ |
$\alpha_{0}^{T_{2}}$ | $0.28\pm 0.03$ |  | $0.27\pm 0.03$ |
$\alpha_{1}^{T_{2}}$ | $0.40\pm 0.18$ | $5.829$ | $0.47\pm 0.13$ | $5.724$
$\alpha_{2}^{T_{2}}$ | $0.36\pm 0.51$ |  | $0.58\pm 0.46$ |
$\alpha_{0}^{T_{23}}$ | $0.67\pm 0.08$ |  | $0.75\pm 0.08$ |
$\alpha_{1}^{T_{23}}$ | $1.48\pm 0.49$ | $5.829$ | $1.90\pm 0.43$ | $5.724$
$\alpha_{2}^{T_{23}}$ | $1.92\pm 1.96$ |  | $2.93\pm 1.81$ |
Table 5: Parameters in the form factors of the $B\to\rho(K^{*})$ processes
with $k_{\rm max}=2$ Bharucha:2015bzk .

##  Appendix E Specific renormalizable models to illustrate a possible origin
of the LEFT operators

Here we illustrate with two simple renormalizable models a possible origin for
the LEFT operators. First for the case of scalar DM generating the operators
with a scalar mediator. Then for the case A of vector DM illustrating a
possible origin for the additional mass factors that we argued should
accompany the operators in Eq. (3) when $X_{\mu}$ is regarded as a dark sector
gauge boson.

In these examples we will employ three dark sector fields: a light real scalar
$\phi$ and a real vector gauge particle $X_{\mu}$ from a $U(1)_{X}$ gauge
group, and a second scalar $\Delta$ that gives mass to $X$. These three
particles are singlets under the SM gauge group. To generate the FCNC in the
quark sector we introduce two Higgs doublets $H_{1}$ and $H_{2}$. Under the
full gauge group $SU(3)_{\rm c}\times SU(2)_{L}\times U(1)_{Y}(U_{X}(1))$,
these new scalars have the following charge assignments,

$\displaystyle H_{i}(1,2,1/2)(0)=\left(\begin{array}[]{c}h^{+}_{i}\\\
{v_{i}+h_{i}+iI_{i}\over\sqrt{2}}\end{array}\right),\quad\Delta(1,1,0)(1)={v_{\Delta}+h_{\Delta}+iI_{\Delta}\over\sqrt{2}}.$
(55)

Scalar case: The relevant Lagrangian is given by

$\displaystyle{\cal L}_{\rm
scalar}\ni\lambda_{ij}^{H\phi}H_{i}^{\dagger}H_{j}\phi\phi+\left(\bar{Q}_{L}Y_{1}H_{1}D_{R}+\bar{Q}_{L}Y_{2}H_{2}D_{R}+{\rm
h.c.}\right),$ (56)

where $Q_{L}$ and $D_{R}$ are the SM left-handed quark doublet and right-
handed down-type quark singlet, respectively. After spontaneous symmetry
breaking and rotating back to the physical states, $h_{1}$ and $h_{2}$ will
mix giving rise to a SM Higgs ($h^{m}_{1}$) and a heavy Higgs ($h^{m}_{2}$),
$h_{i}=\alpha_{ij}h_{j}^{m}$. Integrating out the physical Higgs bosons, we
generate the operators ${\cal O}_{q\phi}^{S}$ and ${\cal O}_{q\phi}^{P}$.

Vector case: The relevant interactions are given by

$\displaystyle{\cal L}_{\rm
vector}\ni(D^{\mu}\Delta)^{\dagger}(D_{\mu}\Delta)+\lambda^{H\Delta}_{ij}H^{\dagger}_{i}H_{j}\Delta^{\dagger}\Delta+\left(\bar{Q}_{L}Y_{1}H_{1}D_{R}+\bar{Q}_{L}Y_{2}H_{2}D_{R}+{\rm
h.c.}\right),$ (57)

where $D_{\mu}=\partial_{\mu}-ig_{X}X_{\mu}$. After $\Delta$ develops a vev,
the $\Delta$ kinetic term will induce a mass for $X$ given by
$m=g_{X}v_{\Delta}$ and a $h_{\Delta}$-$X$-$X$ vertex with a coupling
$g_{X}^{2}v_{\Delta}$. The Higgs quartic interaction leads to a mixing of
$h_{\Delta}$ and $h_{i}$ with a mixing parameter ${1\over
2}\lambda_{ik}^{H\Delta}v_{k}v_{\Delta}$. Then integrating out $h_{\Delta}$
and $h_{i}$ we can obtain the operator ${\cal O}_{qX}^{S}$ and ${\cal
O}_{qX}^{P}$, with coefficients proportional to $m^{2}$.

One should note that if only one Higgs doublet is introduced, the
diagonalization of the mass matrix also diagonalizes the Yukawa matrix $Y$ and
there are no tree-level FCNC interactions. Introducing two Higgs doublets,
where $Y_{1}$ and $Y_{2}$ cannot be diagonalized simultaneously, generates
tree-level FCNC interactions. Since exchange of one of these Higgs bosons
would generate undesirable FCNC interactions purely within the quark sector,
this simple model should be considered as an existence proof. A realistic
model would have to introduce many of the features that resolve this issue in
two Higgs doublet models.

## References

  * (1) S. L. Glashow, J. Iliopoulos and L. Maiani, “Weak Interactions with Lepton-Hadron Symmetry,” Phys. Rev. D 2, 1285-1292 (1970)
  * (2) L. Roszkowski, E. M. Sessolo and S. Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81, no.6, 066201 (2018) [arXiv:1707.06277 [hep-ph]].
  * (3) S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo and L. Vittorio, “Closing the window on WIMP Dark Matter,” Eur. Phys. J. C 82, no.1, 31 (2022) [arXiv:2107.09688 [hep-ph]].
  * (4) S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40, 223-226 (1978)
  * (5) F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40, 279-282 (1978)
  * (6) B. Dasgupta and J. Kopp, “Sterile Neutrinos,” Phys. Rept. 928, 1-63 (2021) [arXiv:2106.05913 [hep-ph]].
  * (7) M. Fabbrichesi, E. Gabrielli and G. Lanfranchi, “The Dark Photon,” [arXiv:2005.01515 [hep-ph]].
  * (8) C. O. Dib, J. C. Helo, V. E. Lyubovitskij, N. A. Neill, A. Soffer and Z. S. Wang, “Probing $R$-parity violation in $B$-meson decays to a baryon and a light neutralino,” [arXiv:2208.06421 [hep-ph]].
  * (9) J. F. Kamenik and C. Smith, “FCNC portals to the dark sector,” JHEP 03, 090 (2012) [arXiv:1111.6402 [hep-ph]].
  * (10) G. Li, T. Wang, J. B. Zhang and G. L. Wang, “The light invisible boson in FCNC decays of B and $B_{c}$ mesons,” Eur. Phys. J. C 81, no.6, 564 (2021) [arXiv:2103.12921 [hep-ph]].
  * (11) X. G. He, X. D. Ma, J. Tandean and G. Valencia, “Evading the Grossman-Nir bound with $\Delta I=3/2$ new physics,” JHEP 08, no.08, 034 (2020) [arXiv:2005.02942 [hep-ph]].
  * (12) G. Li, T. Wang, Y. Jiang, J. B. Zhang and G. L. Wang, “Spin-$1\over 2$ invisible particles in heavy meson decays,” Phys. Rev. D 102, no.9, 095019 (2020) [arXiv:2004.10942 [hep-ph]].
  * (13) T. Felkl, S. L. Li and M. A. Schmidt, “A tale of invisibility: constraints on new physics in b → s$\nu$$\nu$,” JHEP 12, 118 (2021) [arXiv:2111.04327 [hep-ph]].
  * (14) T. Li, X. D. Ma and M. A. Schmidt, “Implication of $K\to\pi\nu\bar{\nu}$ for generic neutrino interactions in effective field theories,” Phys. Rev. D 101, no.5, 055019 (2020) [arXiv:1912.10433 [hep-ph]].
  * (15) F. F. Deppisch, K. Fridell and J. Harz, “Constraining lepton number violating interactions in rare kaon decays,” JHEP 12, 186 (2020) [arXiv:2009.04494 [hep-ph]].
  * (16) X. G. He and G. Valencia, “RK(*)$\nu$ and non-standard neutrino interactions,” Phys. Lett. B 821, 136607 (2021) [arXiv:2108.05033 [hep-ph]].
  * (17) R. L. Workman et al. [Particle Data Group], “Review of Particle Physics,” PTEP 2022, 083C01 (2022)
  * (18) J. Grygier et al. [Belle], “Search for $\bm{B\to h\nu\bar{\nu}}$ decays with semileptonic tagging at Belle,” Phys. Rev. D 96, no.9, 091101 (2017) [arXiv:1702.03224 [hep-ex]].
  * (19) O. Lutz et al. [Belle], “Search for $B\to h^{(*)}\nu\bar{\nu}$ with the full Belle $\Upsilon(4S)$ data sample,” Phys. Rev. D 87, no.11, 111103 (2013) [arXiv:1303.3719 [hep-ex]].
  * (20) C. Hambrock, A. Khodjamirian and A. Rusov, “Hadronic effects and observables in $B\to\pi\ell^{+}\ell^{-}$ decay at large recoil,” Phys. Rev. D 92, no.7, 074020 (2015) [arXiv:1506.07760 [hep-ph]].
  * (21) D. M. Straub, “flavio: a Python package for flavour and precision phenomenology in the Standard Model and beyond,” [arXiv:1810.08132 [hep-ph]].
  * (22) A. J. Buras, D. Buttazzo, J. Girrbach-Noe and R. Knegjens, “${K}^{+}\to{\pi}^{+}\nu\overline{\nu}$ and ${K}_{L}\to{\pi}^{0}\nu\overline{\nu}$ in the Standard Model: status and perspectives,” JHEP 11, 033 (2015) [arXiv:1503.02693 [hep-ph]].
  * (23) A. V. Artamonov et al. [E949], “New measurement of the $K^{+}\to\pi^{+}\nu\bar{\nu}$ branching ratio,” Phys. Rev. Lett. 101, 191802 (2008) [arXiv:0808.2459 [hep-ex]].
  * (24) A. V. Artamonov et al. [BNL-E949], “Study of the decay $K^{+}\to\pi^{+}\nu\bar{\nu}$ in the momentum region $140<P_{\pi}<199$ MeV/c,” Phys. Rev. D 79, 092004 (2009) [arXiv:0903.0030 [hep-ex]].
  * (25) E. Cortina Gil et al. [NA62], “An investigation of the very rare ${K}^{+}\to{\pi}^{+}\nu\overline{\nu}$ decay,” JHEP 11, 042 (2020) [arXiv:2007.08218 [hep-ex]].
  * (26) E. Cortina Gil et al. [NA62], “Measurement of the very rare K+→${\pi}^{+}\nu\overline{\nu}$ decay,” JHEP 06, 093 (2021) [arXiv:2103.15389 [hep-ex]].
  * (27) B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, “Dimension-Six Terms in the Standard Model Lagrangian,” JHEP 10, 085 (2010) [arXiv:1008.4884 [hep-ph]].
  * (28) L. Lehman, “Extending the Standard Model Effective Field Theory with the Complete Set of Dimension-7 Operators,” Phys. Rev. D 90, no.12, 125023 (2014) [arXiv:1410.4193 [hep-ph]].
  * (29) Y. Liao and X. D. Ma, “Renormalization Group Evolution of Dimension-seven Baryon- and Lepton-number-violating Operators,” JHEP 11, 043 (2016) [arXiv:1607.07309 [hep-ph]].
  * (30) E. E. Jenkins, A. V. Manohar and P. Stoffer, “Low-Energy Effective Field Theory below the Electroweak Scale: Operators and Matching,” JHEP 03, 016 (2018) [arXiv:1709.04486 [hep-ph]].
  * (31) Y. Liao, X. D. Ma and Q. Y. Wang, “Extending low energy effective field theory with a complete set of dimension-7 operators,” JHEP 08, 162 (2020) [arXiv:2005.08013 [hep-ph]].
  * (32) B. V. Lehmann and S. Profumo, “Cosmology and prospects for sub-MeV dark matter in electron recoil experiments,” Phys. Rev. D 102, no.2, 023038 (2020) [arXiv:2002.07809 [hep-ph]].
  * (33) A. Badin and A. A. Petrov, “Searching for light Dark Matter in heavy meson decays,” Phys. Rev. D 82, 034005 (2010) [arXiv:1005.1277 [hep-ph]].
  * (34) J. Kumar and D. Marfatia, “Matrix element analyses of dark matter scattering and annihilation,” Phys. Rev. D 88, no.1, 014035 (2013) [arXiv:1305.1611 [hep-ph]].
  * (35) B. Henning, X. Lu, T. Melia and H. Murayama, “Operator bases, $S$-matrices, and their partition functions,” JHEP 10, 199 (2017) [arXiv:1706.08520 [hep-th]].
  * (36) J. C. Criado, “BasisGen: automatic generation of operator bases,” Eur. Phys. J. C 79, no.3, 256 (2019) [arXiv:1901.03501 [hep-ph]].
  * (37) M. Williams, C. P. Burgess, A. Maharana and F. Quevedo, “New Constraints (and Motivations) for Abelian Gauge Bosons in the MeV-TeV Mass Range,” JHEP 08, 106 (2011) [arXiv:1103.4556 [hep-ph]].
  * (38) J. Brod, A. Gootjes-Dreesbach, M. Tammaro and J. Zupan, “Effective Field Theory for Dark Matter Direct Detection up to Dimension Seven,” JHEP 10, 065 (2018) [arXiv:1710.10218 [hep-ph]].
  * (39) J. C. Criado, A. Djouadi, M. Perez-Victoria and J. Santiago, “A complete effective field theory for dark matter,” JHEP 07, 081 (2021) [arXiv:2104.14443 [hep-ph]].
  * (40) 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, 086 (2022) [arXiv:2202.06968 [hep-ph]].
  * (41) N. Gubernari, A. Kokulu and D. van Dyk, “$B\to P$ and $B\to V$ Form Factors from $B$-Meson Light-Cone Sum Rules beyond Leading Twist,” JHEP 01, 150 (2019) [arXiv:1811.00983 [hep-ph]].
  * (42) P. Ball and R. Zwicky, “New results on $B\to\pi,K,\eta$ decay formfactors from light-cone sum rules,” Phys. Rev. D 71, 014015 (2005) [arXiv:hep-ph/0406232 [hep-ph]].
  * (43) A. Bharucha, D. M. Straub and R. Zwicky, “$B\to V\ell^{+}\ell^{-}$ in the Standard Model from light-cone sum rules,” JHEP 08, 098 (2016) [arXiv:1503.05534 [hep-ph]].
  * (44) V. Shtabovenko, R. Mertig and F. Orellana, “New Developments in FeynCalc 9.0,” Comput. Phys. Commun. 207, 432-444 (2016) [arXiv:1601.01167 [hep-ph]].
  * (45) F. Abudinén et al. [Belle-II], “Search for B+→K+$\nu$$\nu$¯ Decays Using an Inclusive Tagging Method at Belle II,” Phys. Rev. Lett. 127, no.18, 181802 (2021) [arXiv:2104.12624 [hep-ex]].
  * (46) J. Gasser and H. Leutwyler, “Chiral Perturbation Theory to One Loop,” Annals Phys. 158, 142 (1984)
  * (47) J. Gasser and H. Leutwyler, “Chiral Perturbation Theory: Expansions in the Mass of the Strange Quark,” Nucl. Phys. B 250, 465-516 (1985)
  * (48) O. Cata and V. Mateu, “Chiral perturbation theory with tensor sources,” JHEP 09, 078 (2007) [arXiv:0705.2948 [hep-ph]].
  * (49) G. Colangelo and S. Durr, “The Pion mass in finite volume,” Eur. Phys. J. C 33, 543-553 (2004) [arXiv:hep-lat/0311023 [hep-lat]].
  * (50) S. Z. Jiang, Q. Wang and Y. Zhang, “Computation of the $p^{6}$ order low-energy constants with tensor sources,” Phys. Rev. D 87, no.9, 094014 (2013) [arXiv:1203.0712 [hep-ph]].
  * (51) E. Cortina Gil et al. [NA62], “First search for $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ using the decay-in-flight technique,” Phys. Lett. B 791, 156-166 (2019) [arXiv:1811.08508 [hep-ex]].