FCNC portals to the dark sector

The most general basis of operators parametrizing a low-scale departure from the SM particle content is constructed. The SM gauge invariance is enforced, and operators of lowest dimensions are retained separately for a new light neutral particle of spin 0, 1/2, 1, and 3/2. The basis is further decomposed into couplings to the SM Higgs/gauge fields, to pairs of quark/lepton fields, and to baryon/lepton number violating combinations of fermion fields. This basis is then used to systematically investigate the discovery potential of the rare FCNC decays of the K and B mesons with missing energy in the final state. The most sensitive decay modes in the s to d, b to d, and b to s sectors are identified and compared for each type of couplings to the new invisible state.


Introduction
Though extremely successful, the Standard Model (SM) is not expected to be valid up to arbitrary high energies. It certainly needs to be amended at the Planck scale, with the advent of quantum gravity. But then, some new dynamics should show up already at a much lower scale, though still above the electroweak scale, to avoid hierarchy problems. With this picture in mind, it seems natural to assume that all the New Physics (NP) degrees of freedom are heavier than the known SM particles, and decouple at low energy. In a fully model-independent way, the impact on the SM of any NP model can then be embedded into nonrenormalizable effective operators, under the provision that these involve only the SM fields and satisfy the SM gauge symmetries. The full classification of these operators has been achieved many years ago, by Buchmuller and Wyler [1] for the baryon-and lepton-number conserving operators, and by Weinberg [2] for those violating these global symmetries.
Still, whether the SM particles are the only dynamical degrees of freedom within the electroweak energy range is far from established. Not only is the existence of new light particles not excluded, since they would evade direct detection when sufficiently weakly interacting, but their presence could even be welcome. Indeed, many NP models are built upon some spontaneously broken symmetries, and do often have remnants at lowenergy in the form of massless or very light Goldstone bosons. A well-known example is the axion [3], introduced to cure the strong CP problem of the SM. More crucially, there are now very strong indications that the universe is filled with dark matter, so there should be at least one new electrically neutral colorless particle, possibly lighter than the electroweak scale. Once opening that door, it is not such a drastic step to imagine a whole dark sector, i.e. a full-fledged set of darkly interacting dark particles only loosely connected to our own visible sector. Further, it should be stressed that adjoining a dark sector to the SM is always possible, does not need to be directly related to dark matter (so one would rather speak of a hidden sector), and is actually quite generic in supersymmetric models. For a recent review, including further physical motivations from string theory or extra dimensional settings, see e.g. Ref. [4].
In the present work, our main goal is to construct the lowest-dimensional effective interactions parametrizing a low-scale departure from the SM particle content. Specifically, we assume that there is a new particle of spin 0, 1/2, 1, or 3/2, neutral under the SM gauge group SU (3) C ⊗ SU (2) L ⊗ U (1) Y , and write down the gauge invariant operators coupling this particle to SM fields. These effective interactions are not yet included and thus complement the NP operator basis of Refs. [1,2]. To our knowledge, such a complete basis has never been presented, though parts of it already appeared in the literature. In particular, those effective couplings between SM and dark particles which are renormalizable, sometimes called portals, have already been investigated [5].
Our second goal is to constrain the effective operators. Since the new state is assumed neutral under the SM gauge group, it looks like a natural dark matter candidate. However, the viability of this hypothesis would require constraining its mass, couplings, and lifetime, and this in general requires more inputs about the dark sector dynamics. Here, we refrain from doing so and rather concentrate exclusively on the quark FCNC transitions s → dX, b → dX, and b → sX, where X collectively denotes any final state made of dark particles. Our focus on these modes, instead of for example leptonic observables or collider signals, is motivated on one hand by their extreme sensitivities to NP (as detailed in the next section), and on the other, by the next generation of experiments currently under construction. Indeed, rare K decays are the main targets of the NA62 (CERN) and K0TO (J-Parc) experiments, while rare B decays could be accessed at the Super-B (Italy) and Belle II (KEK) facilities. So, in the present work, we further assume that the dark particle is light enough to be directly produced in K and/or B meson decays, and sufficiently long-lived to escape detection in flavor factories. For all practical matters, this new particle is invisible, and would show up as missing energy in FCNC-induced rare K and B decays (for recent works along this line, see e.g. Refs. [6,7]).
Our analysis is organized according to the dark particle spin. To keep the discussion focused on the operator basis, we rely on extensive appendices to cover the issues of hadronic matrix elements and differential rates. So, before entering the discussion, the next section summarizes the main features of the observables considered here, i.e. the rare K and B decays. Then, given our focus on FCNC processes, a flavor-based classification of the dark operators is described in the following section, on which we rely throughout the paper.  Table 1: Naive reach, in terms of scales Λ and as a function of the effective operator dimension n, of the rare FCNC-induced K and B decays, as estimated from Eq. (1) with the CKM values of Eq. (6).

Rare FCNC decays
The FCNC-induced decay modes are very suppressed in the SM, where the missing energy is carried away by a νν pair (see Fig. 1). So, even relatively small NP contributions could be evidenced. Specifically, to set the stage and get an idea of the sensitivity of the rare K and B decays, imagine that a NP operator of dimension n contributes to d I → d J X, with I = 2, 3, J = 1, 2 the quark generation indices. If its Wilson coefficient is set to one, then there is a scale Λ such that the NP contribution equates the SM prediction for d I → d J νν, with m 2,3 = m K,B , g the SU (2) L coupling constant, and V the CKM matrix. As shown in Table 1 as a function of the dimension, the SM loop factor combined with the CKM suppression pushes the scales Λ well above the electroweak scale for n 7. On the contrary, the rare decay constraints cease to make sense for SU (2) L ⊗ U (1) Y invariant operators of dimension n 9, since powers of (H † H)/Λ 2 → v 2 /Λ 2 grow unchecked when Λ v, where v ≈ 246 GeV is the SM Higgs vacuum expectation value.
Clearly, these scales are only indicative. The true sensitivity to a given dark operator depends essentially on two additional factors. First, the quark transitions d I → d J X have to be probed through hadronic processes. Hence, depending on the modes, hadronic matrix elements as well as phase-space factors can alter significantly the estimates of Table 1. In the following, we compare the sensitivities of all the leading modes. Specifically, in the K sector, we include the modes with the least number of pions and photons in the final states, i.e. K → X, K → πX, K → γX, and K → ππX, and leave out the K → πππX modes. Similarly, in the B sector, the considered modes are the B → X, B → (K, K * )X, and B → (π, ρ)X decay channels. The γX channel, driven in the K sector by the QED anomaly, is suppressed and difficult to reconstruct experimentally in the B sector, and will thus not be included [6].
The second factor determining the true sensitivity of a given mode is related to the experimental strategies deployed to measure it. Since invisible states are not seen, the kinematical reconstruction is limited. In addition, these modes are so rare (in the SM) that they require very aggressive background suppressions. To this end, the central tool is the differential rate in terms of the kinematical parameters of the visible products. But this differential rate depends on the nature of the dark particle. Currently, most experimental analyses implicitly impose the SM differential rate (for X = νν). This means that the current bounds cannot be directly translated to other types of final state particles. This motivates another goal of the present paper, which is to provide the full dictionary of the differential rates for all the leading effective interactions involving invisible final states.
These spectra should be used by the experimentalists to derive bounds for each type of new invisible state.
For more details on these issues, including kinematics, matrix elements, current measurements or bounds, and experimental prospects for the various modes, we refer to appendix A.1 (B.1) for K (B) decays.

