Flavour anomalies from a split dark sector

We investigate solutions to the flavour anomalies in $B$ decays based on loop diagrams of a"split"dark sector characterised by the simultaneous presence of heavy particles at the TeV scale and light particles around and below the $B$-meson mass scale. We show that viable parameter space exists for solutions based on penguin diagrams with a vector mediator, while minimal constructions relying on box diagrams are in strong tension with the constraints from the LHC, LEP, and the anomalous magnetic moment of the muon. In particular, we highlight a regime where the mediator lies close to the $B$-meson mass, naturally realising a resonance structure and a $q^2$-dependent effective coupling. We perform a full fit to the relevant flavour observables and analyse the constraints from intensity frontier experiments. Besides new measurements of the anomalous magnetic moment of the muon, we find that decays of the $B$ meson, $B_s$-mixing, missing energy searches at Belle-II, and LHC searches for top/bottom partners can robustly test these scenarios in the near future.


Introduction
Several flavour anomalies have been observed in the last few years in various B-meson decay measurements by different experimental collaborations. Some of the anomalous measurements involve semileptonic b → s transitions and include: (1) the suppression with respect to the Standard Model (SM) expectation of the ratios R K and R K * -the branching ratios of the B-meson decay into a K or K * meson and muons, over the decay to the same kaon and electrons -which were observed at LHCb [1][2][3][4] and which imply the possible violation of lepton-flavour universality (LFUV); (2) an enhancement in the angular observable P 5 , first measured by the LHCb [5] and Belle collaborations [6], and later observed also by ATLAS and CMS [7,8]; and (3) a suppression in the branching ratios for the decays B s → φµ + µ − [9] and B → K ( * ) µ + µ − [10,11].
The solutions based on a light dark sector can be divided roughly into two categories, depending on whether the masses involved lie above or below the characteristic window for LFUV observables, which is commonly identified as ranging roughly from the dimuon threshold to √ 6− √ 8 GeV. The first category involves a new light vector with mass m V 2 GeV, coupled to the b − s and the muon currents, interfering negatively with the SM amplitude [18,33]. The light particle can contribute strongly to the physical observables thanks to the proximity of a resonance to the experimental bins [33] and NP effects can be parametrised in this case by Wilson coefficients with an explicit q 2 dependence. Interestingly, the corresponding resonance in the dimuon spectrum from B-meson decay could be hidden due to the sizeable hadronic uncertainty and the presence of the J/Ψ resonance in the same region [33,37]. A second class of models, featuring instead the exchange of light particles below the dimuon threshold, have been considered in refs. [18,32,34,35]. These scenarios are strongly constrained by a number of observations. On the one hand, a light scalar particle with effective couplings to the s and b quarks and leptons yields a positive contribution to the decay rate. One thus requires a sizeable coupling to electrons, a scenario that is in most cases [34] in tension with the measurement of B → K ( * ) e + e − processes in agreement with the SM [38]. On the other hand, a light vector boson with effective couplings to the s and b quark and muons can interfere negatively with the SM process, but is strongly constrained by the measurement of the anomalous magnetic moment (g − 2) µ . This in turn implies a sizeable coupling to the b and s quarks, leading to strong bounds from B s mixing. Finally, in ref. [35] it was pointed out that a resonance coupling to electrons can actually be used to reproduce the low-q 2 bin of R K * . A vector boson very close -but below -the dimuon threshold, where BR(B → K ( * ) V ) can be as small as 1 × 10 −7 , can explain the anomaly and escape the limits from ref. [38].
A common trait of the constructions mentioned above is the presence of an effective offdiagonal coupling g bs to the quark current. Since the quarks carry SU(3) c charge, this must arise from particles with colour integrated out in the UV theory, which must be heavy to avoid exclusion by the strong limits from the LHC. This also means that in realistic scenarios g bs will be suppressed by loop effects and/or a small mixing angle between the SM quarks and these heavy particles.
In this paper we explore in detail the concrete possibility that the flavour anomalies emerge from UV-complete models giving rise to the effective b−s coupling via loop effects involving the particles of both the heavy and light dark sectors. We will show that, under these assumptions, several phenomenological challenges must be faced to generate a g bs large enough to provide a good fit to the LFUV observables once the constraint from (g − 2) µ is taken into account, and at the same time maintain a reasonable agreement with perturbativity of the couplings.
In the specific cases mentioned above, we will show that solutions that belong to the first category, m V 2 GeV, can be made viable with appropriate UV completions, even though they are subject to a combination of increasingly tightening bounds that include, in the UV, LHC direct constraints on hadronic and leptonic new states and Drell-Yan dimuon constraints from the Z lineshape, and in the IR, an appropriate width-and bin-dependent treatment of the bounds from BR(B → K + invisible) and BR(B → Kµµ). Conversely, complete solutions that belong to the second category mentioned above, characterised by m V under the dimuon threshold, are much less likely to survive, as they are strongly constrained by a combination of bounds from BR(B → K + invisible), neutrino trident production, intensity frontier limits on kinetic mixing, and BR(B s → µµ).
More in general, we perform in this work a detailed Monte Carlo scan of the broad m V range and the loop-induced couplings. We identify and characterise the specific properties of the viable models and provide an indication of possible strategies for a timely detection of the associated NP states. The paper is organised as follows. In Sec. 2 we recall expressions for the loop-generated Wilson coefficients from box and penguin diagrams. We provide the quantum numbers of the particles entering the loops and estimate the characteristic size of the couplings required to fit the flavor anomalies. In Sec. 3 we enlist the complete set of constraints we apply to the models. Section 4 is dedicated to the main results, with a description of the fitting procedure and discussion of the allowed parameter space. We present our conclusions in Sec. 5. Appendix A is dedicated to the detailed treatment of the recasting procedure for B → K + inv. limits.

