Flavour anomalies and the muon $g-2$ from feebly interacting particles

We perform a phenomenological analysis of simplified models of light, feebly interacting particles~(FIPs)that can provide a combined explanation of the anomalies in $b\to s l^+ l ^-$ transitions at LHCb and the anomalous magnetic moment of the muon. Different scenarios are categorised according to the explicit momentum dependence of the FIP coupling to the $b-s$ and $\mu-\mu$ vector currents and they are subject to several constraints from flavour and precision physics. We show that viable combined solutions to the muon $g-2$ and flavour anomalies exist with the exchange of a vector FIP with mass larger than $4 \,\textrm{GeV}$. Interestingly, the LHC has the potential to probe this region of the parameter space by increasing the precision of the $Z\to 4\mu$ cross-section measurement. Conversely, we find that solutions based on the exchange of a lighter vector, in the $m_V<1\,\textrm{GeV}$ range, are essentially excluded by a combination of $B\to K +\textrm{invisible}$ and $W$-decay precision bounds.


Introduction
Feebly interacting particles (FIPs) represent a very large and well-motivated category of new physics (NP) scenarios. Loosely defined as new light particles with mass below the electroweak symmetry-breaking (EWSB) scale and with a feeble interaction with Standard Model (SM) fields, FIPs encompass NP particles as diverse as the dark photon and axionlike particles. A particularly exciting possibility [1][2][3][4][5][6][7][8] is that one or more FIPs may be responsible for the anomalies in lepton-flavour universality violating (LFUV) observables recently confirmed in new data from LHCb and for the discrepancy between the measured value of the anomalous magnetic moment of the muon and its SM expectation.
Equally intriguing is the recent measurement of the anomalous magnetic moment of the muon by the E989 experiment at Fermilab [37]. The experimental collaboration reports a 3.3 σ deviation from the value expected in the SM [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57]. When the new measurement is statistically combined with the previous experimental determination, obtained a couple of decades ago at Brookhaven [58], one gets [37] δ(g − 2) µ = (2.51 ± 0.59) × 10 −9 , (1.1) which yields a deviation from the SM data-driven prediction at the 4.2 σ level. If the LFUV anomalies and δ(g − 2) µ are due to the interactions of one or more FIPs in the MeV to GeV range, one expects the accompanying presence of new flavoured physics around the TeV scale (see, e.g., Ref. [8]). On the one hand, higher-dimensional FIP interactions, such as the ones characteristic of axion-like particle models, are simply effective field theories (EFTs) requiring an ultraviolet (UV) completion. On the other hand, even for renormalisable FIP interactions the presence of NP above the EWSB scale is often required, for instance to evade the strong bounds from neutrino trident production, which otherwise drastically constrains a solution to the (g − 2) µ anomaly based on light states [59].
While a model-independent description of heavy physics in flavour-violating meson decay via the Weak Effective Theory (WET) has proved invaluable in identifying the specific operators associated with the anomalies emerging in a large set of observables, the WET fails to account properly for the possibility of one or more extra degrees of freedom that are light but hidden, as it is not by construction equipped to take into account the effects of momentum-dependent couplings, or the presence of possible resonances in the experimental energy bins. We thus consider in this work EFT operators more suited to the study of b → s and µ-related physics under the assumption that the unspecified heavy NP is accompanied by a light FIP [60,61] (see also dark-matter motivated constructions, for instance Refs. [62][63][64][65][66]).
To this end, we construct a set of operators that parametrise the interactions of the FIP with Lorentz-invariant vector bilinears of the b− s and µ −µ current, in agreement with the measurement of the LFUV and (g − 2) µ anomalies. They are characterised by increasing powers of the FIP momentum transferred in the low-energy process. We then perform a comprehensive analysis of the constraints that can be applied on the Wilson coefficients of these new operators. Constraints arise from several sources: direct measurements of the branching ratios for B → K + invisible and B → K * µ + µ − resonant decays, the measurements of BR(B s → µ + µ − ) and B s -mixing, and several precision measurements associated with the W -and Z-boson decay widths.
The purpose of our analysis is twofold. Firstly, we expand significantly on our previous study of a light "dark" U(1) D gauge boson solution to the b − s anomalies with a q 2dependent coupling [8]. On the one hand, we add strong new constraints that were not considered previously in the literature (we compute numerically the best currently available bounds on GeV-scale muonphilic particles from the recent ATLAS analysis of Z → 4µ [67] and the limits based on the violation of flavour universality in W -boson decay [68]). On the other, we perform a global scan of two new simplified models with different momentum dependence of the couplings with respect to Ref. [8]. One of these models, incidentally, overlaps with the case considered in Ref. [2], for which we exclude the region of parameter space highlighted there and point to a slightly heavier preferred solution for the flavour and (g − 2) µ anomalies, with no fine tuning of the vector and axial-vector Z couplings to the muon. The second and main goal of this work is to provide a concise and broad model-independent picture of the current phenomenological status of a light-FIP solution to the muon anomalies. In the process we will establish a "dictionary" allowing the reader to map easily the Wilson coefficients of a FIP-friendly EFT to the WET Wilson coefficients used in the publicly available numerical packages.
The paper is organised as follows. In Sec. 2 we introduce the operators better suited to describing FIP interactions in b → s and LFUV processes. In Sec. 3 we provide the correspondence with the WET basis usually employed in analyses of rare meson decays. The list of applied constraints is introduced in Sec. 4. The result of a global scan in the new operator basis are presented and discussed in Sec. 5, and we conclude in Sec. 6. In Appendix A we show the details required to perform a full calculation of processes involving the B meson in the basis introduced in Sec. 2.