Flavor-based classification of the dark operators
At the electroweak scale, once the whole NP particle spectrum but the X has been integrated out, the lowest dimensional operators can be split into three types according to their quark and lepton field contents: By definition, H int contains only gauge and Higgs fields, while H mat and H ∆B,∆L contain at least one SM fermionic field. The operators of H ∆B,∆L have a non-zero charge under the baryon (B) or lepton (L) number U (1)s. As in Ref. [1,2], all the operators are to be written as manifestly invariant under SU (3) C ⊗ SU (2) L ⊗ U (1) Y , i.e. in terms of the quark (lepton) doublets Q (L) and singlets U, D (E) of each flavor, of the SU (3) C , SU (2) L , U (1) Y field strengths G a µν , W i µν , B µν , of the Higgs doublet H, as well as of covariant derivatives acting on these fields, insofar as these cannot be reduced using the SM equations of motion (EOM).
Due to their different field contents, these three types of operators do not contribute equally to the quark transitions s → dX, b → dX, and b → sX, so let us organize them differently, in terms of four classes of scenarios, as shown in Table 2.
Consider first the operators of H mat involving down-type quarks (those with leptons or up-type quarks are obtained by substituting D, Q → E, L or D, H → U, H * ). Up to possible partial derivatives acting on the quark or invisible fields, and omitting the Dirac structures, the quark currents are Those operators have a flavor structure, and thus can in principle induce d I → d J X. Clearly, when analyzing the physics reach of rare K and B decays in terms of the scale Λ, the assumptions made on the c I =J are crucial. There are two main scenarios: I. The constraints derived from rare FCNC decays are the tightest when the NP flavor structure is generic, As shown in Table 1, the bounds on the NP scale Λ are then often far above the electroweak scale.
II. Since H mat results from integrating out the whole NP particle spectrum, the flavor-breaking character of its operators could originate from dynamical effects not related to the dark sector. In that case, the NP dynamics would also quite naturally correct the visible FCNC operators, on which there are many tight experimental constraints from K and B physics [8]. This is typically the case in supersymmetric settings, where the flavored soft-breaking terms cannot be fully generic. Phenomenologically, a simple way to account for such a non-generic NP flavor structure is to impose the Minimal Flavor Violation (MFV) ansatz [9], i.e. force the quark currents to have the forms In the down-quark mass-eigenstate basis, the diagonal vY d = √ 2m d tunes the chirality flips, while vY u = √ 2m u V parametrizes the flavor change (m u,d denotes the diagonal quark mass matrices, V the CKM matrix, and v ≈ 246 GeV the Higgs vacuum expectation value). So, MFV rescales the Wilson coefficients according to c IJ Upon these rescalings, the accessible scales Λ are then much lower, especially for operators of low dimensions, and for s → d operators involving light right-handed quarks.
If the whole NP dynamics is flavor blind, then X couples only to the flavor-diagonal quark currents and c I =J = 0. A flavor transition is of course still possible but it must proceed through the SM weak interactions, i.e. the flavor-blind quark currents must be dressed by a flavor-changing W interaction, either between the quark lines or as a self-energy on these quark lines.
In particular, all the operators of H int are of this type. Indeed, to contribute to FCNC processes, gauge fields have to be coupled to quarks, while Higgs fields are either coupled to quarks or left as external tadpoles (to be replaced by the vacuum expectation value after the electroweak symmetry breaking). The resulting couplings between quarks and X can then be matched onto H mat , and satisfy 1 c I =J ≈ 0. Thus, the only difference with respect to the flavor-blind H mat operators is that the c II coefficients are initially suppressed by some power of the SM coupling constants (at the scale Λ), by some loop factors, and by quark Yukawa couplings when a Higgs field is coupled to the quark line (so that c 33 ≫ c 11,22 ).
For reasons entirely pertaining to the SM dynamics, it is different to probe the couplings to heavy quarks, c 33 , and to light quarks, c 11,22 , so these constitute our third and fourth classes of scenarios.
III. Let us assume we have an operator coupling X to the top quark current. Dressing it with a W exchange (see Table 2), the necessary GIM breaking arises at the electroweak scale. From the rare decay perspective, this electroweak physics is local, so it can be matched onto the flavor-changing operators of Class II. The Wilson coefficients end up suppressed by the MFV scalings (6) and by a loop and gauge coupling, i.e.
Typically, the bounds on the scale Λ are brought down very significantly, often at around the electroweak scale. Still, the rare K and B decays remain ideal probes for such kind of effective couplings since a direct collider signal in the tt channel is presumably hidden by the large flavor-blind SM backgrounds.
IV. With X coupled to the light quarks (u, d, or s), it is much trickier to derive bounds on the scale Λ for three reasons. Firstly, B physics is unable to provide competitive constraints, given the small CKM factors V * ub V us or V * ub V ud . Secondly, though in the K sector the weak transition is favored by the large V * us V ud for CP-conserving observables 2 , the light quarks are never integrated out and remain dynamical. This renders the theoretical control difficult since the K physics scale is too low to allow for a perturbative QCD treatment. This is true both for the effective operators describing the weak transition and for the effective operators describing the production of the new invisible states. Thirdly, there are already many constraints on the flavor blind production of invisible particles, in particular from π 0 or quarkonium decays. Compared to the other classes, it is a priori not clear whether rare decays offer privileged windows. In the following, we will consider only specific scenarios where competitive bounds can be derived.
The final type of operators is a bit different and does not immediately fit in the above classification. First, in general, ∆B and ∆L dark operators cannot be present simultaneously, since an X exchange would induce proton decay. A way out would be to impose MFV, which forces ∆L operators to be proportional to the tiny neutrino masses [11], rendering them irrelevant. The ∆B operators are not particularly suppressed under MFV, but they never contribute to the FCNC-induced rare decays considered here. Indeed, these modes trivially conserve B since an odd number of baryons would be required in the final state. This is not possible for K decays, while the signatures for B decays should be experimentally clear, but are beyond our scope.
If MFV is not imposed and ∆B operators are somehow disposed of, then ∆L contributions to the d I → d J X processes are possible. But, the only ∆L = 0 invisible final states either involve an odd number of neutrinos, 1 The contributions to c I =J are tiny because, as long as the SU (2) L ⊗ U (1) Y symmetry is exact, the gauge/Higgs fields of the H int operators can couple to a flavor-changing quark current only at the cost of at least two Yukawa insertions, i.e. two Higgs tadpoles, on the quark line. At the high scale, these tadpoles cost a factor Λ −2 compared to the flavor-blind version of the same operator. When Λ is large, flavor-changing effects thus dominantly originate from electroweak scale GIM breaking, even though this necessitates an extra W loop (see Table 2).
2 CP-violating observables, for which Im(V * us V ud ) = 0 (this standard CKM convention is used throughout the paper, see e.g. Ref. [10]), are induced by heavy quarks and fall into Class III. Due to the scaling (7), they are very suppressed.
Heavy quark: q = (c), t Light quarks: q = u, d, s, (c) Bounds on c IJ /Λ n directly derived from the d I → d J X processes. When MFV holds, c IJ ∼ V * tI V tJ times the appropriate chirality flip factors m d I,J /v, see Eq. (6).
Same local operator basis, but with the coefficients rescaled as c IJ → (g/(4π)) 2 V * tI V tJ × c 33 times the appropriate chirality flip factors m d I,J /v, see Eq. (7).
Due to the small V * ub , B decays are not competitive. For K decays, the q = u contributions are dominant but non local, and require controlling long-distance hadronic effects. Class I, II Class III Class IV or some ∆L = 2 neutrino pairs ν L ν L . Since a neutrino field in an effective operator costs Λ −3/2 , these are in general significantly suppressed compared to the operators of H mat and H int . The only exceptions are those contributing to d I → d J ν L ψ or d I → d J ν L Ψ. As will be discussed in the relevant sections, because ν L is part of the lepton doublet, these operators are always accompanied by the charge-current transitions d I → u J ℓ − ψ or d I → u J ℓ − Ψ, which may offer better windows than the rare FCNC transitions.
2 Invisible spin-1/2 fermion When the new invisible fermion is neutral under the SM gauge group and is produced in pairs, imposing the SU (3) C ⊗ SU (2) L ⊗ U (1) Y gauge invariance reduces the basis to the usual ten chiral currents: with the definition σ µν ≡ i(γ µ γ ν + g µν ), and where Q (D) stands for the left-handed quark doublet (righthanded down-quark singlet). Similar operators can be written down for the up-quark right-handed singlet or for leptonic transitions, and the generalization to a two Higgs-doublet model is straightforward. The coefficients are not assumed real, and their flavor indices are understood. For example, written in full for the s → d, b → s, and b → d sectors which concern us here:  are indicative; the bounds are derived for m X = 0 or, if the rate vanishes (signaled by (*)), for m X = 100 MeV. For those which do not vanish, most V − A, S − P , or T −T degeneracies are lifted when m X = 0. Still, the dependences of the Λs on m X are rather weak. The bounds stay within an order of magnitude of their values at m X = 0, except for m X close to the kinematical threshold, where they sharply drop towards zero. This is shown in the last column for K + → π + νν, where Λ (normalized to the value quoted in the table) is plotted against m X over the range [0, m π ].
There are only two tensor operators because the identity 2σ µν γ 5 = iε µναβ σ αβ forbidsQσ µν D ×ψ R σ µν ψ L and Dσ µν Q ×ψ L σ µν ψ R . The dimension-seven operators involving the SM Higgs field are retained because their suppression is only by v/Λ. By contrast, true dimension-seven operators involving an additional derivative are more severely suppressed by m K,B,ψ /Λ, and are thus not included.
After the electroweak symmetry breaking, the gauge-invariant operators are rewritten in terms of the (pseudo)scalar, (axial)vector, and (pseudo)tensor currents. This minimizes the interferences between the currents in physical observables, since these operators dominantly produce the invisible fermions in different states. Table 4: Pair production of invisible states: Scales Λ (in TeV) accessible for the various operators with present (future) measurements (see App. B.1). The operator dimensions d are written with a bar for those involving a Higgs field, and thus a v/Λ factor. We assume m X = 0 everywhere, except when the rate vanishes in this limit (indicated by *). Specifically, the B → XX bound on f AA , the B → π(K)XX bounds on g T , f T S,T P , and the B → K * XX bounds on g T,V,A are derived for m X ≃ 2 GeV. The behaviors of the scales Λ as functions of m X are shown for X = ψ in Fig. 2.
The change of basis is where f XY tunesq I Γ X q J ×ψΓ Y ψ with Γ V,A,S,P,T,T = γ µ , γ µ γ 5 , 1, γ 5 , σ µν , σ µν γ 5 , and flavor indices are understood. The corresponding differential rates are listed in Appendix A.3 (B.3) for K (B) decays, and the physics reach is summarized in Tables 3 (4). Specifically, the entries of these tables are obtained by turning on each f i in turn, while keeping the others to zero. We do not do this at the level of the c i in order to minimize the interferences among the NP contributions. The value of each |f i | is set to one for the vector and axial vector currents, and to v/Λ for the others, in order to keep track of the SU (2) L ⊗ U (1) Y structure of the underlying operators. Since K L,S are approximate CP eigenstates, the CP phase of the f i must be kept arbitrary. Effectively, we turn on Im f i and Re f i to one (or v/Λ) separately, and the tightest bound (largest Λ) is indicated in Table 3.
To derive the entries of Tables 3 and 4, specific experimental bounds on the rare decay branching ratios are used, as detailed in App. A.1 (B.1) for K (B) decays. Note that the scales Λ corresponding to tighter or looser experimental bounds are immediately obtained by a simple rescaling of the numbers in Tables 3 and 4, given the dimensions of the operators (8) and the definitions (10).  Table 4, except for B → ψψ, given at m ψ = 2 GeV. Note that these B → (K ( * ) , ρ, π)ψψ bounds assume flat experimental acceptances and full phase-space coverage, which is not true in the present experimental analyses, see App. B.1.
Finally, in Fig. 3, we compare the sensitivity of the various K and B decay modes for two illustrative examples, c V LL and c S LL , with and without imposing MFV. In the former case, the coefficients are set to (c V LL ) IJ = λ IJ and (c S LL ) IJ = λ IJ m d I /v, see Eq. (6). Actually, to avoid dragging a relative factor m b /v, the scales Λ are plotted taking a two Higgs doublet model of type II at large tan β = v u /v d (where v u and v d are the two Higgs vacuum expectation values), so that m b /v d ≈ 1 and (c S LL ) sd is simply suppressed by m s /m b . As can be seen in Fig. 3, this chirality flip is expensive in the K sector, and pulls the K → πψψ constraint an order of magnitude below those from rare B decays. By contrast, for the vector current, similar constraints are drawn from rare K and B decays when MFV is imposed.
Let us now turn to the operators not involving the SM fermions, but rather the SM gauge and Higgs fields. Using these fields EOM, as well as partial integration and identities like 2σ µν γ 5 = iε µναβ σ αβ , the leading operators reduce to where ← → D µ = ← − D µ − − → D µ and Q G i stand for each of the following gauge-invariant combinations of SM fields But for the last three CP-odd combinations, the Q G i monomials are precisely the dimension-four gauge and Higgs couplings of the SM Lagrangian. Operators with partial derivatives acting on ψ are systematically discarded. The leading dimension-five operators involve either B µν or H † H. The B µν operators effectively assign a (mass-dependent) hypercharge to the invisible fermions, since under partial integration With m ψ /Λ ≪ 1 but otherwise arbitrary, the second term describes millicharged fermions, a scenario to be discussed later on. The H † H operators effectively correct the ψ mass after the electroweak symmetry breaking, with (setting c H LR,RL to one) Taken at face value, naturality would thus prefer Λ > 24 (240) TeV if m ψ ≈ δm ψ < m B /2 (m K /2), as required to produce these states in rare B (K) decays. As seen in Tables 3 and 4, this would push the NP effects from the FCNC operators (8) just beyond accessibility. Of course, these numbers are only indicative, a specific model need not produce these dimension-five operators, and even if generated, a sizeable deviation from naturality cannot be ruled out. All the other Hψ ψ int operators are, for our purpose, only marginally relevant. As explained in Sec. 1, they reduce to the four-fermion operators of Eq. (8) once Higgs and gauge fields are coupled to quarks, with the four-fermion Wilson coefficients obeying the SM scalings (7). For example, the dimension-six operators with H † ← → D µ H reduce, after the electroweak symmetry breaking, to an effective Z coupling to ψ which can be treated in perfect analogy to the SM Z coupling to νν (see Fig. 1). The rare decays proceed through the flavor-changing hadronic Z penguin, and the four-fermion operators scale like c IJ LL(RR) ∼ gc H LL(RR) k IJ with k IJ given in Eq. (7). The final class of operators involves only one dark fermion field, and thus, from Lorentz invariance, violates either lepton (L) or baryon (B) number: For simplicity, the possible Dirac, SU (2) L , and flavor structures are understood, as well as the operators with ψ R → ψ C L . Operators with three dark fermion fields are at least of dimension seven, hence not included. As explained in Sec. 1, the phenomenology of the ∆B or ∆L operators is different from that induced by those of Hψ ψ mat and Hψ ψ int . Since c ∆B and c ∆L i cannot be simultaneously large as a tree-level ψ exchange would induce proton decay, and since c ∆B cannot be probed with the rare FCNC-induced decays considered here, let us concentrate on the ∆L operators. The renormalizable c ∆L 0 interaction, as well as c ∆L 1 after the electroweak symmetry breaking, mix ψ with ν L . As a result, these operators are bounded by neutrino masses, ψ R behaves effectively as a right-handed neutrino, and the c ∆L 2,3 effective interactions include both FCNC and chargedcurrent interactions. Since the renormalizable operator (the so-called neutrino portal) is not suppressed by the NP scale, its coefficient c ∆L 0 must be tiny. So, the effective interactions tuned by c ∆L 2,3 can be sizeable only if the NP dynamics responsible for such a suppression does not apply equally to all the ∆L operators.
In this respect, the semileptonic operator tuned by c ∆L 3 may be the most likely to escape such a suppression. It is the only one directly accessible with the quark FCNC transitions. However, the processes induced by the charged current component of the c ∆L 3 operator,d I (1−γ 5 )u J ×ψ(1−γ 5 )ℓ K , may give tighter constraints. Indeed, since there is no possible interference with the SM contributions, this operator enhances all the kinematicallyallowed (semi-)leptonic decays of the K + , B + , and D + , as well as hadronic τ decays. In particular, the most interesting channels are the K + , B + , D + → ℓ + ψ transitions, whose rates are where P = B, K, D for I, J = (3, 1), (2, 1), (1, 2), r i = m i /m P , f P is the P ± decay constant (defined as in Ref. [12]), and the kinematical function λ is defined in Eq. (95). Indeed, the SM charged-current Fermi operators are vectorial, hence these transitions are helicity-suppressed in the SM (proportional to m ℓ K ). On the contrary, the relative strengths of the P + → ℓ + ψ transitions for ℓ = e, µ, τ depend only on that between (c ∆L 3 ) IJK , K = 1, 2, 3. For example, if (c ∆L 3 ) IJK is independent of K, the NP effects would be most easily seen for eψ final states.
Such an enhancement may be welcome in the B sector. The persistent tension between the Belle and BaBar measurements of B(B → τ ν) = (1.68 ± 0.31) · 10 −4 [13][14][15][16] and the predictions of the global CKM fits within the SM (see Ref. [17] for a recent review) can be addressed provided m ψ < 3.5 GeV. Assuming (c ∆L 3 ) 313 ≈ 1, the reconciliation of the discrepancy of the order ∆B(B → τ X) ≃ 10 −4 would point towards Λ 5 TeV, where the equality is reached for m ψ = 0. Note, though, that a non-universal lepton flavor structure in (c ∆L 3 ) IJK is necessary, in order to circumvent the current bounds from B → eν and B → µν at the level of 10 −6 [12].
Actually, if the c ∆L 3 couplings are fully universal in their quark and lepton indices (i.e. (c ∆L 3 ) IJK ∼ O(1) for all I, J, K), and if m ψ ≪ m K , the best constraints currently come from the K ℓ2 universality test 3 , where the NP effect gets enhanced by the small electron mass: From R SM K = 2.477(1) · 10 −5 [18] and R exp K = 2.487(13) · 10 −5 [19], we require R exp K − R SM K 0.013 · 10 −5 which translates as Λ 82 TeV for m ψ ≪ m K . This is comparable to the scales probed with the FCNC modes.