One loop Wilson coefficients from split dark sector models
Our goal is to investigate the extent as to which the flavour anomalies can be explained by generic loop effects involving light new particles, in association with a heavy sector above the EWSB scale. Possible constructions for the Wilson coefficients of the effective Hamiltonian descending from loops of TeV-scale new particles have been investigated, e.g., in refs. [39][40][41][42][43][44][45][46]. They involve either box diagrams of scalar and fermionic states like in Fig. 1(a), or penguin diagrams like the ones represented in Fig. 1(b).

Box diagrams
In the presence of direct Yukawa couplings between the quarks b, s, the muons, and a NP sector composed of fermions ψ i and scalars φ j , the relevant Lagrangian is given by where a sum over repeated indices is implied.
Since at least one of the fermion or scalar fields in the box diagram must necessarily carry the colour charge, the bounds from LHC searches for colour production with b-tagged jets will contribute to push these states above the 1−1.2 TeV scale [48][49][50][51]. At the same time there exist strong bounds from the measurement of B s mixing [52] that either strongly limit the available size of the Yukawa couplings or push one of the new states coupled to b, s to a very large scale. In the case where there is a single pair of particles per vertex, ordered like in Fig. 1(a), the bounds from B s mixing require |Y where we have suppressed subscript indices in the Yukawa couplings, we define x = (m 1 /m 4 ) 2 , y = (m 2 /m 4 ) 2 , z = (m 3 /m 4 ) 2 , and F is a loop function, which equals approximately 1 when m 1 ≈ m 2 ≈ m 3 ≈ O(100 GeV). As shown, e.g., in refs. [39,40,46,47], it follows from Eq. (2) that C µ 9 −0.6 requires Yukawa couplings |Y , a condition that does not guarantee the validity of perturbation theory up to scales much larger than EWSB.
In line with the recent interest in light dark sectors, one can perhaps envision enhancing the contributions of the muon Yukawa coupling by reducing the mass of the colour-singlet particles in the box. However, this is not possible without incurring severe bounds from multiple sources.
One in fact infers from Eq. (3) that F receives a logarithmic enhancement of a few units if x, y, and z are all at the same time significantly smaller than 1. But at least one among m 1 , m 2 , m 3 cannot be much smaller than m 4 . There are two reasons for this. On the one hand, multi-lepton searches at the LHC via Drell-Yan production constrain the particles carrying SU(2) L ×U(1) Y quantum numbers to masses above the 400−600 GeV range [53][54][55]. On the other, if one assumes that the uncoloured particles in the box are lighter than the mass of the muon, so to possibly result invisible at the LHC, they still necessarily contribute at one loop to the anomalous magnetic moment of the muon. The measured 2σ upper bound, δ(g − 2) µ 4 × 10 −9 [56][57][58] implies that the muon Yukawa coupling be This is too small to fit the flavour anomalies unless, as we said above, m 1,2 ∼ O(100 GeV) at least. We conclude that there is arguably no common parameter space for b → s anomalies and (g − 2) µ with perturbative Yukawa couplings and NP masses below O(100 GeV).