EFT plus light degrees of freedom
We introduce in this section Dark EFT (DEFT) operators of the Lagrangian connecting a generic FIP to the (axial-)vector SM currents relevant for the LFUV and (g−2) µ anomalies.
Spin-1 FIP We parametrise the exchange of a light vector V out of the left-handed b − s current, in agreement with the results of the global fits, in terms of the following operators [5]: where we have used the field strength tensor V ρσ = ∂ ρ V σ −∂ σ V ρ . 1 We complete Eqs. (2.1), (2.2) with the corresponding interactions with the muon current, Since we do not specify the origin of these interactions, Eqs. (2.1)-(2.4) will serve in describing the B-physics for all light vector states. These include, in particular, the GeV-scale top-philic particle in Ref. [69] and the renormalisable model introduced in Ref. [8]. 1 We leave for future work the dimension 5 case, e.g., Q bsV 5 = (sσρσPRb) V ρσ . Since it couples the FIP to the tensor SM current it leads to a starkly different phenomenology.
Spin-0 FIP: the pNGB We consider next the case of a (pseudo-)scalar FIP, a, possibly a pseudo Nambu-Goldstone boson (pNGB). The standard derivative interaction term with the quark current reads and are the corresponding couplings to the muon current. These interactions are strongly constrained by flavour physics and, more importantly for our purposes, do not lead to the vector four-fermion operators O µ( ) 9 and O µ( ) 10 of the WET, which induce the solution to the LFUV anomalies preferred in the global scans. We therefore leave the study of these operators for future work.
FIPs as dark matter: spin 0, 1/2 When the light state is a stable new particle, which can play the role of dark matter, one needs to include operators with two FIP insertions at dimension 6. Using the dark current, we can then define and and equivalent operators can be defined in the scalar case, with χ → S. These effective interactions can be used to parametrise, for example, UV completions to the WET operators O µ( ) 9 , O µ( ) 10 based on dark-matter induced box loops, in the case where the dark matter is relatively light. Note, however, that in such a minimal scenario the dark matter mass is expected to be large enough to evade strong limits from the direct decay B → K + invisible [8].

Relation with the WET
In order to quantify the effects of DEFT operators on the flavour observables, we first create a dictionary between the DEFT and the WET. The procedure is equivalent to considering an explicit momentum dependence of the WET Wilson coefficients, of the type explored, e.g., in Refs. [1,3,6,8]. The specific Lorentz structure of the DEFT operators considered in this work allows us to draw the correspondence with the WET straightforwardly at the level of the spinor fields. We show in Appendix A that the results of this section agree with a full calculation of B-meson decays.

Tree level
Spin-0 FIP The critical ingredient in linking both theories is that we are interested in physical processes in which the muons are on-shell. Given two on-shell external muons of four-momenta k 1 , k 2 , emerging from the exchange of a FIP of momentum q = k 1 + k 2 , the equations of motions yield where M µ is the muon mass. A direct consequence of the above equations is that the derivative scalar current introduced in Eqs. and O µ( ) 10 . By repeating the argument for "on-shell" initial quarks with q = p b − p s , one obtains where M s(b) is the strange (bottom) quark mass. The above equality is useful for deriving the correspondence between DEFT and WET in the remainder of this section.
Spin-1 FIP We write the vector FIP propagator in the form where m V and Γ V indicate, respectively, the mass and total decay width of the light NP vector V ρ . By considering the b → sµµ process via a vector FIP exchange, we can derive the relation between the two effective theories. In the case, for example, of an interaction involving operators Q bsV The above amplitude is expressed in terms of the UV cutoff, Λ, and the DEFT dimensionless coefficients C bsV 6 , C µµV

4
. This leads in turn to replacing the weak effective Hamiltonian for b → sµµ with where we have used the relation (3.3). Following a similar procedure, we can connect to the WET all the operators introduced in Sec. 2a. The WET coefficients C µ 9 , C µ 10 , C µ P , C µ P , can then be expressed in terms of DEFT coefficients. One gets, An explicit UV example We conclude this subsection by matching the DEFT operators to an example of renormalisable UV-complete model in the case of a vector FIP. We focus on the flavour-violating coupling to b and s. Let us consider the SM extended by a "dark" U(1) D gauge group, whose gauge boson is a light V . Let us then introduce a scalar field φ and a multiplet χ of vector-like fermions of mass m χ whose charges with respect to SU(3) c ×SU(2) L ×U(1) Y ×U(1) D are given by φ : (1, 1, 0, Q φ ) χ : (3, 2, 1/6, Q χ ) .
Depending on whether or not the field φ acquires a vacuum expectation value (vev) different DEFT operators of the b − s current are generated at the leading order from the coupling in Eq. (3.12). If φ develops a vev, v φ , the effective coupling of V to b − s, which is generated via the mixing of the SM and NP quarks, does not carry q 2 dependence and thus gives rise to the operator Q bsV 4 . One gets where g D is the dark gauge coupling. The phenomenology of the light gauge boson V in relation to R K ( * ) was analysed, e.g., in Ref. [2]. If, on the other hand, the scalar φ does not develop a vev and U(1) D is broken by other means unspecified here, the coupling of V to the b − s current is generated via a penguin diagram constructed out of a φχ loop. Terms presenting an explicit q 2 dependence give rise to the operator Q bsV 6 . Its Wilson coefficient reads where m φ is the mass of the scalar field and the loop function is [70] F becomes the dominant operator emerging from the φχ loop, thanks to the logarithmic enhancement in Eq. (3.15). The phenomenology of this "split" dark sector was analysed in Ref. [8] and q 2 -dependent operators were considered also in Refs. [1,3,6].
Analogous to this is the case where fermion and scalar multiplet charges are swapped in Eq. (3.11) [8], and UV constructions along similar lines can be considered to generate the DEFT operators involving the muon current in Eqs. (2.3), (2.4).
Note that the DEFT coefficients C bsV 4 , C bsV