Some models with light weakly coupled fermions
The operator basis (8) can be used to describe the SM transitions to neutrino pairs, ψ L = ν e , ν µ , ν τ , ψ R = 0, by setting all the c i to zero but for c V LL . The SM rates for K decays are given in App. A.3.1, and those for B decays in App. B.3.1. With ψ R = 0, all but the vector FCNC operators c V LL and c V RL drop out. Since most NP models do not modify the particle spectrum below the electroweak scale, their effect simply enter as corrections to these two coefficients [21,22], and lead to the same vector-like spectrum as the SM contribution. Note that contrary to true dark fermions, the neutrinos are not neutral under SU (2) L ⊗ U (1) Y , and thus the basis (8) is not complete since it is lacking the charged current semileptonic operators. Said differently, if neutrinos had not yet been discovered, charge-current interactions would offer better windows.
The presence of a light right-handed (sterile) neutrino corresponds to ψ R = ν R , ψ L = 0 in Eqs. (8) and (15). Of course, once ψ R is identified with ν R , it can be given a charge under U (1) L so that all these operators become ∆L = 0. Some NP dynamics is nevertheless needed to couple ν R , a gauge-singlet, to the quark currents. This can arise in the νMSM model [23], based on the c ∆L 0 coupling and with the quark flavor transition induced by the SM weak interaction. Alternatively, in some left-right symmetric models [24], the c ∆L 4 operator can arise from combined W L and W R boxes. Note that these two examples are matched onto the operators of Eq. (15), for which lepton universality tests are competitive with FCNC decays, as our analysis of the previous section shows.
In the MSSM, the lightest neutralino χ 1 is another electrically neutral weakly interacting particle which could be pair produced in rare decays. Though this particle is not neutral under the SM gauge group, it must be produced in pairs when R parity is conserved. Because the neutralino is a Majorana fermion obeyinḡ the corresponding operators disappear. In the general MSSM, there are then tree-level processes contributing to the combinations c V LL + c V LR , c V RL + c V RR , and to the scalar and pseudoscalar currents c S XY , X, Y = L, R. The flavor transitions are tuned by the off-diagonal entries of the down squark LL, RR, and LR mass terms, and the overall scale Λ is set by the exchanged down squark mass. This (Class II) scenario was analyzed in detail in Ref. [25], to which we refer for more information.
Another possible new fermion is the axinoã, the fermionic superpartner of the axion [26]. Depending on the model, it could be light enough to be produced in K and B decays. Note, though, that its Lagrangian couplings are flavor-blind, and further, involve the superpartners of the SM particles when R parity is conserved. So, the effective interactions are not only suppressed by the large scale Λ = f a , but also by the sparticle mass scale, loop factors, and the already tightly constrained flavor-violation occurring in the squark sector. We will not consider this scenario further because, with scales f a above 10 6 TeV [26] (or even much higher in some models), the fermionic operators are far too suppressed, and signals should be more readily accessible using the single scalar axion production discussed in the next section.
A final example is the dark sector millicharged fermion ψ ε [27]. Typically, these fermions arise when there is a new dark U (1) field V µ coupled to the SM U (1) Y through the kinetic mixing εB µν V µν , as well as to some new fermion states initially neutral under the SM gauge group [28]. After diagonalizing the two U (1) fields, the dark sector fermions end up coupled to the SM photon, but with an arbitrary electric charge εe. Alternatively, this scenario could follow from the B µν operator in Eq. (13).
When m ψε m e , various very tight cosmological and astrophysical bounds hold on ε [27], ruling out any signal in meson decays. On the contrary, for m ψε ∼ O(100 − 1000 MeV), the bounds are much less tight, with ε as large as 10 −2 still possible. Since ψ ε couples to quarks through the photon field, its coupling is flavor-blind (Class III or IV of Table 2). The bounds in Table 3 are adapted to this scenario by setting Λ = M W and rescaling the Wilson coefficients as c IJ → εe 2 k IJ with k IJ given in Eq. (7). Alternatively, the physics reach can be more accurately estimated by looking at the K and B radiative modes with a Dalitz pair, and rescaling their branching ratio by ε 2 , Up to simple phase-space corrections accounting for m ℓ = m ψε , this is valid as long as the Z boson does not dominate. From this, no B decay rate appears large enough to reach the interesting ε 10 −3 region.
For example, with B(B → (K, K * )ℓ + ℓ − ) in the 10 −7 range, a bound on B(B → (K, K * )ψ ε ψ ε ) at the 10 −11 level would be needed, far below the experimental prospects. Similarly, in the K sector, assuming about 10 13 kaon decays, the only mode for which one could theoretically reach down to ε ≈ 10 −3 is K L → γψ ε ψ ε , since B(K L → γe + e − ) = (9.4 ± 0.4) · 10 −6 [12]. So, the rare FCNC decays do not appear competitive when ψ ε couples to ordinary matter exclusively through the photon. But, let us stress that if the couplings to matter of the SM and dark U (1)s are not perfectly aligned, the CKM suppression of Eq. (7) may be evaded, and rare decays would become prime sources of information. Also, as we will see in Sec. 4.3, in case the dark photon has a non-zero mass, it may be directly produced and competitive constraints can be derived.