Penguin diagrams
It is more promising to look at another type of loop-induced coupling giving rise to the effective operators O 10 is flavouruniversal, and cannot be used to explain the LFUV anomalies. We will therefore be interested in models realising the penguin diagram topology presented in Fig. 1(b), in which a new vector particle (referred to henceforth as a dark photon for simplicity) couples to the muons, and the b and s quarks.
The particle that carries the colour charge in the loop might be a scalar multiplet Φ or a heavy vector-like (VL) fermion Ψ. In the former case the loop is closed by a light (Dirac) fermion χ, whereas in the latter by a light scalar π. In order to avoid charging the b and s quarks under the dark gauge group U(1) D (with gauge coupling g D ), we assume that the dark charge is confined within the loop, i.e., Q χ(π) = −Q Φ(Ψ) . We thus avoid the strong bounds on V from multi-lepton searches at the LHC.
Without much loss of generality we will focus henceforth on the case where the heavy coloured particle is Φ. The case with Ψ does not present very significant differences, barring a few percent suppression of C µ 9 due to swapping the role of the light and heavy mass in the loop functions (see below). We thus have a scalar doublet Φ = (φ u , φ d ) T and light fermion singlets χ i , whose quantum numbers with respect to SU (3) Note that with the above charge assignment the mass matrix of the dark fermions is not a priori diagonal, unless the off-diagonal entries are forbidden by a flavour symmetry, or suppressed by some other mechanism. We will assume for simplicity that this is always the case, without entering in the specifics of such constructions. We confine ourselves to the treatment of left-handed b − s currents, which give rise to the operators O µ 9,10 . Below the EWSB scale the Lagrangian of the hadronic NP sector reads where a sum over repeated indices is implied, and the Yukawa couplings are related by y ij u = y ik d (V † CKM ) kj . We further confine ourselves to the basis where the only nonzero Yukawa couplings of the down-like type are those of the second and third generation: We finally introduce an effective interaction of the dark photon to the muons: The relative size of the g V µ and g A µ couplings is governed by the UV model-building and so is the eventual size, if nonzero at all, of the coupling of the dark photon to neutrinos. This point is of particular relevance when it comes to the constraints on the model from neutrino trident production [59] at CCFR [60] and CHARM-II [61], to which we come back in Secs. 3 and 4. We shall see that, while it is desirable to embed the effective model in a framework with forbidden or strongly suppressed couplings to the neutrinos, 1 we are able to find in our numerical analysis some viable parameter space lying below the bound from CCFR and CHARM-II.
The penguin diagram contributions to the effective operators O µ 9 , O µ 10 are given by [46] where m V is the mass of the dark photon, we have defined and Γ V is the total width of the dark photon, which reads when all light decay channels are kinematically open. Finally, the loop functions are given by Equation (7) allows one to define the effective coupling g V bs to left-handed quarks as 1 For example, one can construct a Type-I 2-Higgs doublet model where the additional Higgs doublet and extra VL heavy leptons all carry U(1) D charge: Θ : (1, 2, 1/2, Q Θ ), E : (1, 1, 1, Q Θ ). Yukawa couplings with the SM left-handed muon doublet L µ , of the type λΘ † L µ E, generate a left-chiral coupling g L µ of V to the muons once U(1) D is broken. The coupling to neutrinos is absent. Typically one gets g L µ ≈ −g D Q Θ λv 2 Θ /2m 2 E . More than one family of VL singlet fermions, and an additional complex scalar singlet charged under U(1) D can be introduced to generate the right-chiral couplings and to push the mass of the scalars with electric charge above the current bound from LEP and LHC searches [62]. where we have introduced the dimensionful couplingg and left explicit the dependence of the coupling on the momentum of V . This corresponds to the effective operator O 6 defined in ref. [35] and also used in ref. [34].
The typical size of g V bs for representative choices of the input parameters is shown in Fig. 2.
Equation (7) shows that the penguin-generated Wilson coefficients depend on q 2 and the size of the dark-photon width. They enter nontrivially in the calculation of the flavour observables and relative constraints when the dark-photon mass lies in the vicinity of the experimental bins. The corresponding expressions should thus be compared directly to the experimental data, as we do in the numerical scan presented in Sec. 4.
We can nonetheless provide a rough estimate of the typical size required for the q 2 -independent couplingg in the limiting cases that the dark photon mass m V becomes much larger or much smaller than the experimental energy (for the LFUV observables we indicatively consider this to be q 2 ≈ 0.04 − 8 GeV 2 ).
If m V 10 GeV one obtains C µ 9 −0.6 roughly for where we have set q 2 at the mean momentum transfer for the 1 − 6 GeV 2 bin of the R K ( * ) observables.
The two couplingsg and g V µ are independently constrained by two powerful probes. On the one hand, the measured 2σ bound on R BB from B s mixing [52] has a direct impact ong when the vector V is exchanged in the s-channel: 2 where the characteristic experimental scale coincides with the B s mass, q 2 = M 2 Bs . On the other hand, the already mentioned 2σ upper bound on the anomalous magnetic moment of the muon constrains g V µ to The current bounds thus result on a very narrow window of availability for explaining the flavour anomalies with penguin diagrams and a heavy dark photon, m V 10 GeV. At the opposite side of the spectrum, m V < 200 MeV, the dark photon is much lighter than the experimental energy scale. Equation (7) shows that the Wilson coefficients become independent of the dark photon mass and of q 2 . C µ 9 −0.6 requires approximately The most severe constraint in this mass range comes again from the 2σ upper bound on the anomalous magnetic moment of the muon, which requires where the upper value refers essentially only to the region m V 100 MeV, and the lower value to all other masses below that. In light of Eqs. (15) and (16) one needs |g| 10 −6 GeV −2 to fit the flavour anomalies. As one can see in Fig. 2, effective couplings of this size are not easy to obtain in the penguin setup.
It is well known (cf., e.g., ref. [33]) that one can increase the size of g V µ while respecting the 2σ upper bound on (g−2) µ by introducing the axial-vector coupling g A µ and thus creating a negative contribution to (g − 2) µ that has to be fine-tuned. This, however, also induces an extremely strong contribution to the B s → µµ decay rate, and will prove ultimately impracticable, as will appear clear in the next sections.
Altogether, this discussion leads us to conclude that the most natural solutions are likely to be situated inside the window m V ≈ 0.2 − 10 GeV. The remainder of this paper is thus dedicated mostly to this mass range. Note that when the dark photon mass approaches M B , additional resonant enhancement can be obtained to open up the parameter space, albeit at the cost of additional limits from B → K processes, which we will describe in detail in the next section.
3 Constraints on the model 3.1 Flavour constraints B s -mixing Strong limits on the Yukawa couplings to b and s quarks arise from box-diagram contributions to B s -mixing. We recall that for exclusively left-handed couplings, see Eq. (5), the only relevant operator is O 1 [40,46]. The corresponding dimensionful Wilson coefficient is defined as where the sum runs over all possible dark fermions, x i(j) = m 2 χ i(j) /m 2 Φ , and the loop function reads In the small x, y limit, the loop function can be approximated as F (x, y) 1, so that C 1 ∝ 1/m 2 Φ and the limit essentially saturates. This has the unexpected consequence that, in the presence of several dark fermion states, one can readily get a large suppression of the B s -mixing contribution in the limit where ij y i * s y i b y j * s y j b 1, which can be easily obtained with, e.g., two light states χ 1 , χ 2 and y 1 * s y 1 b ≈ −y 2 * s y 2 b . Note that while C 1 receives a strong suppression from the addition of approximately equal and opposite-sign couplings, there is no equivalent suppression for C µ 9 , as in the limit of small x 1 , x 2 the effective couplingg in Eq. (11) becomes proportional to ln(x 1 /x 2 ).
Limits from B s decay When the axial-vector coupling to the muon g A µ is present, it induces a contribution to the operator O µ 10 but also to the axial current operator O µ P , which in turn induces a contribution to the direct decay B s → µ + µ − [63][64][65][66].
By adopting the same convention for the scalar operator as in ref. [46] we obtain The typical bounds on C P are significantly more stringent than those on C 10 . They are likely to have a strong impact on our results so that we include them directly in the full numerical scan.
Given the presence of the heavy scalar doublet Φ and dark fermion χ, one can construct the tree-level decay process B → Kχχ based on the b quark 3-body decay b → χΦ * → sχχ. In the limit where m χ , M K M B , we have the simple expression: where f 2 + ≈ 0.3 is the average value of the form factor over the range of integration of the differential decay rate. This typically leads to the constraint |y * s y b | 10 −2 (m Φ /TeV) 2 on the Yukawa couplings of any new fermion with mass 2m χ < M B − M K . As Eq. (11) shows, the upper bound on the Yukawa couplings strongly limits the available range of g V bs , even when the gauge coupling g D is large. This means that to fit the flavour observables one has to resort to a large g V µ value, which, as we shall see in Sec. 4, is severely constrained by Z-lineshape bounds (in addition to requiring a fine tuning of the axial-vector contribution to avoid exceeding δ(g − 2) µ ). To avoid these problems the minimal particle content will have to include at least on dark fermion with mass m χ > (M B − M K )/2.
On the other hand, there is a case to be made for the presence of additional light states with mass below the (M B − M K )/2 threshold and Yukawa couplings to Φ and the b and s quarks that are small enough to avoid the tree-level B → Kχχ bound.
• Dark matter : for a fermion χ of mass above the (M B − M K )/2 threshold the direct swave annihilation channel, χχ → V → µµ, is strongly constrained by CMB bounds [72]. Introducing one additional lighter state with the same quantum numbers provides instead a viable candidate for forbidden dark matter (we come back to this point in more detail at the end of Sec. 4) • The presence of additional light states directly affects the total width of the dark photon V , opening up additional parameter space for a solution to the flavour anomalies.
In cases where at least one very light state χ 1 appears, two qualitatively different regimes of applicability should be considered for the invisible B → K transition. The first is on-shell decay B → KV , occurring for m V < 2M µ < M B −M K < 2m χ 1 , in which V escapes undetected. It is typically suppressed in the low m V limit due to the momentum dependence of the effective coupling g V bs defined in Eq. (11). We get where we have assumed m V , M K M B , f + (0) 0.3 is a form factor (the full expression, used in the numerical analysis of Sec. 4, can be found in ref. [73]), and λ(x, y, z) is the standard Exchange of a virtual V may dominate the invisible decay width for small m V and large dark coupling g D . The full width reads where Γ KV (s) and Γ V →χ 1 χ 1 are obtained by replacing m V by √ s in Eq. (22) and in the corresponding V decay width to fermions, and Γ V is the total width of V , cf. Eq. (9). Note that with more than one dark fermion in the spectrum one ought to sum over all individual off-shell contributions.
Altogether, the combination of both real and virtual contribution to the B → Kχ 1 χ 1 decay implies a complex kinematic shape in terms of missing energy, which may differ significantly from the SM-like B → Kνν decay. This has the direct consequence that the experimental results of Belle and BaBar, which are optimised for the neutrino process, should not be directly applied to our scenario. We therefore perform a conservative recasting of these analyses, described in detail in Appendix A.
Finally, in the mass regime m V > 2M µ , the light vector state can directly decay into a muon pair. This opens up the resonant channel B → K * V , V → µµ. The typical decay width into K * is given by where an explicit expression for F 1 (x K * , x V ) can be found in Appendix B of ref. [35]. Note that in the limit where m V M B , F 1 (x, y) simplifies to F 1 (x, y) 0.1 √ x(y +1.2x), so that Γ K * V is typically suppressed compared to Γ KV . Furthermore, the branching ratio to muons is inversely proportional to the total width Γ V and can thus be strongly suppressed if the coupling of V to the dark fermions χ i is large.
Note that the limits from LHCb [71] on this process focused on a narrow, or even longlived, new resonance, with limits based on invariant mass bins of a few MeV. This hypothesis is especially problematic when the mediator is around the GeV range, since a large dark gauge coupling implies a very large width for V . We then simply model the resonance via a Breit-Wigner distribution and compare bin-per-bin with the limit of ref. [71], retaining the strongest bin as the main limit. 3 The overall impact of the limits from B → K ( * ) transitions on the size of the effective coupling g V bs is summarised in Fig. 3 as a function of the vector mass m V , for three representative choices of the pair (m χ , g D ). For m χ = 5 GeV, the decay B → Kχχ is kinematically forbidden, so that invisible B decays can only proceed via the on-shell process B → KV . This results in a weak bound in the small m V regime (blue line). Note that when m V > 2M µ , the constraints from resonant searches using B → K * µµ strongly exclude this setup since V → µµ is in this case the only accessible decay channel.
In the presence of dark fermions with a small mass, 2m χ < M B − M K , the decay channel V → χχ opens up and the relative strength of different bounds depends on the size of g D . In the case of large g D (orange line), if 2m χ < m V < M B − M K both the on-shell and off-shell invisible decays contribute to the width. For light vector mediator the off-shell decay takes over due to the q 2 dependence of the coupling and the bound on g V bs saturates. If, on the other hand, g D is smaller (green line) and m V > 2M µ , the constraints from resonant searches using  B → K * µµ typically overcome the invisible decay limits, suppressing g V bs by a factor of the order of 1.
Finally, note that for large m V the limit arises from invisible searches from the off-shell decay B → Kχχ when it is kinematically allowed.