6
-similarly to C µ 9(10) in the WET -do not run at the leading order in QCD, since their colour part is simply a vector current [71]. 2

Loop-level
In the case of FIPs as dark matter particles one typically obtains the Wilson coefficients of the WET from "candy" diagrams constructed out of the operators in Eqs. (2.8), (2.9). By making use of a simple cut-off regularisation, one can derive at one loop the Wilson coefficient C µ 9 from a fermion FIP insertion, Note that the only other scale that enters directly the WET coefficients in Eqs. (3.16) and (3.17), apart from the EWSB scale, is the Λ 2 suppression. The exact value of the coefficients will thus depend on the specific nature of the heavy UV completion and, as a direct consequence, the light dark matter mass will have no impact on the flavour anomalies as long as its couplings to the b − s and µ − µ current are mediated by states that can be integrated out of the low-energy theory. 2 Conversely, a DEFT coefficient corresponding to a tensor operator of the type Q bsV 5 = (sσρσPRb) V ρσ , which can also be generated at the loop level by the UV completion, has to be renormalised before being confronted with the low-energy constraints. We leave the treatment of this and other less straightforward cases to future work.

Flavour physics constraints
Global scans combining the constraints from the LFUV ratios R K and R K * with the full set of b → sl + l − branching ratios and angular observables convincingly show the emergence of NP effects in the WET Wilson coefficient C µ 9 , which can be taken alone or in combinations with others. When parametrising the global fit with one NP degree of freedom, the expected size of the Wilson coefficient is  (3.17) show that the numerical value of C µ 9 is quite insensitive to the presence of light dark matter-like spin-1/2 and/or spin-0 states with dimension 6 interactions. We can thus choose to limit our attention to the only case in which the presence of a light FIP has direct impact on the b → s l + l − global fit: the vector FIP.
The constraints on C µ 9 and other WET coefficients derived from the flavour anomalies lead in this case to a typical order-of-magnitude estimate for the DEFT couplings, obtained from a spin-1 FIP exchange in b → sµµ, see Eqs. (3.7), (3.8). One gets where 4 + δ 1,2 indicates the dimension of the corresponding operator, and q is the typical energy exchange captured in the experimental bin. For example, the biproduct of operators Q bsV 6 × Q µµV 4 develops q 2 dependence in the couplings, while this is not the case for a combination of operators of dimension 4.
The presence of a new light vector particle in the spectrum has far reaching consequence for a variety of SM processes. As is discussed in the literature (see, e.g., Refs. [72][73][74]), at FIP masses far below the typical energy scale E of a given SM process, the phenomenology is dominated by the FIP longitudinal polarisation V L as V µ → ∂ µ V L /m V . While this mode cancels out for dimension-6 interactions like Eqs. (2.2) and (2.4), it dominates for the dimension-4 cases, Eqs. (2.1) and (2.3), as they do not correspond to conserved SM fermionic currents. More precisely, the interactions presented in Eq. (2.1) and Eq. (2.3) lead to tree-level flavour violation, weak-isospin violation (since no coupling to neutrinos where included) and an axial-coupling interaction to the SM fermions. We examine in this section the most relevant of the corresponding limits.