Invisible spin-0 boson
If there is a light scalar particle neutral under the SM gauge group and under any dark gauge symmetry, it can be produced alone. The simplest effective operators are then of dimension five No operator involving the tensor current or the SM field strengths arises, as these would require more derivatives and thus would be suppressed by additional factors of O(m K,B /Λ). Actually, the first two operators (as well as (Q DQ)φ, not explicitly included) collapse to the third and fourth upon using the tree-level quark EOM [1] i i.e., at the cost of the chirality flips Though in the rest of the paper, the EOM are always enforced, we prefer to keep all three operators here because they correspond to well-defined scenarios. On one hand, derivative couplings are characteristic of nonlinearly realized symmetries, while on the other hand, the H †D Q × φ operator is effectively a dimension-four Yukawa-like coupling after the electroweak symmetry breaking. If the scalar field is charged under some dark sector symmetry, or if the Z 2 symmetry under φ → −φ is imposed, it must be produced in pairs. The simplest effective operators are then of dimension six: The minus sign of ← → ∂ in the first two operators is required as the plus combinations reduce to the third and fourth operators by partial integration and use of the quark EOM (21). Tensor currents are not included since they are suppressed by a factor of O(m K,B /Λ).
Bounds are derived on the operators involving quark currents of definite C and P , tuned by the couplings where flavor indices are understood 4 . The corresponding differential rates are listed in Appendix A.4 (B.4) for K (B) decays. The physics reach is summarized in Tables 5 and 6 for one scalar in the final state, and in Tables 3 and 4 for two scalars. As expected from Table 1, the lower dimensionality of the operators (20) and (23) translates as much higher accessible scales compared to spin 1/2 final states. 4 Note that the c φ(φ) V R,V L couplings are hermitian matrices in flavor space (see Eq. (9)), while (c Table 5: Production of a single invisible particle: Scales Λ (in TeV) accessible for the various operators, assuming bounds on the branching ratios of 10 −10 for each mode. As for Table 3, the ranges of accessible invisible masses indicated in the first line are indicative, as the bounds are derived setting m X = 0 (except for those channels which vanish, denoted by (*), for which m X = 100 MeV). These scales naively decrease when m X increases due to the phase-space suppressions, though the experimental acceptances need to be taken into account (see App. A.1). For the production of an invisible vector, see also Fig. 4. Table 6: Production of a single invisible particle: Scales Λ (in TeV) accessible for the various operators with present (future) measurements (see App. B.1). We assume m X = 0 everywhere, except for those channels which vanish, denoted by (*), for which m X = 2 GeV. For the production of vector states, see also The simplest effective operators involving Higgs and gauge fields are where Q G 1−8 are defined in Eq. (12). Only operators of dimension less or equal to those in Eqs. (20) and (23) are considered. The possibility to write renormalizable couplings between φ and H embodies the so-called Higgs portal. One consequence is a mixing between φ and H, which has already been investigated in details, see e.g. Ref. [5], and will not be considered further here (but for a comment on λ ′ in the next section). The other operators are a priori subleading compared to the Higgs portal, simply because of their higher dimensions. Note, though, that when φ is pseudoscalar, most of the H φ,3φ int couplings drop out in the CP limit, leaving only those to F µνF µν with F µν any one of the SM field strengths. Finally, the presence of a neutral scalar field does not open new possibilities compared to the SM to construct ∆B and ∆L couplings. The simplest operators are obtained by multiplying by φ (or by φ † φ if φ is not neutral) those of Ref. [2], and are either of dimension seven and ∆B = ∆L = 1, or dimension six and ∆B = 0, ∆L = 2: The ∆L = 2 operator can produce an invisible ν L ν L φ final state. However, not only is this operator of higher dimension compared to those of Eq. (20), but its contribution to the rare decays proceeds through a neutral Higgs penguin (Class III in Table 2), and is thus extremely suppressed. The ∆B = ∆L = 1 operators have no impact on the ∆B = 0 rare FCNC decays considered here. In addition, if m φ < m p,n , those involving light flavors can induce ∆B = 1 proton or neutron decays. In that case, without highly non-generic flavor structures for the c ∆B=∆L i , the scale Λ must be close to the Planck scale. On the other hand, if m φ > m p,n , these operators could still induce exotic B decays into an odd number of baryons plus missing energy.

Some models with light weakly coupled (pseudo)scalars
Typical non-linear symmetry realizations lead to derivative couplings of the type c φ V L,V R to the dark scalar field, with Λ given by the symmetry breaking scale F . Well-known examples of such light scalar states are the axion, resulting from the breaking of the PQ symmetry [3], or the familon, originating from the breaking of some family symmetry [29]. Only the latter naturally lead to FCNC couplings, since by design the family symmetry relates the three generations. For the axion, the dominant effect comes from its flavor-blind coupling to light quarks (Class IV in Table 2). This is dominated by long-distance effects, specifically by the mixing of the axion with the light neutral mesons, as analyzed e.g. in Ref. [30] to which we refer for more details. Let us mention also that in some axion models, there is no direct coupling to light quarks, but rather couplings to the QED or QCD field strengths of the form φF µνF µν or φG a µνG µν α , as included in Eq. (25). Similar derivative couplings arise from models of meta-stable supersymmetry breaking [31], where a light pseudo Nambu-Goldstone boson φ = P may be present. As discussed e.g. in Ref. [32], through the exchange of three gauge bosons, axion-like effective couplings to quarks are generated with c φ V L,V R ∼ α 3 i (Λ) and α i either the weak, strong, or hypercharge coupling. The scale Λ is the supersymmetry breaking scale, which can be around 10 − 100 TeV in gauge-mediated supersymmetry breaking scenarios. These couplings are dominantly flavor blind, so according to our general rule (7), the flavor-changing vertices scale 5 as c φ,IJ V L,V R ∼ α 3 i (Λ)k IJ . Thus, rare K and B decays seem unable to reach the interesting range Λ > 10 TeV. If supersymmetric particles are allowed to propagate in the loop(s), the operator may be enhanced, though the tight constraints on the flavor violation in the squark sector probably limit the accessible scales to Λ 10 TeV range.
Another scenario involving dark light scalar particles is the Next to Minimal Supersymmetric Standard Model (NMSSM). But, in this case, the scalar mass is often larger than m B [33], the scalar may decay too fast (e.g. to ℓ + ℓ − ) to be considered as an asymptotic state in K and B decays, and the flavor transitions need to be induced by the (supersymmetrized) weak interaction. When MFV holds, the effective operators are then very suppressed (Class II scenario, see Eq. (6)). For all these reasons, the NMSSM does not appear as a likely scenario where our effective operators could play a role. Still, our formulas for the differential rates can be directly adapted to that case, and improve on the analysis of Ref. [34] by including more observables, and by a better treatment of the hadronic matrix elements.
A more generic example is the singlet scalar model for dark matter, denoted S. In its simplest form, it enforces the Z 2 symmetry and includes only the renormalizable coupling λ ′ of Eq. (26), with φ = S a real field. The capabilities of rare decays in probing this scenario were discussed in Ref. [35]. As described there, being flavor blind, the weak interactions are needed to induce the c φφ SR,SL coupling from the λ ′ one (which couples SS to tt through a tree-level Higgs exchange). So, in this case, the key to interpret the numbers in Table 3 is the identification c φφ,IJ SR,SL /Λ 2 → λ ′ k IJ /m 2 h , with k IJ given in Eq. (7). See Ref. [35] for more details. As a final example, the new invisible scalars could be the sgoldstinos S and P , the scalar superpartners of the goldstino [36], for which the scale Λ corresponds to the fundamental scale of supersymmetry breaking. If light enough, these scalars could be produced alone or in pairs in K and/or B decays. In the former case, the coupling is of the form of c φ SR,SL , and is tuned by the LR squark mass insertions. In the latter case, the derivative interactions are of the form of c φφ V L and c φφ V R with φ = S and φ † = P , and are tuned by LL and/or RR squark mass insertions. Note that since these squark mass insertions are rather tightly constrained by visible K and B observables, the accessible scales are limited by the MFV rescalings (6).

Invisible spin-1 boson
Compared to scalar and fermionic invisible particles, the presence of an invisible vector particle (denoted V ) is significantly more difficult to parametrize. A consistent description of a neutral massive vector boson coupled to SM fermions is notoriously delicate. At the same time, adding a U (1) gauge group to the SM is one of the simplest possible extensions [28,37,38]. Furthermore, various theoretical models, whose motivations originate for example from strings [39], extra dimensions [40], or dark matter theories [41], predict such new long or medium range forces, i.e. with masses in the MeV to a few GeV range. Such a dark vector boson would have many phenomenological implications, and has been intensely investigated recently [42]. It could also be produced in rare B and K decays, where it would show up as missing energy if sufficiently long-lived. In this respect, most models do also induce a coupling to leptons, but as long as the V → ℓ + ℓ − vertex is not significantly larger thanqq → V , our bounds should hold since producing a Dalitz pair through a virtual V exchange would push the rates well beyond the experimental reach.
There are several ways to deal with light vector states, which we organize into three scenarios. For the first, the simplest FCNC operators involving the dark vector field are constructed. Those lead to decay rates diverging in the m V → 0 limit. This singularity is then treated assuming some kind of Higgs mechanism takes place in the dark sector (in a way similar to Ref. [38]). For the second scenario, the vector field is supposed to couple to SM matter fields only through its field strength, in a (dark) gauge-invariant way. This automatically ensures a safe m V → 0 limit, but significantly increases the dimensionality of the effective FCNC couplings. Finally, the third scenario considers only low-energy effective couplings of V to conserved flavor-blind quark currents [37], as can arise for example from the couplings of V to SM gauge fields through the kinetic mixing [28].

Simplest FCNC operators
For the first scenario, we consider the lowest-dimensional flavor-changing operators involving a vector field neutral under the SM gauge group, which are simply Thanks to the Lorentz condition ∂ µ V µ = 0, the leading correction starts at O(1/Λ 2 ). The flavor-changing quark currentsQ I γ µ Q J andD I γ µ D J are not conserved when I = J (see Eq. (21)). As a result, a naive rate computation with the polarization sum diverges as m V → 0. This well known phenomenon is related to the impossibility of defining consistently the massless limit without an active gauge symmetry. It is interesting to compare this to the behavior of the dimension-four operators originating from the Z penguin in the SM (see Fig. 1). There, the Ward identity is violated only once the electroweak symmetry breaking takes place, and With Λ SM ∼ v, this effective interaction is simply proportional to g 3 , up to loop factors. However, the mass of the Z is also of O(v), and thus it can never be on-shell once using the effective operator formalism. Instead, dimension-six four-fermion operators proportional to g 4 /M 2 Z ∼ g 2 /v 2 are relevant at low energy, see Eq. (1). In our case, V is light enough to be produced in meson decays. If it also gets its mass through some spontaneous symmetry breaking, then 6 m V ∼ ε V L,R v dark for some v dark presumably similar or larger than v.
But as long as v dark > 0, the limit m V → 0 requires ε V L,R → 0, which never diverges for physical processes. Actually, enforcing ε V L,R ∼ m V /v dark in the m V → 0 limit, any decay rate behaves as with M µ the hadronic matrix element P ′ |Q I γ µ Q J |P or P ′ |D I γ µ D J |P and dΦ P ′ V the phase-space integrations. As expected from the equivalence theorem, this is precisely the rate one would derive from the axionic For the pair-production of dark vectors, the simplest operators are of dimension six: which are orthogonal to those tuned by c V V DL,DR (which are missing in Ref. [6]). These operators are given for completeness, but will not be considered further for two reasons. First, in the present minimal theoretical setting, there is no reason for the renormalizable operators of Eq. (28) to be absent, and those would clearly offer better windows. Second, the leading operators of H VV mat [I] are those tuned by c V V SL,SR since the others are comparatively suppressed by m K,B /v. But these operators also arise from the H †D Q × V µν V µν and HQD × V µν V µν operators considered in the next scenario (albeit in a rescaled form, see Eq. (38)). As explained there, these operators could become leading in the presence of a non-abelian dark gauge invariance, which would forbid single vector production.

Phenomenology
Let us consider separately the vector When m V = 0, both terms of the polarization sum (29) contribute, so experimental bounds translate as constraints in the (m V , ε V ) or (m V , ε A ) planes. As explained above, we identify the regions where m V /ε V,A > v as physical (for a similar reasoning, see e.g. Ref. [38]).
As shown in Figs. 4 and 5, the pair (m V , ε V ) is very constrained in both the K and B sectors, with typically m V /ε V ≫ v. Indeed, because the second term of Eq. (29) is growing as 1/m 2 V in the m V → 0 limit, ε V has to be correspondingly tiny to pass the bounds on the branching ratio. Said differently, the dark spontaneous symmetry breaking scale v dark in Eq. (31) has to be much larger than the electroweak scale. Interestingly, this remains true even when the flavor-violating part of ε V satisfies MFV (green regions in Figs. 4 and 5), i.e. is rescaled as ε V → |V * tI V tJ |ε V , see Eq. (6). So, even if there is no theoretical requirement for the hidden symmetry breaking scale to be very different from that of the visible sector, current FCNC constraints nonetheless require v dark ≫ v.
For comparison, the right panel in Fig. 4 shows the constraints one would get from a similar bound at the 10 −10 level on the K → γV channel. As this process is induced by the anomaly, the coupling is of the form with F µν the QED field strength. So, the current is effectively conserved, the 1/m 2 V term in Eq. (29) drops out, and the m V → 0 limit becomes smooth.

Gauge invariant FCNC operators
For the second scenario, we impose that only the dark field strength V µν ≡ ∂ µ V ν − ∂ ν V µ and its dualṼ µν ≡ 1 2 ε µνρσ V ρσ occur in the flavor-changing couplings. This restores gauge invariance and ensures a smooth m V → 0 limit, but increases the dimension of the simplest operators. Depending on whether the vector is assumed abelian or non-abelian, the lowest-dimensional operators are either dimension six, or dimension eight, To reach this minimal basis, a number of identities and approximations were used. First, for the operators with a single vector field, the dual field strengthṼ µν is absent from H V mat [II] since it can always be reduced to operators involving V µν using the quark EOM, integration by part, and the Chisholm identity [1].
For the operators with two field strengths, we first note that summation over the generators of the adjoint representation of the dark gauge group may be needed to enforce gauge invariance. However, to keep things simple, we consider that after the dark sector symmetry breaking, only one (possibly complex) vector field is light enough to be produced at low energy, and simply discard the adjoint index. In addition, the nonabelian terms of the field strengths are systematically removed since the signatures we are after involve only two vectors. These restrictions, together with D µ V µν = −m 2 V V ν and D µṼ µν = 0, permit to reduce all the operators with derivatives acting on the field strengths (as well as those involvingQγ µ (

Reduction and phenomenology
Not all the operators are relevant phenomenologically. First, we can discard all those involving a derivative acting on the quark fields, since they are comparatively suppressed by O(m K,B /v). Then, consider the operators This is thus an explicit realization of the H V mat [I] operators, though with a major difference compared to the previous section. In the m V → 0 limit, the coupling ε V L,R here effectively scales at least like m 2 V instead of m V , and the P → P ′ V rates vanish when m V → 0 (compare with Eq. (31)). The precise scaling is not fixed though, because that between c ′V L,R and m V is not known. If the vector gains its mass at the scale Λ from some spontaneous symmetry breaking, all we can say is that m V ∼ gΛ for some g, and c ′V L,R ∼ g n for some n ≥ 0, is shown by the light yellow regions in Figs. 4 and 5. Comparing these regions with those corresponding to the first scenario (in blue), the rescaling (36) is very expensive in terms of accessible scales.
It is interesting to compare the reduction (36) to that of the magnetic operator of H V mat [II]. After integrating by part and using the quark EOM (21), Figure 4: Above: Exclusion regions drawn from a bound on the K + → π + V (left) and K L → γV (right) branching ratios at the 10 −10 level. The plots on the first line show the coupling g as a function of m V (in MeV) for the scenario I (blue, g = ε V ), scenario I with MFV (green, g = ε V |V * ts V td |), scenario II from the tensor operatorssσ µν d × V µν (yellow) and II from the vector operatorssγ µ d × ∂ ν V µν (light yellow), and scenario III (g = εe). The grey area represents the region where m V /g = v dark < v ≃ 246 GeV. The plots in the second line show the same, but replace g with Λ = m V /g (in TeV).
The first two terms match those of H V mat [I], Eq. (28) after electroweak symmetry breaking. Since we started from a gauge-invariant operator, gauge invariance is now hidden in the quark-mass dependent relationships among the couplings ε V L,R of H V mat [I], as well as of those of the higher-dimensional operators in Eq. (37). It is only upon imposing these relationships that all the terms originating from the k µ k ν /m 2 V piece of the polarization sum (29) cancel out. Though this is another physically sound interpretation of the H V mat [I] operators, it is much easier phenomenologically to consider directly the operator H †DI σ µν Q J × V µν of H V mat [II] from which they derive.
Note, finally, that a similar reduction starting with the last two operators of H VV mat [II] can also be done, Figure 5: Constraints on the coupling g against m V from B → KV (left) and B → K * V (right) for the scenario I (blue, g = ǫ V ), scenario I with MFV (green, g = ǫ|V * tb V ts |), scenario II from the tensor operatorsbσ µν s × V µν (yellow) and II from the vector operatorsbγ µ s × ∂ ν V µν (light yellow). The gray areas represent the regions where m V /g = v dark < v ≃ 246 GeV. Lower plots: Same as above, but in terms of Λ = m V /g (in TeV).
leading to the leading dimension-six non-gauge invariant operator of Eq. (32) As explained before, in the present work, we consider this operator exclusively in the gauge invariant form, since the dimension-four couplings of Eq. (28) dominate if no dark gauge invariance is enforced. After the reductions described above, the only operators relevant for rare FCNC decays are those involving the Higgs field, which we write after the electroweak symmetry breaking as (I > J) with, omitting flavor indices for simplicity: The rates and differential rates for K (B) decays are in Appendix A.5 (B.5), while the physics reach are in Tables 3 to 6. Also, the constraints on the tensor operators are shown by the dark yellow regions in Figs. 4 and 5, setting h T,T = (v/Λ)g and Λ = m V /g.

Flavor-blind operators
For the last scenario, the invisible vector boson is allowed to couple to conserved quark currents only, so that the singular term in the polarization sum (29) automatically cancels out in decay rate computations. To implement this, we have to drop the SU (2) L ⊗ U (1) Y gauge invariance requirement (see below), and couple V to quarks as The c q cannot be completely arbitrary but must reflect the flavor structures present in the SM. If that was not the case, then the FCNC couplings (28) of the first scenario are in general present [43], and those would completely dominate in B and K decays.
Using the MFV language, there are two possible flavor structures. Either the c q are universal or they are proportional to the Yukawa couplings. In the latter case, the top current would completely dominate, corresponding to Class III in the nomenclature of Table 2. In the former case, the remaining freedom corresponds to with c U and c D a priori arbitrary. Note that this two-parameter freedom means that it is always possible to write J µ [c q ] in terms of the electromagnetic and the baryon number currents only, as was pointed out in Ref. [44]. When the invisible vector is aligned with the photon, rather tight constraints are set by flavor-blind observables, e.g. the muon g − 2, quarkonium decays, beam dump experiments (see e.g. Ref. [4] for a recent analysis). While we will compare our constraints with those, let us stress that consistency does not require V µ to couple to leptons 7 . If it is leptophobic, many of these limits can be evaded. Hence, only the purely hadronic production of V in K or B decays will be considered here, and should be compared e.g. with those from π 0 or quarkonium decays with missing energy. As explained in Sec. 1, it is very different to probe the couplings to heavy (Class III) or to light quarks (Class IV), so we will analyze each situation in turn in the next two sections. But before that, let us return to the issue of an SU (2) L ⊗ U (1) Y gauge invariant origin for the conserved currents in Eq. (41).
In Sec. 4.1, the 1/m 2 V singularity of the decay rates arising from the FCNC couplings of Eq. (28) was interpreted in terms of a dark symmetry breaking scale v dark v. Specifically, we enforced that the coupling constant ε of V to quarks and its mass satisfy m V ∼ v dark ε, so that it is always the finite combination ε 2 /m 2 V → 1/v 2 dark which occurs in decay rates. But, generic vector or axial-vector quark currents are conserved for massless quarks, which is the case when SU (2) L ⊗U (1) Y is exact (for simplicity, we disregard the chiral current anomalies, which do not concern us here). So, the presence of the apparent 1/m 2 V singularity could be due to the electroweak symmetry breaking instead of to the dark symmetry breaking.
To make this statement more concrete (see also the discussion in Ref. [46]), let us write down the flavor-blind renormalizable couplings between V and SM gauge, Higgs, and quark fields (the flavor-blind summation over I = 1, 2, 3 is understood) These dimension-four couplings presumably dominate over the higher dimensional flavor-blind operators, which we do not list here (See Ref. [45]).
Only the ε 1 , ε 2 , ε θ , and ε B couplings are compatible with a dark sector gauge invariance associated to V , since they vanish under V µ → ∂ µ φ (see Eq. (31)) upon partial integration, and using the free Higgs boson and quark EOM. The ε B coupling involves the (conserved) baryon number current, and is thus directly matched onto Eq. (41). The other couplings, ε H , ε D , and ε U , would break the dark gauge invariance and are thus discarded. Note that ε H must in any case be tiny if V is to be light enough to be produced in rare decays. Indeed, the (H † H) × V µ V µ coupling gives a mass to V after the electroweak symmetry breaking, so barring a large cancellation between the dark and visible Higgs sectors, ε H m 2 K(B) /v 2 ≈ 10 −6 (10 −4 ). We can further discard the ε θ term, which is a total derivative and is relevant only for magnetic monopoles [47].
The two remaining couplings are ε 1 and ε 2 . The first one is the celebrated kinetic mixing [28]. Since The piece proportional to the (conserved) electromagnetic current ∂ µ F µν = −J em ν can be matched onto Eq. (41), while the other one mixes Z and V . The ε 2 coupling also generates a direct mixing Z µ ×V µ after the electroweak symmetry breaking, since The Z − V mixings induced by ε 1 and ε 2 are very different. The kinetic mixing Z ν × ∂ µ V µν is safe in the m V → 0 limit, because it is insensitive to the electroweak symmetry breaking. Actually, once the Z is integrated out, this Z ν × ∂ µ V µν vertex together with the SM flavor-changing hadronic vertex of Eq.
However, in parallel to the above FCNC operator, ε 2 also corrects the V mass as δm 2 V ∼ ε 2 2 v 2 . As a result, m phys V → 0 requires ε 2 → 0, ensuring again the safety of all decay rates in the massless limit. This shows that indeed, the visible symmetry breaking scale can play the same role as a dark symmetry breaking scale. The only difference, besides v = v dark in general, is that the former relies on the SM dynamics to drive the FCNC, and thus brings in the loop and CKM suppression factors, see Eq. (45). The allowed range of ε 2 values can be derived from the green areas in Figs. 4 and 5, up to the rescaling by g/16π 2 ∼ 10 −3 . Note that ε 2 values acceptable for rare decays ensure that δm 2 V ∼ ε 2 2 v 2 can be safely neglected. The above electroweak mechanism ensuring a safe massless limit may render the extension to a two Higgs doublet model desirable. Indeed, there would then be two different ε H couplings, and an additional conserved current can be constructed [37], whose charges are aligned with those of the Peccei-Quinn (PQ) symmetry [48]. Combined with the conserved lepton number current, this allows in principle to make V completely leptophobic [49]. The cost being the presence of a flavor-blind axial-vector quark current, not conserved at low energy. The corresponding 1/m 2 V singularity is nevertheless under control thanks to the additional sources of Z − V mixing present in the PQ current, i.e. tuned by the same coupling constant. We will not further elaborate on this construction, but just retain that a leptophobic setting is in principle possible, allowing to evade the many low-energy constraints based on theēγ µ e × V µ andμγ µ µ × V µ couplings [4,42].
Apart from the Z −V mixing effects, matched onto Eq. (28) or Eq. (34), the ε B and ε 1 couplings are genuine SU (2) L ⊗ U (1) Y invariant realization of the two conserved quark currents of Eq. (41). Let us now see how to constrain them from rare decays.

Phenomenological constraints on the couplings to heavy quarks
Once the heavy quarks are integrated out along with the weak bosons, the presence of V in K and B physics is felt through the operators of the second scenario, in particular H V mat [II]. For example, the last operator in Eq. (34) is induced in complete analogy to the electromagnetic operators describing b → sγ and s → dγ in the SM. This situation is thus simple to account by adapting the coupling of H V mat [II] according to Eq. (7), and setting the scale Λ at M W . Alternatively, a more precise estimate can be obtained when the new invisible vector boson is very light and aligned with the photon (in the quark sector). If we set c U = 2/3εe and c D = −εe/3, the branching ratios for b → sV and s → dV are obtained by rescaling by ε 2 the SM predictions for the b → sγ and s → dγ processes, up to simple phase-space corrections.
Specifically, in the B sector, the branching ratio for b → sV is when m V ≪ m B and E γ,V > 1.6 GeV. This cut on the photon energy is actually at the opposite end of phase-space compared to those set for b → sX. But even without a definite prediction, it is clear that the expected sensitivity of about 10 −5 in the B → (K, K * )X channels would at best probe ε down to a few percent. For comparison, typical bounds on ε derived from flavor-blind hadronic observables are currently down to the 10 −3 range [4,42].
The situation is worse in the K sector, where only CP-violating observables are sensitive to the shortdistance (c and t) magnetic operator. As analyzed in Ref. [51], those are beyond experimental reach even in the SM case, and thus cannot be used to set constraints on ε. This is actually clear from Table 5: rescaling by k sd ∼ 10 −6 , the scale Λ ends up well below the electroweak scale.
So, rare K and B decays are rather ineffective at constraining the presence of a new flavor-blind vector coupled exclusively to heavy quarks. Fortunately, in many cases, as e.g. from Eq. (44), universality holds and this vector must also couple to light quarks, where the situation is much better.

Phenomenological constraints on the couplings to light quarks
In this case, the CKM factors strongly favor the K sector to derive competitive bounds. At the K mass scale, only the u, d, and s quarks are active quark degrees of freedom. Adopting a matrix notation in the q = (u, d, s) flavor space, H V ef f [III] takes the form with 1 = diag(1, 1, 1), Q = diag(2/3, −1/3, −1/3), and e the QED coupling constant. So, from the point of view of low energy physics, there are only two possibilities: either V µ is effectively aligned with the photon (ε term) or its charges are proportional to baryon number (ε ′ term) [44]. This H V ef f [III] coupling must be directly embedded within Chiral Perturbation Theory (ChPT) [52]. At the leading p 2 order, the V µ field enters only through the covariant derivative acting on the meson fields The ε ′ term cancels out in the commutator, leaving V µ coupled exactly like the photon A µ . This ensures the absence of a direct K → πV coupling at leading order, relegating them to O(p 4 ). Such a direct leading order coupling only exists when the d and s charges are different. Indeed, in that case, the generator Q ′ would no longer commute with that of the weak interaction. This is another way to see that when the universality (42) fails, the dimension-four FCNC couplings of H V mat [I] should be allowed.  To get bounds on ε is rather immediate since the phenomenology is completely analogous to that of the radiative K decays (see Ref. [51] for a recent review). It suffices to consider the dominant radiative modes and replace one photon by V . When it is massless, the rates are obtained from those in QED simply by rescaling by ε 2 (or by α ′ /α, if one defines α ′ ≡ ε 2 α). When massive, the amplitudes are essentially the same, the polarization sum is identical (since the QED Ward identity holds), so the main impact is a reduced sensitivity due to the truncated phase-space. Assuming that about 10 13 kaon decays will be analyzed in the next generation of experiments, and bounds in the 10 −12 range are set, the reach in ε is thus naively for massless (or very light) vector boson V , and n + m > 0, n < 4. Since the current bounds on ε are down to the 10 −3 range, competitive bounds could be obtained from all the modes with B(K → nπ + (m+ 1)γ)) 10 −7 . Those are listed in Table 7.
The K → πV channel is special because K → πγ is forbidden due to gauge invariance. Further, even when off-shell, K → πγ * vanishes at leading order in ChPT, and so does K → πV . The leading contribution thus starts at O(p 4 ), from loops and local counterterms, and is approximately given by [53] A(K + (P ) → π + V (q)) = ε eG F 8π 2 a + q 2 P µ − q µ P · q ε * µ (q) , with a + an O(1) constant. The K S rate is expressed similarly in terms of the O(1) constant a S , while that for K 2 ≈ K L is CP-violating and thus driven by heavy quarks (since Im(V * us V ud ) = 0). From this amplitude, we get the rate with λ(1, r 2 π , r 2 V ) the kinematical function defined in Eq. (95), and r i = m i /m K . Normalizing with the K + → π + e + e − process to get rid of a + , with Φ e ≈ 0.145 the K + → π + e + e − phase-space factor, and B (K + → π + e + e − ) = (3.00 ± 0.09) · 10 −7 [12]. Reminiscent of K πγ, the rate vanishes in the m V → 0 limit. This seriously hampers the reach in |ε|, as indicated in Table 7. Even in the most favorable window m π m V 2m π , K + → π + V is less sensitive to |ε| than K L → γV by about an order of magnitude, see the red regions in Fig. 4.
To get bounds on ε ′ is more difficult because it cancels out from the O(p 2 ) Lagrangian, and thus also from the O(p 4 ) meson loops (see Eq. (48)). It can thus occur only in some local counterterms involving the V µν field strength (of vanishing anomalous dimension since there are no divergent loop contributions), and in the odd-parity anomalous O(p 4 ) Lagrangian. The former are very suppressed compared to the loop contributions induced by QED and by the ε piece, and will be neglected [51].
Concentrating on the odd-parity sector, only the K L → γV mode appears useful to constrain ε ′ (its anomalous amplitude, driven entirely by the up quark [54], is sensitive to both ε and ε ′ ). The magnetic direct emission amplitudes in K → ππγ are significantly more suppressed and difficult to access experimentally. Rescaling K L → γγ according to Eq. (49) shows that a bound on K L → γV at the 10 −12 level would probe couplings down to at least |ε, ε ′ | 10 −4 . Interestingly, this is more than an order of magnitude better than using the flavor-blind transition π 0 → γV , for which the best limit is 3.3 · 10 −5 , i.e. |ε, ε ′ | 6 · 10 −3 [12]. Further, the range of accessible V masses is evidently larger in K decays.

Baryon and lepton number violating operators
As for the dark scalar, Lorentz invariance requires an even number of SM fermion fields. However, the vector field index allows for alternative chiral structures compared to the Weinberg operators [2]. Specifically, keeping only operators of leading dimensions, where the tensor b µν stands for g µν or σ µν , the set a = QQQL, QQU E, DU U E, DU QL denotes the gaugesinglet combinations of fields of the Weinberg operators [2], the corresponding currents are defined as and where the SU (2) L triplet contraction for the QQQL current is understood. Note also that c II 1 is antisymmetric in flavor space, c I 3 is symmetric, while c I 1,2 have no particular symmetry, though one of them is redundant when flavor-diagonal. For the tensor current, which of the three Dirac structures does exist depends on the chiralities of the fermions involved (and may require reordering the fields using Fierz identities). The only other possible vector current is J µ D C ULL , but it vanishes upon Fierzing due to the antisymmetric SU (2) L contraction of the two lepton doublets and thus requires two more Higgs fields. The superscripts I and II refer to the scenarios discussed previously, i.e. separates those operators which explicitly break a dark gauge invariance associated with V from those which do not.
All these interactions have high dimensions, especially compared to the renormalizable couplings to SM particles of Eqs. (28) and (43). Phenomenologically, they may be relevant only when Λ is not too large, which requires the absence of direct FCNC couplings with V , see Figs. 4 and 5. At the same time, the ∆B = ±∆L interactions can induce nucleon decay when m V < m n (note that ∆B = −∆L interactions contain at least two d quarks and do not contribute to p + decay at the leading electroweak order), and thus require either the scale Λ to be extremely high, or the Wilson coefficients to have highly non-generic flavor structures [11].
A scenario with m B > m V > m n is interesting since astrophysical, leptonic, and nucleon decay bounds are essentially circumvented. In that case, rare B decays into an odd number of baryons (together with an invisible V ) may offer the best windows for the ∆B = ±∆L interactions. Note, though, that if these interactions occur concurrently to those of Eq. (28), (34), or (43), the V may occur as an intermediate state, bringing back the tight proton decay constraints. It remains to be seen whether in that case, signals in B decays are nevertheless possible. Such virtual exchanges are beyond our scope, since in the present work, we require the dark particle to be sufficiently long-lived to escape as missing energy in rare decays.
Finally, as for the dark scalar scenario, the ∆L = 2 interactions can produce three-body invisible ν L ν L V final states, but are not particularly interesting for FCNC decays. Indeed, they are unable to induce quark flavor transitions, and the FCNC decays would thus proceed through an extremely suppressed hadronic Higgs penguin.

Invisible spin-3/2 fermion
Spin-3/2 particles are described by Rarita-Schwinger fields, denoted Ψ µ , which transform as spinors with a vector index under the Lorentz group. The corresponding Lagrangian kinetic term can be written as [55] In addition, these fields are also subject to the conditions Ψ = 0 (spin-3/2 projection), (i/ ∂ − m Ψ )Ψ µ = 0 (Dirac equation), and ∂ µ Ψ µ = 0 (Lorenz condition). For external states, the spin summation is performed as [55] Π The sum over spin for v(p) s µ spinors is given by −Π(−p) µν . We distinguish the possible operators for the pair-production of these fields by their Lorenz structures, which can be scalar, vector, or tensor-like (even though Eq. (55) is written down for a Majorana field, Ψ will be taken as complex from now on). Taking into account the above-stated conditions reduces the possible leading operators to where Ψ [µ ΓΨ ν] = i(Ψ µ ΓΨ ν − Ψ ν ΓΨ µ )/2. Only the leading operators of each kind are kept; operators with additional derivatives are systematically discarded. For the vectorial couplings, the operators involving ε µνρσ Ψ ν γ ρ Ψ σ and ε µνρσ Ψ ν γ 5 γ ρ Ψ σ have been reduced to the others using the Chisholm identity. Finally, tensor structures are similarly reduced using the Chisholm identity together with σ µν ε µνρσ = −2iσ ρσ γ 5 , which in particular permits to get rid of the ǫ µνσρ Ψ σ γ 5 Ψ ρ and Ψ ρ σ µν γ 5 Ψ ρ structures.
The effective couplings of dimensions up to six involving gauge or Higgs fields are easy to construct from those in Eq. (11), After the electroweak symmetry breaking, the situation is similar as for spin 1/2 fields. The dimension-five operators in the first line generate a correction to the Ψ mass (upon enforcing the Ψ = 0 constraint), those in the second, third, and fourth line couple Ψ to the photon and to the Z boson. Note, though, that Ψ does not become millicharged in the usual sense, as these effective operators do not match those derived from the minimal substitution principle (which, in any case, is not consistent for spin 3/2 particles [56]). The last class is made of operators involving a single Ψ field. As for the spin-1/2 case, Lorentz invariance requires an odd number of SM fermion fields, so these operators break either baryon or lepton number: Notably, compared to the spin-1/2 case (15), no renormalizable coupling can be constructed, and thanks to the extra derivative in the dimension-five operator, there is (of course) no direct mixing between Ψ and ν L after the electroweak symmetry breaking. Phenomenologically, the signatures for the ∆B operator would again require specific searches in B decays, while those for the ∆L operators are to be found in semileptonic decays. Note, however, that the interactions in Eq. (59) are more difficult to access than those for spin 1/2 invisible particles, Eq. (15), because the tensor matrix elements vanish, 0|d I σ µν u J |P + = 0. This means that c ∆L 3 does not contribute to the P + → ℓ + Ψ decays, but only enters in the P → P ′ ℓΨ decays for which the helicityallowed SM contribution P → P ′ ℓν is large. Purely leptonic processes are only sensitive to higher-dimensional operators involving covariant derivatives acting on the quark or lepton fields, for exampleDD µ Q × Ψ µ L. So, the ∆L operators do not appear promising in searching for dark spin 3/2 particles, and will not be further considered here.

Reduction and phenomenology
The non-conserved quark flavor-changing neutral currents break the gauge symmetry appearing in the Lagrangian (55) when m Ψ → 0. As a result, the 1/m Ψ terms in the spin sum (56) are not projected out, the massless limit is singular, and we can force Λ up to arbitrarily high values simply by decreasing m Ψ . The situation is analogous to that encountered for the massive vector case in Sec. 4.1, and may resolve itself in a fully dynamical theory in a similar way. To get physically meaningful bounds on the scale Λ, there are two possible routes. The first procedure is inspired from the supergravity setting [57], where the spin 3/2 gravitino mass is related to the supersymmetry breaking scale as Λ SUSY = ( √ 3m Ψ M P lanck ) 1/2 with M P lanck = (8πG N ) −1/2 . In some sense, this can be understood as the fermionic equivalent of the constraint m V ∼ gv dark enforced for the vector bosons (which would here be insufficient given the harder (m −2 Ψ ) 2 singularities occurring when there are two spin-3/2 particles in the final state). Indeed, when Λ SUSY ≪ M P lanck , m Ψ is very small and only those terms originating from the m −2 Ψ singularity of the spin sum (56) are relevant [55], This projects out the ±3/2 helicity states, leaving the ±1/2 goldstino helicity states in a way similar to Eq. (31) for the massive vector boson. Specifically, the spin-3/2 operators become equivalent, in the m Ψ → 0 limit, to the spin-1/2 derivative operators obtained by replacing Given that the effective operators are at least of dimension six, there are at least two powers of Λ = M P lanck which can be eaten away by enforcing Λm Ψ → Λ 2 SUSY .
The supergravity scenario is thus characterized by the rescaling Λ SUSY = ( √ 3m Ψ M P lanck ) 1/2 . So, even if Ψ is here not necessarily identified with a light gravitino [58], let us assume that whereΛ may not be related to Λ SUSY in any way but could denote some dark sector symmetry-breaking scale. Phenomenologically, this rescaling permits to derive sensible bounds onΛ from the rare P → P ′ ΨΨ decays, even when m Ψ ≪ m P . A second route would start from a basis made entirely of gauge-invariant operators, i.e. involving the field-strength Ψ µν ≡ ∂ µ Ψ ν − ∂ ν Ψ µ and its dualΨ µν ≡ ε µνρσ ∂ ρ Ψ σ . Though the m Ψ → 0 limit would always be smooth, we do not perform this construction explicitly because with two field strengths, there are too many operators for the basis to be useful phenomenologically. Instead, we simply remark that starting from such a basis of gauge-invariant operators, it must be possible to generate the H ΨΨ mat operators by partial integration and use of the EOM, exactly as for the massive vector boson.
Specifically, when only the spin-3/2 EOM is used, iγ µ Ψ µν = m Ψ Ψ ν and γ µΨ µν = −m Ψ γ 5 Ψ ν , an extra factor m 2 Ψ /Λ 2 is generated. For example, This is analogous to the reduction (36) for the vector boson. By contrast, whenever the reduction involves the quark field EOM (21), the gauge invariance ends up hidden in relationships among the c i of H ΨΨ mat , exactly like in Eq. (37). The resulting operators are then suppressed either by the light quark masses, or by derivatives acting on the quark fields.
Enforcing some cancellations among the operators is incompatible with our procedure of turning on one operator at a time. So, we consider only the situation where the P → P ′ ΨΨ decay rates are finite in the m Ψ → 0 limit thanks to the rescaling withc i ∼ O (1). Note that this cures the singularity (60). Away from the strict m Ψ = 0 limit, the other terms of the spin sum (56) also contribute and tend to suppress the rates.
To derive the bounds on the scale Λ from the rare decays, we must impose one of the above two prescriptions (62) or (64) to make sense of the m Ψ → 0 singularities. Comparing them, these rescalings appear precisely equivalent for the dimension-six operators of H ΨΨ mat . For the dimension-seven operators, the gravitino-like rescaling (62) leads to an additional suppression by m Ψ /Λ compared to (64), making them completely irrelevant (remember that m Ψ is assumed smaller than m K,B ). Thus, the bounds we quote for these dimension-seven operators are understood to hold only for the second scenario.
As for the other types of invisible particles, we rewrite the various operators in terms of currents of definite C and P by introducing the fourteen complex couplings (for each s → d, b → s, and b → d operators) (65) The rates and differential rates are in App. A.6 (B.6) for K (B) decays, and the corresponding bounds are shown in Tables 3 and 4. As these numbers show, the rescaling prescription pushes the dimensionality of the operators to eight or nine. For such dimensions, the accessible scales Λ are at or even below the electroweak scale, as expected from Eq. (1) and  Table 8: References in the text for the various pieces of the basis of effective operators. Dimensions are denoted with a bar when the leading operator involves a Higgs field reducible to its vacuum expectation value after the electroweak symmetry breaking. In the last column, the dimensions are indicated for ∆B operators, irrespective of their ∆L components, while the dimensions in parenthesis are those of the purely ∆L operators. For dark vector fields, we distinguish between direct couplings, for which the m V → 0 limit is formally divergent, from those where a dark gauge invariance effectively survives (even though full invariance is not imposed, since we allow for m V > 0). Finally, the dimensions of the operators involving spin 3/2 fields increase when ensuring a sensible m Ψ → 0 limit, see Eqs. (62) and (64).

Conclusions
In this paper, we presented a complete basis of SU (3) C ⊗ SU (2) L ⊗ U (1) Y invariant operators involving SM fields together with a yet undiscovered light invisible spin 0, 1/2, 1, or 3/2 state, neutral under the SM gauge group. As summarized in Table 8, the operators are organized into three classes: couplings to SM fermions, couplings to SM gauge and/or Higgs fields, and baryon/lepton number violating couplings. We retained the operators of lowest dimensions separately for each class. As a result, most of them do not strictly qualify as portals since they are suppressed by the NP scale Λ. However, it makes sense to extend this denomination to those operators for which the experimental constraints push Λ far above the electroweak scale. For example, the typical scale for a dimension-five FCNC operator is greater than 10000 TeV (see Table 1), while it can even be close to the Planck scale for those inducing proton decay. For this reason, in the present paper, we systematically investigated the FCNC operators, and derived bounds from the rare FCNC transitions. Our results can be split into two parts: those concerning the basis of operators itself, and those related to its phenomenological impact on the rare decays. Starting with the basis, some of the main features are: 1. First and foremost, it must be stressed that even though we concentrated on the rare K and B decays involving missing energy to derive bounds on the operators, our basis is completely general and could be used equally well to investigate signals e.g. in lepton flavor violating transitions, flavor-blind quark or lepton observables, or at high energy colliders (see e.g. Ref. [46]).
2. For all spins, there is in the basis an operator involving only the SM Higgs field coupled to the dark state. Though specific models may not generate such operators, their presence would have two consequences. First, in general, the invisible state can no longer be naturally massless, though it can be very light, since a mass shift arises after the electroweak symmetry breaking. Second, the NP scale should always be greater than the electroweak scale, Λ > v ≈ 246 GeV, otherwise these corrections grow unchecked as powers of H † H are inserted (this is actually true for all our effective operators). It should be noted though that for spin 1 and spin 3/2 particles, enforcing a dark gauge invariance on the couplings to SM fields explicitly forbids such mass terms.
3. The contribution to FCNC transitions of the flavor-blind operators, i.e. either those involving SM gauge and/or Higgs fields or those involving SM fermion fields of the same flavor, has been clarified. Specifically, for scales Λ much higher than the electroweak scale, we pointed out that it is always advantageous to dress flavor-blind operators at the low scale with a W exchange (see Table 2). Indeed, such a low-scale GIM breaking does not necessitate additional Higgs tadpoles, which would each bring in a 1/Λ suppression.
4. For spin 1 and 3/2 dark particles, special care was devoted to maintaining a sensible massless limit, or at least to interpret the seemingly divergent limit. Indeed, the flavor-changing neutral currents are not conserved in general, so that the 1/m 2 V (1/m 2 Ψ ) term of the polarization (spin) sum is not projected out in physical observables. This is particularly relevant for massive vector states, for which renormalizable couplings to non-conserved quark currents can be constructed. Several mechanisms were discussed, and the corresponding NP scales derived from experimental bounds on the rare decay branching ratios were compared. As shown in Figs. 4 and 5, these scales strongly depend on the assumed dark sector dynamics.
5. All the leading operators producing a single dark fermion, whether of spin 1/2 or 3/2, violate either baryon (B) or lepton (L) number, but not both simultaneously, and have dimensions smaller or equal to that of the FCNC operators. By contrast, most of the ∆B and ∆L violating operators involving a dark vector or scalar particle directly derive from the dimension-six ∆B = ∆L = 1 Weinberg operators, or from the dimension-five ∆L = 2 operator. The former are negligible when they induce proton decay, i.e. when m φ,V < m p , but could induce exotic B decays into an odd number of baryons plus missing energy when m B > m φ,V > m p , since then proton decay is kinematically forbidden. The latter ∆L = 2 operators do not induce the quark flavor transitions, hence have a negligible impact on rare decays.
6. To implement the SU (2) L ⊗ U (1) Y invariance, the FCNC operators are constructed in terms of the SM chiral fermions. Though theoretically sound, and particularly convenient to implement the MFV flavor restrictions, such a basis is not convenient phenomenologically because many of these operators interfere in physical observables. So, our bounds are always derived using the alternative basis of operators obtained by projecting on currents of definite C and P . This minimizes interference terms, since in most cases these currents produce different final states.

7.
Although many examples of NP models involving new light states were mentioned, no attempt was made to precisely match them onto our basis of operators. Indeed, this would require dwelling into the detailed dynamics and parameters of each model, and would bring us too far from our main objectives, which were to construct the most general basis and constrain its operators from rare FCNC decays.
Concerning the rare FCNC decays of the K and B mesons, let us remind that provided the non-standard light states are neutral and sufficiently long-lived, they would show up as missing energy. Those modes cannot be experimentally distinguished from the SM processes producing a neutrino pair in the final state. However, these SM processes are very suppressed, and thus in principle, the rare decays could permit to identify even tiny NP effects. The main outcomes of our detailed phenomenological analysis of such effects are: 8. First, we stressed the importance of including the correct kinematical dependences for probing NP operators. This is crucial because experimentally, the rare decay modes with missing energy do not allow for a complete kinematical reconstruction, and require aggressive background suppressions. In practice, most experimental analyses implicitly assume at various stages that the differential rates have the shapes predicted by the SM, and this seeps through down to the final bounds on the branching ratios. Note, importantly, that these dependences are not always accounted for by simply enforcing the various experimental kinematical cuts. For these reasons, the specific kinematical dependences of each NP effect may have to be implemented by the experimentalists (those are detailed in App. A for K decays, and App. B for B decays). This provision has to be kept in mind when interpreting our bounds.
9. In the K sector, the sensitivities of the K L → X, K → πX, K L → γX, and K → ππX channels were compared, see Tables 3 and 5, with X a single or a pair of dark particles. The two-body K → X mode turns out to be the most sensitive, though for a very limited number of operators, but it is also the most difficult to deal with experimentally. At the other extreme, the K → ππX channels are sensitive to nearly all possible operators, but do not appear competitive given their phase-space and chiral suppressions. This leaves the K → πX and K L → γX channels, whose sensitivities to NP operators are in general comparable. Still, it should be noted that the latter, not yet considered experimentally, has some advantages. First, its SM contribution K L → γνν is at the 10 −13 level (see App. A.3.1), and thus cannot obscure even a tiny NP contribution. Second, for massive flavor-blind dark vector bosons or for millicharged fermions, as derived e.g. from a dark U (1) kinetically mixed with U (1) Y [28], the K L → γX mode is significantly superior to K → πX, whose relevant matrix elements vanish at leading order in Chiral Perturbation Theory.
10. In the B sector, the sensitivities of the B s,d → X, B → (K, K * )X, and B → (π, ρ)X decay channels were compared, see Tables 4 and 6. As for K, the fully invisible decays are both the most sensitive and the most difficult to probe experimentally. The main new feature compared to the K sector, besides an extended kinematical range allowed for m X , is that the modes with two light mesons in the final states can resonate, so that B → K * X, and B → ρX are competitive. This is particularly interesting since these modes are sensitive to all the quark currents but the scalarbs andbd. It should be said also that the present sensitivity of b → s and b → d decays, in terms of NP scales, is very similar. So, if the NP operators do not follow an MFV-like scaling, a small NP effect could be easier to identify in the latter.
11. The relative sensitivity of K and B decays was also compared. As expected, if the flavor structures of the NP operators involving X are generic, K decays are far more sensitive than B decays. However, it is well-known that in the visible sector, generic NP flavor structures are at odds with current experimental constraints. If MFV is imposed on both the visible and dark sector operators, the constraints from B decays become often tighter than from K decays (see Fig. 3), especially for chirality flipping currents q I (1, γ 5 , σ µν )q J , relatively suppressed by m s /m b , and for low-dimensional operators. Indeed, the impact of MFV on the scale Λ for an operator of dimension n decreases as n increases, since it is approximatively given by Λ MF V /Λ ≈ (V * tI V tJ ) 1/(n−4) for the d J → d I transitions, and with the CKM coefficients given in Eq. (6). Note, however, that n cannot be too large, since rare decay constraints give Λ v when n 8, see Table (1). In other words, for Λ v, the impact of such operators on the rare decays is beyond reach.
12. The ∆B and ∆L operators have low dimensions only for X = ψ or Ψ. The ∆B = 1 operators can only be probed with specific searches in B decays involving an odd number of baryons plus missing energy in the final state, and should certainly be included in future experimental programs. For the ∆L = 1 effects, the low-dimensional operators are accessible only for X = ψ, which contribute to P + → ℓ + ψ (P = K, D, B), and would thus apparently enhance the purely leptonic transitions P + → ℓ + ν. If the flavor structure of these NP operators is non-universal, this could resolve the persistent discrepancy in B → τ ν while remaining consistent with the B → (e, µ)ν bounds, as well as, if m ψ < m K , with the tight lepton universality constraint derived from K ℓ2 decays, see Eq. (17).
In conclusion, the presence of a light invisible state weakly coupled to SM particles is not only far from excluded, but is even compelling in many NP models. To find such states, a host of experimental facilities currently available or in planning are called in, from high-intensity meson and lepton factories to high energy colliders, neutrino detectors, earth or space-based direct or indirect dark matter searches, high intensity lasers,... In this big picture, the very rare FCNC decays of the K and B mesons, with their unique sensitivities and kinematical ranges, could play a crucial role in the very near future thanks to the leap in luminosity expected at the next generation of experiments, namely NA62 at CERN and K0TO at J-Parc dedicated to these K decays, and Super-B in Italy and Belle II at KEK aiming for the B decays.
A Differential rates for K decays A.1 Experimental observables in rare K decays Let us start by reviewing the kinematics and current experimental limits for the various K decays induced by neutral currents and involving missing energy.
A.1.1 K → π+ missing energy When the missing energy consists of two invisible particles, the differential rate depends only on the invariant mass of these particles, z ≡ q 2 /m 2 K , or equivalently, on the pion momentum P π ≡ |p π |/m K = √ λ/2, with λ ≡ λ(1, z, r 2 π ) defined in Eq. (95), r π ≡ m π /m K . The phase-space integral is then with r X = m X /m K . In the SM, the only available invisible particles are the neutrinos. The SM spectrum for K + → π + νν and K L → π 0 νν then derives entirely from the vector current matrix element π|sγ µ d|K , and involves the corresponding form-factor (see Eq. (81) in the next Section) slopes λ ′ + and λ ′′ + : where i = +, 0. Translated in terms of the pion momentum, i.e. using z(P π ) = 1 + r 2 π − 2 r 2 π + P 2 π , this becomes The slopes λ ′ + and λ ′′ + are conventionally normalized by the charged pion mass, and are equal for the K + and K 0 decays to an excellent approximation (see Ref. [59]). They are extracted from K ℓ3 decays as λ ′ + = r λ (24.82 ± 1.10) · 10 −3 , λ ′′ + = r λ (1.64 ± 0.44) · 10 −3 , r λ = 0.990(5) , with r λ a rescaling factor accounting for the K * + − K * 0 mass difference. These (highly correlated) errors are negligible compared to the experimental and theoretical errors on the integrated rate, and are neglected in Table 9. Currently, only the charged decay has been observed at Brookhaven [60], in two momentum windows separated by the K + → π + π 0 peak, and with the lower end corresponding to the K → πππ threshold (see Table 9 and Fig. 6). The proposed charged K experiment at J-Parc would use the region above the K + → π + π 0 peak [61], while the two windows planned at NA62 [62] and proposed at Fermilab [63] are similar. Note that these experiments use very different techniques (stopped vs. in flight), but in both cases, the momentum of the initial and final charged particles are in principle accessible. It is important to stress that not only the combination of the measurements done for each specific window (see Table 9) assumes the SM spectrum, but also that within each window.
For the neutral mode, the K S → π 0 νν mode is CP-conserving but difficult to access given the very short K S lifetime, so we concentrates on the CP-violating K L → π 0 νν mode. The best limit [64] B(K L → π 0 νν) < 2.6 · 10 −8 , was obtained by the E391a experiment at KEK, and will be further improved using the same techniques at J-Parc [61]. In these experiments, the K L momentum is not fixed. So, the high hermeticity of the detector is essential to ensure sufficient suppression of the backgrounds. Though less effective in this case, kinematical |p π | (M eV )  Table 9: Experimental measurement of K + → π + νν and SM prediction within each momentum window [60], in units of 10 −10 . Figure 6: The experimental windows in π + momentum used for controlling backgrounds in the K + → π + νν measurements, with the seven events seen at Brookhaven. The SM spectrum corresponds to a vector couplinḡ sγ µ d ×ν L γ µ ν L , and is implicitly implied in computing the branching ratios from the events.
cuts are still useful. In particular, the transverse momentum P T of the reconstructed π 0 is required to be large, between 120 and 240 MeV. This does not cut away the background from K L → π 0 π 0 , but rather ensures that the two extra photons have high momentum, and are thus difficult to miss. Since the momentum spectrum of the π 0 cannot be directly measured at KEK or J-Parc, and since the SM decay spectrum is implicitly assumed in the analysis, it is far from immediate to translate the current limit (70) into bounds on non-standard currents involving other types of invisible particles. So, for both the charged and neutral modes, it is not currently possible to deconvolute the SM spectrum from the experimental numbers. To proceed and derive the bounds quoted in the text, we require that the predicted branching ratio for the production of new invisible states does not exceed 10 −10 when integrated over the momentum windows of the charged mode. This is a rather loose approach, which could significantly underestimate the experimental reach in case the spectrum is very different than the SM one. To illustrate this, note that the bounds for two-body decays are already slightly tighter [65], To improve our naive bounds, either the specific modulations of the spectrum in the presence of NP have to be included throughout the experimental analysis 8 , or the true momentum spectrum must be measured (maybe using TOF techniques for the neutral mode [63]). Finally, independently of the NP spectrum, it should be noted that the sensitivity to the K → πX(X) processes is ultimately bounded at around 10 −12 by the theoretical error on the SM predictions for the K → πνν branching ratios.

A.1.2 K → ππ+ missing energy
The phase-space integration for the K → π(K 1 )π(K 2 )X(q) decays is with y = (K 1 + K 2 ) 2 /m 2 K the invariant mass of the pion pair. Compared to K → πX, these modes are suppressed by the smaller hadronic matrix elements and by phase-space. Further, the kinematical range is much reduced. Currently, the best limits are (see the respective papers for different m X values) and are thus very far from Eq. (71). Note that the hadronic matrix elements π + π 0 |sΓd|K + , π 0 π 0 |sΓd|K 0 , and π + π − |sΓd|K 0 are related in the isospin limit, see Eq. (94) below, so that these experimental constraints suffice to completely bound the K → ππX system. The K (P ) → π(K 1 )π(K 2 )X(p 1 )X(p 2 ) decays are similarly suppressed, and the experimental information is less precise. The phase-space integrals reduce to that over the invariant mass of the invisible pair T 2 = (p 1 + p 2 ) 2 = zm 2 K and of the pion pair K 2 = (K 1 + K 2 ) 2 = ym 2 K , over the range Here again, the SM spectrum critically enters and the current experimental bound cannot immediately be translated into bounds for invisible particles of a different type. Though probably optimistic given the limited phase-space and complicated signatures, we assume bounds of 10 −10 on each of these K → ππX and K → ππXX modes are achievable to derive the numbers in Tables 3  and 5. In any case, if the bounds are different, it is a simple matter to rescale the numbers accordingly. Further, using the same 10 −10 branching ratio bounds as for the K → πX(X) modes permits to clearly illustrate the reduced sensitivity of the K → ππX(X) channels.

A.1.3 K → γ+ missing energy
Compared to K → ππX and K → ππXX, the modes K (P ) → γ(k)X(p 1 )X(p 2 ) and K (P ) → γ(k)X(T ) are less suppressed and could offer simpler experimental signatures. The phase-space integral for the three-body decay is with z = (p 1 + p 2 ) 2 /m 2 K the invariant mass of the invisible pair. Despite their theoretical sensitivity, there is currently no experimental limit on these modes. So, for now, we assume that the next generation of experiments will reach B(K L → γX) < 10 −10 . Note that with about 10 12 − 10 13 K L decays, as required to measure K L → π 0 νν, this may be pessimistic. Further, while the K → πνν processes ultimately limits the sensitivity to K → πX(X) at a few 10 −12 given the current theoretical errors, the SM rate for K L → γνν is at the 10 −13 level, so bounds at or even below that level are in principle achievable.
Other modes with photon and missing energy will not be considered, as the K → nπ + γ + X processes are either too suppressed and difficult to access experimentally when fully neutral (given the many photons from the π 0 s), or superseded by the non-radiative processes K → nπ + X when some mesons are charged (since the photon of K → nπ + γ + X is essentially a bremsstrahlung radiation off the charged meson [68,69], the amplitude for K → nπ + γ + X is actually driven by that for K → nπ + X [70]).

A.1.4 K → missing energy
The simplest decays are those where the K L or K S simply disappear. Though difficult to probe experimentally, the simpler matrix element together with the minimal number of final state particles strongly enhance their sensitivity to NP effects. In the case of the π 0 → XX process, the best bound is [71] B(π 0 → XX) < 0.27 · 10 −6 .
But this measurement is actually a by-product of the study of the K + → π + XX decay, since a bound on K + → π + XX indirectly constrains K + → π + π 0 [→ XX]. Doing the same for K L → XX would require very tight bounds on some B or D decays with missing energy, well beyond current capabilities (see Sec. B.1). Alternatively, a direct bound on K L → XX may be obtained from φ factories, where the other K can be tagged. In deriving NP scales in the text, we will use B(K L → XX) < 10 −10 to simplify the comparison with the other modes, but it should be kept in mind that such a bound appears extremely challenging.

A.2 Matrix elements for K decays
In the K sector, the quark currents are represented within Chiral Perturbation Theory (ChPT) [52]. For simplicity, only the leading chiral order is kept. Specifically, the vector and axial-vector currents start at O(p): with, at leading order, F = F π = 92.4 MeV. Thanks to the QED gauge invariance, there is no unknown low-energy constant in these currents. The scalar and pseudoscalar currents start at O(p 0 ): with the low-energy constant B 0 related to the quark masses, with r π = m π /m K , and using m s (2 GeV) = 100 ± 20 MeV (so when deriving bounds on c i /Λ n , c i ≡ c i (2 GeV) is understood) [12]. The scalar and vector currents are related by the EOM. Specifically, the most general matrix elements for the scalar or vector K → π transitions have the form (z = q 2 /m 2 K , q = P − K) up to some simple Clebsch-Gordan coefficient. Taking the divergence of the vector current produces q 2 f − (z) = (m 2 K − m 2 π ) (f 0 (z) − f + (z)). So, at the leading chiral order, f +,0 (z) = 1 and f − (z) = 0. Refinements are only needed for a precise prediction of the SM rates, but are not numerically relevant for the bounds on the production of new invisible states. Note that in practice, the scalar and vector currents do not need to be parametrized as external couplings, but can be directly introduced through the ChPT source terms. Doing this using the leading O(p 2 ) Lagrangian reproduces Eq. (78) and Eq. (79). We do not consider the next-to-leading O(p 4 ) meson loops and local terms, except for the odd-parity contact interactions obtained by introducing the vector and axial vector sources in the O(p 4 ) anomalous Wess-Zumino-Witten (WZW) action [52]. Indeed, owing to their opposite parity, these interactions drive the leading order contributions for some amplitudes.
Finally, from Lorentz and chiral symmetry, combined with parity and charge conjugation (valid for the strong interactions), the most general parametrization for the tensor current starts at O(p 2 ), where it is given by [51]q Two new low-energy constants a T and a ′ T occur, for which we use the Lattice estimates (see the discussion in Ref. [51]) B T (2 GeV ) = 2m K a T = 1.21(12) [72] , B ′ T (2 GeV ) = 2F π a ′ T = 0.6(2) [73] .
Note that a more recent lattice estimate B T = 0.65(2) [74] is two times smaller, and thus suppresses the sensitivity of K decays to the tensor currents. Still, in terms of the NP scales Λ given in the Tables 3 and 5, the precise value of B T is not that relevant at present since these numbers are to be understood as order of magnitude estimates. So, altogether, and defining Γ ≡ c S + c P γ 5 + c V γ µ + c A γ µ γ 5 + c T σ µν + cT σ µν γ 5 , the matrix elements in the isospin limit and to the leading chiral order are 9 and − π + (K 1 )π 0 (K 2 )|sΓd|K + (P ) = M − (K → ππ) , √ 2 π 0 (K 1 )π 0 (K 2 )|sΓd|K 0 (P ) = M + (K → ππ) , with where Terms proportional to the number of QCD colors, N C = 3, come from the WZW action. The m 2 K − T 2 − denominators arise from the kaon pole topologies, K → ππK 0 followed by K 0 → 0 (from Eq. (84)). Note that M + (K → ππ) is even under K 1 ↔ K 2 , while M − (K → ππ) is odd, hence these amplitudes describe two-pion states with even and odd orbital angular momentum, respectively. The π 0 π 0 state is purely even due to Bose statistics, while the π + π 0 state has total isospin one, hence is purely odd.
Only the tensor currents contribute to the M(K → γ) amplitude at tree-level. So, it may seem that together with the O(p 4 ) WZW amplitude, we should include also the even-parity meson loops along with their counterterms. However, for neutral current sources, i.e., in terms of Gell-Mann matrices, for v µ , a µ , s, p ∼ λ 6 ± iλ 7 , the only allowed even-parity K → γ matrix elements vanish when the photon is on-shell γ(q, ν)|sγ 5 where T = P − q, and Φ(q 2 , m 2 ) the loop functions occurring for K → πγ * (see e.g. Ref. [68]), defined in terms of the standard scalar one-loop integrals as Contrary to K → πγ * , the FCNC matrix elements are finite since the UV divergences cancel in the difference between the K ± and π ± loop contributions, and there are no counterterms. Though in principle, the K → ℓ + ℓ − X modes could thus offer precise probes, their rates are far too suppressed by α, the loop factors, and the cancellation between the K ± and π ± loops. So, only K → γX will be considered here, which is thus induced exclusively by tensor currents and O(p 4 ) anomalous interactions. Finally, the operator basis also includes vector and scalar currents with a covariant derivative. Though these operators are never retained in deriving bounds on the new physics scale in the main text, for the sake of completeness, let us nevertheless write down the relevant chiral realizations. For the vector currents, extending the analysis of Ref. [75], we writē Most of the terms are fixed by taking divergences and imposing the EOM, but for the constant a V , a priori of O(1). For the scalar currents, the chiral realizations start at O(p), These currents are completely fixed by imposing parity and charge conjugation, together with Note, however, that the above chiral representations lead toq . So, terms at that order would be required to get a correct divergence.

A.2.1 Decay rates in the isospin limit
The strong matrix elements (84)(85)(86) are derived in the isospin limit. As a result, all the differential decay rates can be reconstructed entirely from those of the K L ≈ K 2 . The contributions coming from the εK 1 piece of the K L , suppressed by ε ∼ 2 · 10 −3 , are neglected here.
For the K → (γ)X modes, the K S ≈ K 1 rates are obtained from those for K L ≈ K 2 by interchanging the real and imaginary parts of the couplings x = f i , g i , h i : The K + → π + X decay rates are proportional to the sum Γ(K S → π 0 X) + Γ(K L → π 0 X), hence are obtained from Γ(K L → π 0 X) through the substitutions (x, Finally, the whole set of K → ππX decay rates can be reconstructed from the K L → π + π − X rate as follows. Denoting Γ (K L → π + π − X) ± ∼ |M ± (K L → ππ)| 2 from Eqs. (84,85), and noting that these two amplitudes do not interfere, we get In the following sections, the differential rates are given for the various scenarios adopting the notations of Eq. (66) for K → πXX modes, Eq. (74) for K → ππXX modes, Eq. (72) for K → ππX modes, and finally, Eq. (76) for K → γXX modes. For K → πX and K → γX, the total integrated rate is directly written down. Given their regular occurrence, let us also introduce specific notations for the usual kinematical functions. First, λ αβ = λ(1, α, β) , λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc) , with α,β standing for the reduced variable y, z, or the reduced mass r π , r X (in which case we simply denote α, β = π, X). Similarly, we define

A.3 Spin 1/2 invisible particles in the final states
The rate for the fully invisible decay is: The differential rates for the decays into a pion plus invisibles are with the normalization and kinematical functions In the massless limit, this expression simplifies a lot because the fermion helicity states do not mix. The interference terms drop out while the parity of the current becomes indistinguishable; J i = J ′ i , i = 1, 2, 3. The decays K → ππψψ also receive contributions from all the currents: with the normalization and kinematical functions Again, the interference terms drop out and F ′ i = F i , i = 1, ..., 6, when m ψ → 0. The rate K L → γψψ is driven either by the anomalous vertices at O(p 4 ) for vector currents, or directly by the tensor currents, and has the differential rate where N C is the number of QCD colors. The O(p 4 ) loops and normal-parity counterterms cancel out for the (axial-)vector and (pseudo-)scalar currents when the photon is on-shell, see Eq. (87).
The corresponding branching ratios are B K L → π + π − νν = 1 · 10 −13 , B K L → π 0 π 0 νν = 6 · 10 −14 , B K + → π + π 0 νν = 5 · 10 −15 , (107b) Adopting the usual chiral counting, the typical error on these LO estimates is expected to be of about 30% at the amplitude level. Note, in this respect, that the K → πνν rates given above are just indicative, as higher order corrections as well as isospin-breaking effects are known and more precise estimates have been obtained [59]. For K → ππνν, the anomalous term gives negligible percent-level contributions, so that ℑ(y ν ) < ℜ(y ν ) implies B (K L → π + π − νν) ≈ 2B K L → π 0 π 0 νν ≫ B K + → π + π 0 νν , in fair agreement with Ref. [79]. For K L → γνν, previous estimates incorrectly rely on K ℓ2γ for the matrix elements [80], hence included a parityeven contribution in contradiction with Eq. (87). Numerically, this is however without consequences since these contributions are CP-violating hence subleading for the rate. Finally, the rates for K L → ℓ + ℓ − νν are dominated by the Dalitz emission from the purely anomalous K L → γνν. The rates from the even-parity contributions arising from the matrix elements Eq. (87) are in the 10 −19 range. Note that we did not consider the tree-level process , which may actually be competitive given these strong suppressions. In any case, these rates are far too small to be accessible any time soon.

A.4 Spin 0 invisible particles in the final states
The rates for the production of a single invisible scalar from the H † (DQ)φ operator of Eq. (20) are The matrix element γ|sd|K L occurring for K L → γφ vanishes at O(p 4 ). The rates corresponding to the derivative couplings (sγ µ d)∂ µ φ and (sγ µ γ 5 d)∂ µ φ are obtained through the replacement (22), i.e.
Note that if simultaneously present, the g S,P and g V,A currents obviously interfere. The differential rates for two-scalar final states can be written in terms of the same kinematical functions as for fermionic final states, but for the obvious replacement m ψ → m φ everywhere. They are significantly simpler though, because the angular momentum of the two-scalar states is purely orbital. Specifically, with the kinematical quantities defined in Eqs. (97, 99, 102, 103). Note that if φ =φ, Bose statistics has to be enforced and these rates should be divided by two.

A.5 Spin 1 invisible particles in the final states
From the single dark vector production from the operators of H V mat [I] and H V mat [II], Eqs. (28) and (34), we find , Note that the H ′ 3 term is finite when r V → 0 since it is induced by the WZW anomaly, while by construction, all the Γ [II] are finite in that limit. For the two-vector modes induced by the scalar and pseudoscalar couplings of H VV mat [II], the rates are Γ [II] K L → VV = 4Γ 0 m 6 K Λ 6 I 1 (ℜ(h P P ) 2 + ℑ(h P S ) 2 ) + 6r 4 V I ′ 1 ℑ(h P S ) 2 , dΓ [II] dz K L → π 0 VV = 4Γ π 0 m 6 K Λ 6 zJ 3 (ℜ(h SS ) 2 + ℑ(h SP ) 2 ) + 3r 4 and the kinematical quantities defined in Eqs. (97,99,102,103) but for m ψ → m V everywhere. The matrix element γ|sd|K L occurring for K L → γV V vanishes at O(p 4 ). Note that if the vector field is real, V =V , then these rates have to be divided by two to account for Bose statistics.
with the definitions in Eq. (103) together with polarization of the K * . The experimental information that can be obtained from the process B → K * (→ Kπ)X with an on-shell K * is completely described by the double differential decay distribution in terms of the two kinematical variables s = (p B − p K ) 2 , corresponding to the invariant mass of the final state invisible particles, and θ, the angle between the K * flight direction in the B rest frame and the K flight direction in the Kπ (K * ) rest frame. The spectrum can be expressed in terms of B → K * transversity rates Γ L,T corresponding to longitudinally and transversely polarized K * final states (see App. B.1.4), while the double differential spectrum can be written as d 2 Γ ds d cos θ = 3 4 dΓ T ds sin 2 θ + 3 2 dΓ L ds cos 2 θ .
Thus, dΓ L /ds and dΓ T /ds can be extracted by an angular analysis of the K * decay products. Finally, the total invisible mass distribution is dΓ ds As an illustration of the potential impact of such future precision measurements, in Tables 4 and 6, we also present the accessible NP scales assuming a 20% relative precision on the B → K ( * ) X rates compared to their SM predictions.
A potential NP signal in these modes would manifest itself via differences in the measured B(B → τ ν) values in the leptonic and hadronic decay modes of the τ -something that future dedicated experimental searches could employ to reduce systematic uncertainties. Although the SM signal shape in this case is well determined and largely free from theoretical form factor uncertainties, the appearance of two neutrinos in the final state means that the same experimental caveats in extracting bounds on NP contributions apply as in the B → K ( * ) X case [84]. Again the polarization states of the ρ in the B → ρX mode, which are well predicted within the SM, could be reconstructed using the angular distributions of the ρ → 2π system and aid in discriminating SM contributions from possible NP effects.
The neutral modes B 0 → π 0 (ρ 0 )X are free from long distance SM contributions, but the purely neutral final states make them more challenging experimentally. The present bounds of B(B 0 → π 0 νν) < 2.2 · 10 −4 [83] , B(B 0 → ρ 0 νν) < 4.4 · 10 −4 [83] , are less constraining than the charged modes analysis. Consequently, in setting our bounds on invisible particles in the massless limit in Tables 4 and 6 and away from this limit in Figure 2 we tentatively allow NP contributions to saturate the experimental uncertainties in the charged modes, Eq. (126).

B.1.3 Other modes with missing energy
Both Belle [85] and Babar [86] have searched for the B → XX decay mode, which is helicity suppressed in the SM. While unresolved soft photons can partially lift this suppression [6,87], the SM predictions for the branching ratio remain at the order of 10 −9 [6]. Being a two body decay process in the scenarios we consider, the kinematical distributions are trivial and no model-dependent efficiency corrections are needed. The latest experimental upper limit reads B(B → XX) < 1.3 · 10 −4 @ 90% C.L. [85] , and can be employed directly to constrain the relevant interactions of invisible particles, as given in Tables 4  and 6 and Fig. 2. In the future, both B s,d → XX modes could potentially be probed to greater accuracy at the super flavor factories [88].
On the other hand, while B s,d → γX with reconstructed final state photons allow to lift the helicity suppression suffered by the NP fermionic contributions proportional to f AA , the additional α EM suppression only makes them competitive with the B → ρ(K * )X modes in the opposite region of large invisible particle masses (m B − m ρ(K * ) < m X < m B ), where the helicity suppression is least and the B s,d → X channels are effective in constraining NP effects as well. In addition, precisely in this region the SM backgrounds, most notably from misreconstructed B q → γX events, limit the experimental reach of B s,d → γX [86]. For a more detailed discussion of the expected NP sensitivity, we refer to Ref. [6] and do not consider these modes any further.

B.4 Spin 0 invisible particles in the final states
The rates for the production of a single invisible scalar from the H * D Q φ operator of Eq. (20) are where Γ BH ( * ) 1i = m B λ 1/2 H ( * ) i /16π, λ ij =z ≡ λ(1, r 2 i , r 2 j ). The vector operator contributions are related to these via quark EOM, see Eq. (22). The rates for the production of two scalars are where Note that if φ =φ, Bose statistics has to be enforced and these rates should be divided by two.

B.5 Spin 1 invisible particles in the final states
The production of a single vector using the simple FCNC operators of Eq. (28) or the gauge-invariant operators of Eq. (34) are B → HV : where Note that in the m V → 0 limit, Γ [II] L (B → H * V ) is well-defined while the tensor operator contributions actually vanish thanks to the form factor relation T H L (0) = 0. The rate and differential rates for the production of two vectors are I P P ψ ′ |f P P | 2 + I P S ψ ′ |f P S | 2 − I AA,P P