Precision physics constraints
Muon anomalous magnetic moment The couplings of V to muons can be constrained by the measurement of the anomalous magnetic moment of the muon. A contribution to (g−2) µ is in this case given by [74,75] where As was mentioned in Sec. 2, given the limited range achievable in penguin constructions for g V bs , when only vector-like couplings to the muons are present it becomes difficult to find a g V µ value large enough to allow for a reasonable agreement with the flavour anomalies and at the same time not too large a deviation from the measured value of (g − 2) µ . A certain level of cancellation with the contribution from the axial-vector coupling must take place in most situations [33]. For a GeV-scale vector mediator, this occurs for g A µ ≈ −0.44 g V µ . Note, however, that including the axial-vector contribution triggers the strong bounds from B s → µµ discussed above.
Z physics and intensity frontier limits The coupling of the Z boson to the muon is modified at the one-loop level [76,77] within UV constructions of the type as in Footnote 1. However, due to the smallness of the Z-boson coupling to charged leptons in the SM, the limit is typically subdominant with respect to the (g − 2) µ bound.
A powerful method for discerning light resonances through precision measurements of Drell-Yan dimuon production was proposed in ref. [78]. For m V = 1 − 5 GeV an upper bound can be derived Finally, the Belle-II Collaboration recently provided a bound on the final state radiation process e + e − → µ + µV , V → invisible, based on 0.28 fb −1 of data from the 2018 run [79], which applies directly to our model. While the current limit can hardly compete with the Drell-Yan bound, the 2019 run has stored ∼ 10 fb −1 and moreover a few ab −1 should be obtained in 2020, so that future data will become rapidly relevant.
Notice that more generically, for dark photon mass above ∼ 10 GeV, the phenomenology of the light fermions in intensity frontier experiments can be obtained by integrating it out and considering the fermion portal four-fermion operators [80]. Similarly, the limit from the tree-level b → χΦ * → sχχ decay can also be obtained by integrating out the heavy scalar Φ and using the existing bounds on the fermion portal operatorbγ µ sχγ µ χ. While we cover directly the relevant limits in this section, the latter approach could be particularly fruitful to study and constrain the possible couplings between new light fermions and the other SM generations.
Neutrino trident production If the gauge boson features a coupling g ν to muon neutrinos (cf. Sec. 2.2), one expects a strong enhancement in the neutrino trident production from scattering on atomic nuclei, N : N ν → νN µ + µ − [59].
The cross section for this process has been measured by the CCFR [60] and CHARM-II [61] collaborations to be in agreement with the SM prediction. In the range m V > 1 GeV it results in the generic bound which roughly saturates for smaller mass to about g ν 0.001.
When the bound appears in the plots of Sec. 4 it is obtained under the assumption g ν ≡ g V µ (we repeat that whether or not the bound is relevant depends on the UV completion).
Kinetic mixing In presence of states charged both under the U(1) Y and U(1) D symmetry groups, kinetic mixing between the photon and the vector V will be generated at the loop level. The corresponding 1-loop contributions from fermions and scalars are given by where Y f,s is the fields' hypercharge, Q f,s is the dark charge, and the coefficients N 3 f,s indicate the dimension of the SU(3) c representation. The fields that contribute to the kinetic mixing are Φ, µ L , and µ R which, when g V,A µ g D , results in Such a kinetic mixing is at the limit of exclusion given the current intensity frontier searches (see, e.g., ref. [81] for a recent review), especially when invisible decay channels are not available for the vector mediator. However, the precise value of the kinetic mixing is strongly dependent on the UV physics, and additional VL fields can modify the prediction of Eq. (29) although not by many orders of magnitude.
LHC constraints on t → cµµ We work in this paper under the assumption that the only nonzero Yukawa couplings of the down-like type are y i s , y i b . However, as Eq. (5) shows, the corresponding Yukawa couplings of the up-like type, y i c , y i t , do not receive CKM suppression. They can thus generate non-negligible contributions to processes involving t → cµµ transitions.
Effective operator analyses of LHC bounds from rare top decays [82,83], derived originally for the very-high mass regime, m V ∼ O(TeV), impose a fairly weak bound on the coupling product when m V M t : A rough comparison with Eq. (12) shows that this is not likely to be constraining for our scenarios. However, given that Eq. (7) presents a nontrivial q 2 dependence, for the light mass range investigated in this paper one should rather perform a detailed recast of the experimental searches. This tasks exceeds the purpose of the present paper in view of the fact that, as we will show in Sec. 4, the flavour and intensity frontier experiments discussed above provide already a set of powerful and often inescapable constraints on the dark sector.