B-meson constraints
We first consider the constraints from B-meson physics. B s → µµ Based on a physical process related to the one generating B → K ( * ) µ + µ − transitions, BR(B s → µ + µ − ) can provide a strong constraint on axial muon interactions. After performing a statistical combination of the full LHCb Run 2 data with the global average of Ref. [75], Ref. [33] finds (4. 2) The ratio R Bs between Eq. (4.2) and the SM prediction, BR(B s → µ + µ − ) SM = (3.67 ± 0.15) × 10 −9 [76], can be parametrised in terms of the DEFT. For a vector FIP, exchanged between the currents Q bsV 4 and Q µµV 4 , we replace the WET Wilson coefficients C µ 10 , C µ P , C µ P like in Eqs. (3.8)-(3.10). Substituting into standard expressions [77] one gets where C SM 10 = −4.31 and we have taken the small coupling limit. Equation (4.2) excludes the SM prediction at 2σ, we thus choose to impose the bound at 3 σ to avoid cutting out all the points with a zero or extremely small NP contribution. On the other hand, we observe that Eq.
> 0 , which implies that NP contributions can potentially bring R Bs closer to the measured value. As we shall see in the next section, this sign choice is indeed preferred for m V larger than a few GeV. Conversely, we shall see that in the low mass regime (m V 1 GeV) one gets C bsV 4 C µµV 4 < 0, which produces a strong upper bound at 3σ: Fundamentally different are the cases of the pairs of operators Q bsV 6 , Q µµV 4 and Q bsV 6 , Q µµV 6 for which we find that the NP contribution to B s → µµ vanishes exactly, R Bs = 0, since the term contributing to the axial component is cancelled by the one feeding into the pseudoscalar current. Interestingly, this implies that increasing the sensitivity of the measurement of BR(B s → µ + µ − ) may provide a handle to distinguish scenarios characterised by operators of different dimension. Similarly to the B s → µµ case, the severity of the bound strongly depends on which DEFT operator is generated in the UV completion yielding the b−s coupling. The operator Q bsV 4 leads to an m 2 V -rescaled modification of the WET Wilson coefficients C 2 , C 2 , and C 4 . We obtain, in the small coupling approximation, g., Ref. [78] for an updated value), all computed at the b-quark mass scale, µ b . As was mentioned above, we have not run the SM Wilson coefficients since we compare operators directly at the B s scale. By combining Eq. (4.6) and Eq. (4.5) one gets a strong constraint in the m V 1 GeV range: |C bsV 4 | 2.5 · 10 −6 m V /GeV. Conversely, with UV interactions giving rise to the operator Q bsV 6 , Eq. (4.5) becomes One obtains |C bsV B → K ( * ) X Strong limits arise from the direct measurement of the branching ratios of the V , both in the case of a visible resonance and that of an invisible decay. Resonant decays to visible particles are subject to extremely strong constraints from LHCb in the m V range between 2M µ and M B − M K * [79]. We apply this bound following the numerical recast described in detail in Ref. [8].
To impose bounds on the invisible decay width we use a combination of BaBar results [80,81]. 3 The adopted numerical procedure is described in detail in Appendix A of Ref. [8]. Assuming a large Γ V to invisible final states in the m V 1 GeV range, the bin-dependent bounds on BR(B → K + inv.) induce a strong constraint on the coupling to the hadronic current. In the dimension 4 case we get |C bsV 4 | 10 −8 m V /GeV, whereas for the operator of dimension 6 we get |C bsV 6 /Λ 2 | 5 · 10 −9 GeV −2 .

Lepton sector constraints
The second set of constraints relies instead on probing the FIP interactions with muons, the most relevant of which arise from precision measurements of the W and Z decays.
W and Z decays In order to further probe the parameter space, we perform a simple recast of the result from the ATLAS Collaboration [67] on pp → (referred to as the Z → 4µ search hereafter) in the Z-boson mass window. We implement our effective Lagrangian via FeynRules/UFO [83][84][85] files, then generate hadron-level events pp → µµµµ within the MadGraph5 aMC@NLO platform [86], including the selection cuts from Ref. [67]. 4 Our simple simulation chain leads to a SM fiducial cross-section of 22 fb, in good agreement with the predicted SM result from Ref. [67], 21.2 ± 1.3 fb. We therefore use the final measured result of 22.1 ± 1.3 fb to constrain our parameter space. Note that interference with the SM plays an important role in the final cross-section computation. 5 When the FIP is produced on-shell, an additional constraint on its couplings to muons arises from precision measurements of Drell-Yan dimuon production, as was shown in Ref. [88]. We find, however, that this limit is of the same order as the Z → 4µ ATLAS bound, or even subdominant, for the mass range m V = 1 − 5 GeV in which it was quoted.
Note that the on-shell V production, Z → µµV , participates directly in the Z → 4µ search and introduces a dependence of this limit to the invisible decay width of V .
In the low m V regime, the presence of the massive vector FIP longitudinal mode leads to a E 2 /m 2 V enhancement of various SM decay widths. In particular, in the limit and We then use the world experimental average on the ratio Γ(W → νµ)/Γ(W → νe), 0.996 ± 0.008 [68], to derive the 2 σ limit which holds as long as V decays mostly invisibly. A more conservative limit (which would also apply in presence of electron couplings) is simply to require Eq. (4.8) to be smaller than the total uncertainty on the measured W width. Using Γ W = 2.085 ± 0.042 GeV from Ref. [68], one obtains the 2 σ bound which agrees with Ref. [89].
Resonance search in B-factories The BaBar Collaboration has searched for the process e + e − → µ + µ − V , V → µ + µ − , when the FIP is assumed to have a small width and a mass up to around 10 GeV [90]. In the relevant mass region, our model requires in any case a large invisible width to escape resonant B → K * µµ searches, so that this constraint is subdominant. 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 [91]. While the current limit can hardly compete with other bounds the 2019 run has stored ∼ 10 fb −1 of data so that the search is expected to become rapidly relevant in the near future.
Coupling to electrons We choose not to consider in this work explicit interactions of the FIP with electrons. Nevertheless, a coupling to electrons is generated via the vector kinetic mixing when at least one of the DEFT operators is at dimension 4 [92]. Thus, even below the di-muon threshold, we expect V to decay visibly into e + e − in the absence of an invisible decay channel. Constraints on a light vector FIP coupled to electrons were discussed, e.g., in Refs. [3,5]. In particular, resonance searches at low q 2 in LHCb [13,93] exclude the parameter space relevant for the b → s anomalies, thus requiring the V to decay mostly invisibly (in which case the constraints on B → K + inv. described above apply). 6