Fitting procedure and results
Fitting procedure We perform a multidimensional fit of the following free parameters: with only one light dark fermion contributing with Yukawa coupling y s and y b ). r is employed here as a proxy for the effective couplingg once the mass parameters are fixed. We choose m χ = 2.5 GeV and m Φ = 1 TeV; under these assumptions r relates tog asg = 1.1 × 10 −8 GeV −2 r. Since the fitted flavour observables depend ong rather than the couplings composing it, our results can be extended straightforwardly to the case with more light fermions.  The prior ranges of the fitted parameters are presented in Table 1. Separate fits are performed whether m V lies above or below the relevant bins for B anomalies. This is required by the fact that, in order to reproduce the desired sign for C µ 9 , the product r · g V µ needs to have a different sign in these two regions. We observe that m V = 2 GeV works as a satisfactory threshold to separate the two regimes.
The first step of our numerical analysis is to perform a fit of the free parameters of the model to the available experimental data reporting anomalies in B-meson decays, namely the LFUV ratios R K , R K * , the angular observables in the B → K * µ + µ − decay and the branching ratio of B s → µµ. 4 In the high-mass regime we include in the likelihood function the experimental 2σ upper bound on the anomalous magnetic moment of the muon, δ(g − 2) µ 4 × 10 −9 .
To carry out the fit, we employ the HEPfit package [84], performing a Markov Chain Monte Carlo (MCMC) analysis by means of the Bayesian Analysis Toolkit (BAT) [85]. A set of O(500K) points is generated for the two scenarios described in Table 1, and for each case the subset of points reproducing the B anomalies, BR(B s → µµ), and the upper (g − 2) µ bound (high-mass only) at the 2σ level, is stored. These initial subsets of points are then subjected to additional constraints coming from B s mixing, Drell-Yan production and, when applicable, B → K + inv. searches.

Results
We start with addressing solutions in the range m V = 2 − 15 GeV, in cases where the dark photon features a non-negligible invisible width, which could stem from the presence of an additional light fermion χ 1 in the spectrum which do not participate significantly to the loop diagram. The results of the scan in the (m V , g V µ ) plane are presented in Fig. 4. The yellow points are obtained in the scanning procedure described above and correspond to the models in which the B anomalies, BR(B s → µµ), and δ(g − 2) µ are fitted at the 2σ level. The green 4 Due to the explicit q 2 dependence of C µ 9(10) it is not possible to match the Wilson coefficients to the modelindependent bounds obtained by the global fits. Instead, C µ 9(10) have to be fitted directly to the experimental data. On the other hand, if the dark photon mass is in the MeV range, it is possible to directly match the Wilson coefficients to the model-independent bounds obtained by the global fits.  2) µ are fitted at the 2σ level. Superimposed points are allowed after the limits from B → K + inv. are applied, as reported by CLEO [86] (purple) or using our recast procedure (green), cf. Appendix A. Grey shaded region is excluded by precision measurements of Drell-Yan at the LHC [78]. Purple dotted line indicates the exclusion bound from the neutrino trident production [59]. Dotted black line indicates the upper bound from Belle II [79]. In light blue we mark the region in which the (g − 2) µ constraint is satisfied at 2σ with g A µ = 0.
points are those that remain allowed after the limits from B → K + inv. transitions are applied, following the recasting procedure outlined in Appendix A. The on-shell process B → KV , V → χ 1 χ 1 typically proceeds unsuppressed for a GeV-scale mediator, so that all the yellow points below the M B − M K threshold are excluded. On the other hand, the limits ong are dramatically weakened above the M B − M K threshold (see Fig. 3), so that in this regime one can easily fit simultaneously both (g − 2) µ and the flavor anomalies. Incidentally, we exclude m V ≈ 2.5 GeV, corresponding to the solution found in ref. [33], due the q 2 dependence of g V bs , which induces an (m V /GeV) 2 enhancement in Eq. (22) with respect to the effective coupling of ref. [33].
An additional bound on the parameter space of the model is derived from the Z lineshape from Drell-Yan at the LHC, which strongly affects the maximal allowed value of g V µ [78]. The corresponding exclusion region is depicted in Fig. 4 in dark grey. The shading is obtained under the assumption g A µ = −0.44 g V µ , a relation that induces destructive interference in the calculation of (g − 2) µ (cf. Eqs. (25) and (26)). It is well known [33] that the above relation between the vector and axial-vector coupling requires some level of fine tuning. The grey dashed lines in Fig. 4 trace the value of g V µ that corresponds to the indicated level of fine tuning in g A µ ≈ −0.44 g V µ required to avoid exceeding the 2σ upper bound from the measurement of δ(g − 2) µ .
Note, however, that the tuning of the vector and axial-vector muon couplings is a priori not needed for the dark photon in the considered mass regime, as confirmed by the large number of green points within the blue shaded band, corresponding to the region satisfying the (g − 2) µ constraint at 2σ with g A µ = 0. This is an attractive feature of our model, in which we can obtain relative low values of g V µ and subsequently avoid the Drell-Yan limit. In Fig. 4 the limits from the neutrino trident production derived in ref. [76] for g A µ = 0 are shown as a dashed purple line. The constraint applies only if the mediator couples directly to neutrinos, cf. discussion in Sec. 2.2. It should be stressed that, even when the neutrino trident bound applies, solutions that escape the experimental limit exist, with m V = 5 − 7 GeV and r at the upper end of the scanned range.
Finally, it is instructive to compare the results obtained within our UV-complete setup to those derived in a simplified Z model with an effective coupling to the b − s current, g bs . A solution to the b − s anomalies (without the (g − 2) µ constraint) was found in ref. [18] with m Z = 10 GeV, g V µ = −g A µ ≈ 10 −2 , and g bs ≈ 5 × 10 −6 . Given the bound from the BR(B → K + inv.), which limits the value of g D Q Φ , and the requirement of perturbativity of the couplings, we never reach these high values of g bs ≡ g V bs in our scenario, see Fig. 2.
We conclude the discussion of Fig. 4 by pointing out that the picture does not receive substantial modifications if the light fermion χ 1 is not introduced in the theory. If there remains only one dark fermion of mass m χ = 2.5 GeV the B → K + inv. bound does not apply. However, all points with mass m V < M B − M K become in that case subject to the strong resonant B → K * µµ limit, which cuts drastically the parameter space and induces solutions not dissimilar to the area delimited in green in Fig. 4.
Note that the LFUV and angular observables in the fit depend only indirectly on the NP Yukawa couplings, via the effective couplingg. As discussed in Sec. 3.1, the former are thus constrained only by the tree-level contribution to B → Kχχ, which is kinematically forbidden for m χ ≥ (M B −M K )/2, and B s mixing, which dominantly proceeds via box diagrams involving dark fermions and heavy coloured scalars. In the presence of several light fermions, a possible way of suppressing B s mixing in the limit m χ i /m Φ 1 is obtained when ij y i * s y i b y j * s y j b 1, even if one of the χ i is relatively heavier than the others.
The dependence of B s -mixing bounds on the amount of fine tuning among the NP Yukawa couplings and particle masses is illustrated in Fig. 5. In the (γ V , g V bs ) plane, green points correspond to the solutions also marked in green in Fig. 4, which were consistent with the flavour anomalies and also escaped B → K + inv. limits (assuming that the invisible decay of the mediator is ensured by the presence of a very light fermion with small Yukawa couplings, as discussed above). Yellow points trace instead the solutions respecting the resonant B → K * µµ bound when there is no extra light fermion in the spectrum. The parameter space is very similar, showing perhaps a slight preference for small width when the bound B → K + inv. applies. A third dark fermion with large Yukawa couplings is introduced to cancel the strong The results of the scan in the (m V , g V µ ) plane for the low-mass range, m V = 0.6 − 2 GeV are presented in Fig. 6. This region of the parameter space is obtained by performing a sign switch in the product r · g V µ , which allows one to fit R K ( * ) correctly by means of destructive interference with the SM value of C µ 9 below the experimental bin. The colour code is the same as in Fig. 4. Note that the mediator must have a sizeable invisible width in this region to avoid stringent constraints from a visible dimuon resonance in the B → K * spectrum, as discussed in Sec. 3.1.
A few takeaways emerge from the scan in the low-mass region. The first is that all solutions with mass m V 1.4 GeV are cut out by the BR(B s → µµ) constraint, which is directly implemented in the likelihood function. As a consequence, the surviving points are all characterised by fairly large values of R K ( * ) , which lie 2 − 3σ away from the central value measured at LHCb Recast % % Drell-Ya n ( ) ± Figure 6: Results of the scan in the low m V region. The colour code is the same as in Fig. 4. and closer to the SM expectation. It is for the same reason that the plot also appears much sparser than in the high-mass case: there are few model points within 2σ of the measured values of LFUV observables and BR(B s → µµ). We thus identify a mild tension in this part of the parameter space. We then notice that the models surviving the bound from B → K + inv. searches (green points) require large coupling g V µ to the muon (g is directly constrained by the invisible search) and thus a large level of fine tuning in the corresponding g A µ value necessary to cancel δ(g − 2) µ . Overall, the whole region is under siege from a combination of complementary bounds but at present it is not entirely excluded.
Let us finish this section by mentioning the case of a very light mediator, with a dark photon with mass in the MeV range and GeV-scale dark fermions. An interesting property of this regime is that, in the limit where 2m χ i > M B − M K , the mediator V is essentially long-lived since is does not have any available tree-level decay channel. The invisible B decay is then driven exclusively by the B → KV process, which is strongly suppressed at low m V . Furthermore, the dependence of the Wilson coefficients on q 2 /(q 2 −m 2 V ) converges to a constant when m V q and it closely resembles the standard electromagnetic penguin contribution to the flavour anomalies. However, given that the limit from B s → µµ decay forbids such a light vector mediator to have a significant axial-vector coupling to the muon, the upper limit on δ(g − 2) µ leads to a stringent bound on the vector coupling, g V µ 7 × 10 −4 . Concretely, in order to obtain C µ 9 ≈ −0.7 a couplingg 10 −6 GeV −2 is required, which can only by achieved while satisfying B → K + inv. limits if m V 5 MeV, cf. Fig. 3. A vector mediator this light is already excluded by the standard searches for long-lived dark photon. We conclude that no such solution is available in penguin-generated scenarios.
Dark matter As an interesting aside, the lightest dark fermion can provide a good example of forbidden dark matter candidate when its mass is below the muon mass. The dominant annihilation channel for such a dark matter candidate would be χχ → V * → µμ. Such a process leads to a typical relic density of where x f ≈ 20 for the relevant masses, g D = √ 4π and ∆ = 1 − m χ /M µ . When m χ drops below the muon mass threshold the relic density is exponentially enhanced since the annihilation process can only occur due to the thermal velocity of the dark matter particle in the early universe [87,88]. This ensures that the thermal target is matched for one coupling-dependent mass below M µ , typically around ∼ M µ /2. Furthermore, all other annihilation processes are exponentially suppressed when the universe temperature decreases, ensuring that the CMB limits on late-time annihilating sub-GeV dark matter are automatically escaped.