Muon anomalous magnetic moment
If the flavour anomalies are explained by a vector FIP exchange, the same FIP can potentially contribute to the anomalous magnetic moment of the muon. The recent confirmation at Fermilab [37] of a 3.3 σ discrepancy between the observed value of (g − 2) µ and the SM renders this possibility all the more enticing and timely. The computation of (g − 2) µ can be enhanced at the 1-loop level if the vector FIP V interacts with the muon current via Q µµV 4 and Q µµV 4 . One gets [95,96] where x µ = M µ /m V and the loop functions read , (4.13) (4.14) Conversely, no significant enhancement is obtained if the coupling of the FIP to the muon proceeds through the operators Q µµV 6 and Q µµV 6 . In those cases the value of (g − 2) µ is suppressed by the UV cut-off Λ yielding, for example, in the case of Q µµV 6 and a similar expression for Q µµV 6 . Its exact numerical value depends entirely on the specifics of the UV completion.
As we have described in Sec. 4.1, the viable range of the DEFT coefficients C bsV 4 , C bsV 6 , corresponding to the operators of the b − s current, is bounded by the constraints from B s -mixing and B → K ( * ) X searches. As a direct consequence, for some points in the numerical scan the typical values of C µµV 4 required for a reasonable agreement with the flavour anomalies is very large, leading to a deviation in δ(g − 2) µ widely exceeding the measured value. A cancellation with the contribution from the axial-vector coupling C µµV 4 may therefore be necessary [2]. In particular, as we will show in the next section, such cancellation is always required for light m V whereas it is not needed for a GeV-scale FIP. Note that including the axial-vector contribution can trigger the bounds from B s → µµ, also discussed in Sec. 4.1, which are particularly strong in the m V 1 GeV range.

Numerical results
In the rest of this work we will consider somehow arbitrarily FIP masses up to 20 GeV, so that m 2 V /Λ 2 EW 1. While there is no specific upper bound on the mass of the vector FIP from the constraints listed in the previous section, we also do not include in our simplified models the interactions between V and the electroweak sector. We leave the complete study of the "transitional" regime for larger masses, up to the electroweak scale Λ EW , for future work.