Conclusions
We have presented in this work a solution for the b → s flavour anomalies based on the presence of a split dark sector with a light vector mediator as well as new light Dirac fermions which may constitute all or part of the dark matter. The interaction with the b and s quarks is generated at the loop level via the addition of a coloured scalar particle, resembling a supersymmetric squark. We analysed numerically and analytically the resulting low-energy effective theory, which in particular possesses a q 2 -dependent interaction of the vector mediator with b and s quarks. Varying the mass of the vector mediator from the MeV scale to the tens of GeV, we find two scenarios satisfying all experimental constraints while providing a good fit to the anomalies. In particular, the region with a GeV-scale mediator above the B − K mass threshold appears particularly promising, requiring little to no tuning in the low-energy effective parameters, and to the best of our knowledge it has not been considered previously.
Since our model is partially embedded in a UV completion, we additionally pointed out several constraints that can challenge its viability. We highlighted the constraints from the B → K + inv. decay rate and B s mixing and, in the latter case, provided an example of a mechanism to escape it. We did not make any assumption in this paper on the nature of the dark Dirac-fermion interactions with the neutrinos. Indeed, since the former are complete SM singlets, it would be very interesting to investigate whether or not they could behave as righthanded neutrinos (for instance via the coupling to a dark charged new Higgs doublet), and in that case investigate their relationship with the strong neutrino trident limits.
While the experimental constraints on models addressing the flavour anomalies with light mediators are already quite stringent, we have identified several observables that can easily exclude these scenarios entirely or provide smoking-gun proof of their detection. Chief among those are the limits from B → K and B → K * transitions. While the former play a critical role via the B → K + inv. bounds, experimental searches are typically optimised for the SM process B → Kνν. Including an analysis based on a light, and potentially broad, invisible resonance B → KV could likely strengthen significantly the existing limits, especially when the mediator is light. Similarly, the latest search for B → K * µµ has focused on a very narrow resonance, and should be properly recast for the case of a large and invisible width. Finally, it is important to note that limits on a light dark photon are due to improve in the next few years, and will further constrain the case of a mediator at and below the GeV scale.

A Appendix: Invisible decay limits
We present in this appendix a more detailed treatment of the recasting procedure performed to extract conservative B → K + inv. limits on our model. While the CLEO Collaboration searched explicitly for an on-shell light particle mediating the B → K decay [86], their 2σ limit BR(B → K + inv.) < 4.9 × 10 −5 is relatively weak compared to the ones from B factories. In the following we will base our limit both on the BaBar result [67] -which provides a differential branching ratio limit in bins of q 2 -as well as on an older analysis from the same collaboration [68], which also had some differential limits, albeit on a much larger range for q 2 . Note that the current bounds from the Belle Collaboration [69,70] are typically of the same order as for BaBar. They are however strongly optimised for the SM-like B → Kνν signal and only present their bounds in the total integrated branching ratio. We therefore concentrate on the two BaBar analyses.
First, using the BaBar hadronic-tagging analysis [67], we calculate the branching ratio BR(B → Kχ 1 χ 1 ) in s B = q 2 /M 2 B bins, where for low m V most of our NP signal is concentrated in the lowest bin, s B < 0.1. We then compare it with the 2σ limits from Fig. 6a of ref. [67]. While this approach leads to a strong bound when the real process B → KV , V → χ 1 χ 1 dominates, these limits can be significantly weakened when the virtual process dominates, since the branching ratio accounts for a broader spread in s B bins.
We therefore also include partially integrated limits from the BaBar semileptonic-tagging analysis [68], which combined the world-leading limit on B + → K + νν with detailed information about the signal efficiencies as function of the momentum of the K + (and hence on the missing energy). We select the q 2 ranges [3.4 2 , 4 2 ] GeV 2 and [0, 2.4 2 ] GeV 2 (corresponding to p K in the range [1, 1.5] GeV and 2 GeV, respectively), where the Boosted Decision Tree (BDT) efficiencies presented in Fig. 3 of ref. [68] are larger than ∼ 0.3. This ensures that the signal efficiencies for our NP kinematics are of the same order of magnitude or higher than the ones for the SM signal. We then compare both regions with the low-q 2 and high-q 2 95% C.L. limits, BR(B → K + inv.) < 1.1 × 10 −5 and BR(B → K + inv.) < 4.6 × 10 −5 , respectively.