Fit to the b → s anomalies
We perform a multidimensional fit to the b → s anomalies, including in particular the LFUV ratios R K [9] and R K * [10], the most updated B → K * µ + µ − angular-observable data [21], and the finely binned results for the B → Kµ + µ − and B → K * µ + µ − branching ratios [12,17]. We scan over the following free parameters: m V , γ V ≡ Γ V /m V , C bsV 4,6 , C µµV 4,6 , and C µµV 4,6 . The vector mass m V , expressed in GeV, is flatly distributed either in the [0.01, 2] range, or in the [2,20] one. For the ratio between the V width and its mass, γ V , we employ a logarithmically-flat prior in the range [10 −3 , 0.5]. All the absolute values of the DEFT coefficients have a logarithmically-flat distributed prior in the range [10 −10 , 1]. Since the observables included in these fits depend only on the products C bsV i · C µµV j and C bsV i · C µµV j , but not on the single coefficients, the fit results presented in this subsection are given in terms of DEFT biproducts rather than as a function of individual coefficients.
As was described in detail in Ref. [8], we perform separate fits depending on whether m V lies above or below 2 GeV. In fact, in order to obtain an adequate fit, one has to require that the product C bsV i · C µµV j (C bsV i · C µµV j ) assumes a different sign in each of these two regions. We note that this leads to a negative (positive) C µ 9 (C µ 10 ), in agreement with the global WET fits. Following Ref. [8], we refer to these distinct cases as the high-mass fit and the low-mass fit.
Due to the q 2 dependence of Eqs. (3.7)-(3.10), it is not possible to refer explicitly to the results of the global WET fits that can be found in the literature. One needs instead to confront the DEFT parameter space directly with the experimental data. In order to do so, we employ a custom version of the HEPfit package [97], performing a Markov Chain Monte Carlo analysis by means of the Bayesian Analysis Toolkit (BAT) [98].
We present in Fig. 1(a) the results of the fit in the case of DEFT operators of dimension 4 in both the b − s and µ − µ currents (we dub this as the 44 model hereafter). The plots in red are relative to the high-mass fit, where C bsV 4 · C µµV 4 is required to be negative, the plots in green are relative to the low-mass one, where C bsV 4 · C µµV 4 is required to be positive, and the contours correspond to the smallest regions of 68%, 95%, 99.7% probability. In particular, in Fig. 1(a) we are showing the marginalised 2-dimensional posterior probability density function (2D pdf) in planes of the vector mass m V , width-to-mass ratio γ V , and the products of DEFT coefficients C bsV 4 · C µµV 4 and C bsV 4 · C µµV 4 . In accordance with the respective priors, the posterior pdf's for the width and the coefficients are reported in logarithmic scale. The 2D pdf's involving γ V are only shown for the low-mass fit since they are flatly distributed in the high-mass case and hence not particularly informative.
In the high-mass region of Fig. 1(a) (in red) one can observe a strong correlation between the vector mass m V and the coupling product C bsV it is also possible to observe that the latter can grow from negligible values up to the size of the former (at the 3 σ level) but never get larger. This is consistent with the results of global WET fits [23][24][25][26][27][28][29][30][31][32][33][34][35][36], where NP contributions to C µ 10 -which, as Eq. (3.8) shows, in the DEFT arise precisely from  -are, if non-vanishing, usually smaller in size than the ones to C µ 9 . 7 Note that the solution identified in Ref. [2], corresponding to a mediator with mass m V ≈ 2.5 GeV is disfavoured in our study. This is due to primarily two reasons. On the one hand, the set of observables in Ref. [2] included only the LFUV ratios and P 5 , neglecting other measurements, e.g., the branching ratio BR(B → K * µ + µ − ) [17], which generate large chisquared values for m V ≈ 2.5 GeV. On the other hand, the recently updated values for the angular analysis of B → K * µ + µ − , and hence for P 5 , also push the favoured values for m V above the region highlighted in Ref. [2].
In the low-mass region (in green) C bsV The fit results in the case of DEFT operators of dimension 6 in the b − s current and dimension 4 in µ − µ are presented in Fig. 1(b) (we dub this as the 64 model hereafter). The colour scheme is the same as in Fig. 1(a). C bsV 6 · C µµV 4 is required to be negative in the high-mass region (in red), and positive in the low-mass one (in green). Quantitative differences with respect to the 44 model can be found in this case, driven by the different q 2 dependence of the NP contributions, but we can also highlight a few qualitative similarities, for example, there remains a strong correlation between m V and the coupling product C bsV 6 · C µµV 4 in the high-mass region. However, the high-probability region of the pdf is peaked now at higher m V . In the low-mass region one can see that C bsV 6 · C µµV 4 is bound to a narrow range whereas m V follows a distribution similar to the 44 model. The region of high probability at m V ≈ 1 GeV appears more pronounced in the 64 model than in the 44 model, and it seems to require a comparatively smaller width. This is due to the explicit q 2 dependence of the bsV coupling in the 64 model, which facilitates a good fit to many of the BR(B → Kµ + µ − ) data bins.
Note that a fit performed in the case of DEFT operators of dimension 4 in the b − s current and dimension 6 in µ − µ (46 model ) yields results perfectly analogous to the ones shown in Fig. 1(b), given the equal q 2 dependence of the two scenarios, presented in Eqs. (3.7), (3.8). 7 Exceptions to this statement apply if one advocates for a very conservative treatment of the hadronic uncertainties or, alternatively, for an additional universal NP component in both the electron and muon vector currents. Indeed, if one allows for either one of these possibilities it becomes easier to accommodate NP contributions to the muon axial current of the same size as the ones in the muon vector current (and hence |C µµV [24,25,28,30,33,35,36]. Finally, we present in Fig. 2 the marginalised 2D pdf of the fit in the case of DEFT operators of dimension 6 in both the b − s and µ − µ currents (66 model ). The colour scheme is the same as in Fig. 1. The fit shows many qualitative similarities to the 64 case, with notable differences pertaining mostly to the favoured values of the couplings, due to the different q 2 dependence, and a pronounced high-probability peak at m V ≈ 1 GeV, γ V ≈ 10 −1.8 − 10 −1 , with lower probability when the V mass becomes smaller. Note how the transition to the region where the angular observables exert the greatest relative pull, at very low m V , is much smoother in this model than in the previous cases, indicating that it is not difficult to fit the BR(B → Kµ + µ − ) bin data over a broader mass range, due the additional q 2 -dependence of the muon coupling.

Including all flavour constraints
We can now apply the constraints introduced in Secs. 4.1, 4.2 to the high-probability regions of the fits described above. We present the results for the 44 model in Fig. 3(a), in the plane of the vector muon coupling versus the FIP mass m V . The yellow points in the figure are within the 2 σ regions of the fits to b → s anomalies described in Sec. 5.1. Green points satisfy, additionally, B s −B s mixing and B s → µµ constraints, as well as ATLAS Z → µµV ( * ) at the 2 σ level, as described in Sec. 4.2. We have also ensured that all of them pass the constraints on invisible final states from the BaBar searches [80,81], numerically recast as described in Ref. [8], and the bounds from resonant B → K * µ + µ − anomalies, all green points satisfy the constraints from B → K + inv., B → K * µµ, B s −B s mixing, B s → µµ and Z → 4µ (off-shell+on-shell) at 2σ, applied following the procedure described in the text. We show overlaid the bounds that do not depend on the hadronic sector and/or on Γ V,inv . The blue region covers the parameter space excluded by the multilepton ATLAS search [67] (off-shell), whereas the grey regions show bounds from LFUV in W -decay (see main text). The green band is consistent with the (g − 2) µ measurement [37] at 2 σ, with no fine tuning of the couplings. All points above the green band require a non-zero C µµV 4 , with the dashed grey line indicating a 10% fine tuning with C µµV 4 . The dotted purple line in both panels shows the upper bound from neutrino trident production [59], obtained under the assumption that the coupling of V to neutrinos is of the same size as C µµV 4 . at LHCb, also described in Ref. [8]. We finally overlay the remaining limits, which are independent of Γ V,inv . The green band corresponds to C µµV 4 being consistent at 2 σ with the (g − 2) µ measurement [37].
Some of the constraints applied to the parameter space depend explicitly on both C µµV 4 and C µµV

4
. This is the case, particularly, of the W -decay bound of Eq. (4.8) and forward, and the numerical recast of the Z → 4µ cross-section bound. When these bounds are applied to regions of the parameter space for which C µµV 4 is too large to yield a value of δ(g − 2) µ in agreement with the experimental determination, we tune C µµV 4 as required to bring δ(g − 2) µ down to the measured 2 σ region, cf. Sec. 4.3. 8 The dashed grey line in the figure indicates a fine tuning at the 10% level.
FIPs with mass above the B-meson can easily yield an excellent fit to the experimental data while escaping all current constraints. This conclusion holds for the 44 model, shown in Fig. 3(a), as well as for the 64 model, shown in Fig. 3(b). The favoured parameter space in the 64 model is in agreement with the results obtained in the UV model of Ref. [8]. Incidentally, we find it remarkable that these very compact simplified-model solutions to the b → s anomalies feature also excellent prospect to solve the (g − 2) µ anomaly with no fine tuning required. In the high mass region, the lower limit on the green points distribution arises from the B s -mixing constraints on C bsV 4 and C bsV 6 and the upper limit from the ATLAS Collaboration [67] search for Z → 4µ.
We observe in Fig. 3(b) that for the 64 model the fit yields a large number of points compatible with the various B-physics constraints in the very low-mass region of the parameter space, m V 2 M µ . On the other hand, the required couplings to the muon are typically quite large, to overcome the strong bounds on C bsV 6 /Λ 2 from B → K + invisible. This leads to δ(g − 2) µ exceeding the experimental value and, consequently, a unified solution to all anomalies requires the aforementioned fine tuning of the axial-vector muon coupling against the vector one. The overwhelming majority of these points are excluded by flavour-universality tests in W decays, as measured at the LHC [68].
We show in Fig. 4 the corresponding results for the 66 model. There appear to be extremely few points that can pass all the cuts simultaneously, with the limit from the multilepton ATLAS search [67] dominating the constraints on the muon coupling. This exclusion reflects the fact that the FIP coupling to the muon, C µµV 6 , is of dimension 6 and thus more sensitive to high-energy processes than in the 44 and 64 cases.
We point out, finally, that in the presence of a neutrino coupling, stringent constraints from neutrino trident production apply [59], potentially excluding a common solution for the (g − 2) µ and LFUV anomalies in the high-mass region. For indicative purposes we show as a dotted purple line in both panels of Fig. 3 the upper bound obtained under the assumption that the coupling of V to neutrinos is of the same size as C µµV 4 . This would be indeed the case for SU(2) L -conserving interactions of the V with the lepton doublet, like in the well-known L µ − L τ gauge group. However, it is possible to envision more elaborate UV constructions that could allow one to suppress the coupling to the neutrino 8 As the sign of C µµV 4 is not fixed by this procedure, we choose it in agreement with the requirements from the b → s fits, i.e., opposite to the sign of C µµV with respect to the corresponding charged lepton, in which case the purple line in Fig. 3 would shift upwards. 9 Alternatively, the presence of extra particles and interactions in the UV may introduce contributions to the (g − 2) µ beyond the ones computed in our bottom-up approach. The green band in Fig. 3 would in that case shift downwards.

Summary and conclusions
In this paper we have performed a comprehensive phenomenological analysis of a set of simplified models providing a combined solution to the anomalous magnetic moment of the muon and the flavour anomalies emerged at LHCb in b → sl + l − transitions (both flavour-conserving and lepton-flavour violating). Our simplified models are based on a set of operators of the Dark EFT, or DEFT, where with this term we encompass a broad range of interactions between the SM and one or more light, feebly interacting particles (FIPs) that find their origin in heavy degrees of freedom integrated out of the theory. The most promising models for fitting both anomalies are based on vector FIPs, with the FIP exchange at low energy giving rise to the anomalies in flavour and (g − 2) µ .  1, 1, QD). Yukawa couplings with the SM left-handed muon doublet, Lµ, of the type λEΘ † Lµ, generate a left-chiral coupling g µµ L of V to the muons once U(1)D is broken. The coupling to neutrinos is absent at the tree level. Typically one gets g µµ L ≈ −gDQDλ 2 v 2 D /2m 2 E . The spectrum will have to feature additional scalar singlets/multiplets of SU(2)L, whose scalar potential and dark charges have to be adjusted to cancel the direct V /Z mixing and to raise the mass of the charged scalars above direct LHC bounds.
We have divided our simplified models in categories determined by the explicit momentum dependence of the FIP interaction with the b − s and µ − µ vector currents. For all categories, after performing a global fit to the experimental data we have applied a large set of constraints to the favoured parameter space, extracted from direct measurements of the B → K + invisible branching ratio, B → K * µ + µ − resonant decays, the measurements of BR(B s → µ + µ − ), and B s -mixing. Moreover, we have computed for the first time the strong constraints arising from several precision measurements associated with W -and Z-boson decay, namely the measurement of the Z → 4µ cross section at ATLAS and the decay W → µ + invisible. We find that unified explanations of the (g − 2) µ and b → sl + l − anomalies exist and pass all the constraints in a model based on vector-mediator exchange with DEFT operators of dimension 4 in both the b − s and µ − µ vector currents (model 44, with no momentum dependence of the couplings), and in vector-mediator exchange with dimension 6 in the b − s current (q 2 -dependence of the coupling) and dimension 4 in µ − µ (we call it model 64 ). In both models the typical range of the vector mediator mass is m V 4 GeV and the solution to the (g − 2) µ anomaly does not require any fine tuning of the vector and axial-vector couplings of the V with the muon. Despite a similar mass range for the flavour anomalies in models 44 and 64, the observables B s → µµ and R ∆Ms are sensitive to the momentum dependence of their couplings and thus will be able to discriminate between them.
FIPs with m V > 20 GeV are also likely to provide a good fit to the b → sl + l − anomalies. However, we have worked here under the assumption that the mixing of V with the electroweak sector of the SM can be neglected. As this hypothesis may break down when approaching closely the Z mass from below, possibly leading to additional limits from processes such as Z → γV and others, we leave a proper study of the "transitional" regime between light and heavy NP to future work.
While in model 44 the only viable unified solution lies in the m V 4 GeV range, in model 64 the global analysis of b → sl + l − observables highlights additionally a region of good fit at m V ≈ 0.01 − 1 GeV passing all flavour constraints. However, this region of the parameter space is characterised by a large muon coupling and requires large fine tuning of the vector and axial-vector contributions to satisfy the (g − 2) µ measurement. More importantly, we find that an upper bound on the muon coupling imposed by the Γ(W → µνV ) measured width entirely excludes this region of the parameter space.
In summary, we have highlighted in this work viable solutions for a combined explanation of the (g − 2) µ and b → sl + l − anomalies, based on the simple exchange of a light particle and characterised by a minimal set of assumptions on the heavy UV completion. Our goal is that of providing a self-standing and fairly broad compendium of the low-energy experimental bounds affecting these scenarios. The emerging parameter space regions and coupling strengths can then be used at face value to guide the model-building efforts in the high-energy sector.

A B-meson decays in the DEFT basis
This appendix collects the most important ingredients for estimating the FIP contribution to B-meson related observables, as derived directly from the DEFT operators. We eventually compare the results to the standard WET in order to provide a simple dictionary which can be used to re-adapt existing codes to the case of light particles. We use the convention σ µν = i[γ µ , γ ν ]/2 and as customary define B mesons as containing the b quark instead of the antiquark.

A.1 B → K process
We start with a simple B → K ¯ process. Following the standard procedure, which relies on factorisation of the full amplitude into a hadronic and a leptonic part, one can express theB(p) →K(k) (k 1 ) (k 2 ) amplitude as [99] where N −1 = (4G F / √ 2) V tb V * ts and θ is the angle between the B and the exiting lepton (which is assumed to be a muon in this work) in the B-meson frame. We have furthermore denoted with brackets the contraction of the leptonic current on the vacuum di-muon final state, such as for instance, The coefficients F V , F A , F S , F P , F T and F T 5 are given in Ref. [99] as functions of the Wilson coefficients of the WET basis. The B → K matrix elements of the vector current are usually parametrised as [100] where q = p − k and f 0 (q 2 ), f + (q 2 ) are the form factors. The complete amplitude then follows from contracting the effective vertices corresponding to the DEFT operators of Eqs. (2.1)-(2.4). One obtains where the propagator Π is defined in Eq. (3.4) and we have introduced We can now compare Eqs. (A.4)-(A.7) with the corresponding F i functions expressed in terms of the WET coefficients, shown, e.g., in Ref. [99], to define a dictionary between the DEFT and WET approaches: In the end, we obtain the same results as in Eqs. (3.7)-(3.10), although the coefficients C µ P and C µ P arise in a combination that cannot be disentangled in the WET framework for this particular process.
A.2 B → K * process A similar comparison for B decaying into a vector meson is relatively more involved, due to the larger phase space of the final-state hadrons. Our strategy will be to estimate the relevant amplitude in the helicity basis, following the conventions of Ref. [101]. Employing again the factorisation of the amplitude, one can write where a µ V , a µ A , a S , a P , a µ T R , a µ T L , are coefficients containing the hadronic part of the amplitude specific to the B → K * ll process and we have used once more the bracket notation of Eq. (A.2) to denote the bracket-contracted leptonic currents. Repeating the same procedure as in Appendix A.1 and introducing the helicity amplitudes defined in Ref. [101], one finds [101] to obtain the same correspondence between the WET and DEFT coefficients C µ 9 and C µ 10 as in Eqs. (A.10)-(A.11). Interestingly, we are now able to lift the degeneracy between the C µ P and C µ P , obtaining two additional relations, The above can be compared with the SM prediction to obtain the ratio defined in Eq. (4.2). One can then calculate the R Bs ratio as and check using the WET expression [104] that the link between the DEFT and WET operators can be obtained from the B → K one, as advertised in the main text. Note that the amplitude is suppressed when the FIP mass m V nears M Bs . By expanding in the small coupling limit, one then recovers the main text result Eq. (4.3).
For the case of B s mixing, the relevant four-quark currents correspond directly to the operators of the WET basis with colour indices contracted in pairs (we follow the conventions of Ref. [105]): After simplifying the effective vertices, one obtains We would like to stress that, compared to the standard WET approach, the bag parameters should not be run from a high scale since the amplitude of the process is estimated directly at the B s scale. Note also that the DEFT coefficients used in our results do not run, as the currents charged under colour are (axial)-vector. The experimental observable of interest is the ratio R ∆Ms between the SM-predicted and NP contribution to the mass difference of the neutral B s mesons [68]: where we have defined the ratios of bag parameters at the b scale by In the small coupling approximation, one finally finds (A.42)