The Rumble in the Meson: a leptoquark versus a $Z^\prime$ to fit $b \rightarrow s \mu^+ \mu^-$ anomalies including 2022 LHCb $R_{K^{(\ast)}}$ measurements

We juxtapose global fits of two bottom-up models (an $S_3$ scalar leptoquark model and a ${B_3-L_2}$ $Z^\prime$ model) of \bsll\ anomalies to flavour data in order to quantify statistical preference or lack thereof. The leptoquark model couples directly to left-handed di-muon pairs, whereas the $Z^\prime$ model couples to di-muon pairs with a vector-like coupling. $B_s-\overline{B_s}$ mixing is a focus because it is typically expected to disfavour $Z^\prime$ explanations. In two-parameter fits to 247 flavour observables, including $B_{s/d} \to \mu^+ \mu^-$ branching ratios for which we provide an updated combination and LHCb $R_{K^{(\ast)}}$ measurements from December 2022, we show that each model provides a similar improvement in quality-of-fit of $\sqrt{\Delta \chi^2}=3.6$ with respect to the Standard Model. The main effect of the $B_s-\overline{B_s}$ mixing constraint in the $Z^\prime$ model is to disfavour values of the $s_L-b_L$ mixing angle greater than about $5|V_{cb}|$. This limit is rather loose, meaning that a good fit to data does not require `alignment' in either quark Yukawa matrix. No curtailment of the $s_L-b_L$ mixing angle is evident in the $S_3$ model.


Introduction
Tensions between Standard Model (SM) predictions and measurements of some neutral current flavour-changing B−meson decays persist. We collectively call these tensions the b → sµ + µ − anomalies. For example, various lepton flavour universality (LFU) observables like the ratios of branching ratios R K ( * ) = BR(B → K ( * ) µ + µ − )/BR(B → K ( * ) e + e − ) were, prior to December 2022, measured to be lower than their SM predictions in multiple channels and several bins of di-lepton invariant mass squared [1][2][3]. Measurements of other similar ratios in B 0 → K 0 s ℓ + ℓ − and B ± → K * ± ℓ + ℓ − decays [4] (where ℓ ∈ {e, µ}) are compatible with a concomitant deficit in the di-muon channel as compared to the di-electron channel, although we note that the statistics are diminished in the K 0 s and K * ± channels and the tension is not significant in them (it is around the 1σ level only). SM predictions of all of the double ratios mentioned above enjoy rather small theoretical uncertainties due to their cancellation between numerator and denominator in each case for the di-lepton invariant mass squared category n obs SM quality-of-fit as calculated by smelli2.3.2 with an updated BR(B s → µ + µ − ) constraint (see §4.1). Only the values in bold include the 2022 LHCb measurements of R K ( * ) [7] and values in parentheses are obtained without updating the BR(B s → µ + µ − ) constraint. n obs shows the number of observables in each category. χ 2 SM denotes the χ 2 statistic within each category, p is the p−value of the category, and s is the equivalent two-sided 'number of σ' away from the central value of a univariate normal distribution. The category 'LFU' contains lepton flavour universality violating flavour changing observables such as R that the preselected 247 measurements include several observables which do not involve the (bs)(μµ) effective vertex and which also do not display large tensions, diluting (i.e. raising) the p−value and therefore decreasing s/σ. In an analysis which concentrates more specifically on the fewer b → sµ + µ − anomaly observables on data in 2021, an estimate of s = 4.3σ was made [24]. We parenthetically observe from the table that in the context of our SM global fit to many observables, the update to the BR(B s → µ + µ − ) constraint has only a small effect on the quality-of-fit.
Prior to December 2022, it was well known in the literature that two-parameter new physics models can decrease χ 2 by some 30 or so units [25][26][27][28] resulting in a significantly better fit than the SM. One can take the two parameters to be weak effective theory (WET) Wilson coefficients (WCs), defined via the WET Hamiltonian density where C i here denotes the beyond the SM contribution to the WC, G F is the Fermi decay constant and C SM i denotes the SM contribution to the WC. The two dimension-6 operators O i that can best ameliorate the b → sµ + µ − fits are [25][26][27][28] O 9 = e 2 16π 2 (sγ µ P L b) (μγ µ µ) , where P L is a left-handed projection operator in spinor space and e is the electromagnetic gauge coupling. Here s, b and µ are 4-component Dirac spinor fields of the strange quark, the bottom quark and the muon, respectively. In Ref. [29] (Fig. 1), it was shown that the LHCb 2022 reanalysis of R K ( * ) introduces a mild 1 − 2σ tension between the 'quarks' category and 'LFU' category of observable, if they are interpreted only in terms of C 9 and C 10 . Despite this, we find that models which predict non-zero C 9 and C 10 still provide an improved fit to the combination of current flavour data as compared to the SM. There are two categories of beyond the SM fields that can explain the b → sµ + µ − anomalies at tree level in quantum field theory: leptoquarks, which are colour triplet scalar or vector bosons (with various possible electroweak quantum numbers), and Z ′ s, which are SM-neutral vector bosons. In each category, the new state must have family non-universal interactions, coupling to b L and s L . The observables in which the measurements are in tension with SM predictions involve di-muon pairs. The new physics state should therefore certainly couple to di-muon pairs, but in order to agree with the December 2022 LHCb measurements R K and R K * , one may also couple it to di-electron pairs with a similar strength, although then the scenario may become strongly constrained by LEP constraints [29]. For the coupling to di-muon pairs, fits to flavour data by different groups agree that the new physics state should couple to left-handed di-muons 4 . There is a preference from the data for a coupling to muons and fits to flavour data by different groups agree that the new physics state should couple to left-handed di-muons 5 . (A purely left-handed coupling to muons realises the limit C 9 = −C 10 , which will be a prediction of the benchmark leptoquark model that we examine.) However, depending upon the treatment of the predictions of observables in the 'quarks' category, the different fits may or may not have a mild preference [30] for the new physics field to couple to right-handed di-muons in addition with a similar strength to the coupling to left-handed di-muons (e.g. realising the limit where C 9 ̸ = 0 but C 10 = 0; this will be the prediction of the benchmark Z ′ model which we shall choose) [31]. We display example Feynman diagrams contributing to b → sµ + µ − processes for a scalar leptoquark and for a Z ′ in Fig. 1.
Similarly, many Z ′ models have been constructed which also significantly improve the SM's poor quality-of-fit. Many are based on spontaneously broken anomaly-free U (1) X gauge extensions of the SM that yield such a Z ′ , for various choices of X. Examples include gauged muon minus tau lepton number X = L µ − L τ [46][47][48][49][50], for which the required Z ′ coupling tos L b L is mediated via mixing with heavy vector-like quarks, third family baryon number minus muon lepton number X = B 3 − L 2 [31,51,52], third family hypercharge X = Y 3 [53-  56], and other assignments . An attempt to fit the post-December 2022 data using a lepton-flavour-universal 3B 3 − L data was made in Ref. [29], but the model was found to be in significant tension with LEP Drell-Yan constraints coming from e + e − → ℓ + ℓ − .
The primary aim of this paper is to compare global fits to current data of a leptoquark model with those of a Z ′ model. We pick two simple and comparable benchmark models that have been studied in previous literature. They shall be defined in the next Section. While current tensions in semi-leptonic b → sµ + µ − transitions between measurements and SM predictions can be ameliorated by either a leptoquark or a Z ′ , it is anticipated that the ∆M s observable [78], which parameterises B s −B s mixing and does not exhibit a significant tension with the SM prediction, has the power to discriminate between leptoquark and Z ′ explanations. Because the Z ′ model predicts tree-level contributions to ∆M s , while the dominant contributions from the leptoquark are at one-loop order (see Fig. 2) and therefore suppressed, one generically expects the constraints from ∆M s to be stronger for Z ′ models. However, the extent of the impact of this observable on the fits should be evaluated quantitatively in each model for comparison. A secondary goal of the present paper is to examine the impact of the December 2022 LHCb measurements 6 of R K ( * ) .
For both the leptoquark and the Z ′ , we choose to focus on TeV-scale masses so that the models can avoid the bounds originating from direct searches at the LHC. For TeV-scale masses, the new physics state's contribution to the anomalous magnetic moment of the muon (g − 2) µ is negligible without augmenting the model with additional fields [79] and so we exclude this observable from our fits.
In each model, we do expect there to be additional fields that we (usually) shall assume have a negligible effect upon the b → sµ + µ − observables either because they are too weakly coupled or because they are too heavy. For example, additional heavy fermionic fields are expected in order to generate SM-fermion mixing. If additional fermionic fields are added in vector-like representations of the gauge group 7 , then the model remains anomaly free and large masses for said fermions are allowed by the gauge symmetry. Additional fields could be added in order to explain the current tension in (g − 2) µ , for example additional vector-like leptons [79], where a one-loop diagram with a Z ′ along with additional leptons in the loop can resolve it. Other heavy fields associated with dark matter could also be added. Here, our effective quantum field theory is supposed to be valid at and below the scale of the mass of the Z ′ or leptoquark state, i.e. the TeV scale, and we assume that the effect of other heavy fields on the pertinent flavour phenomenology is negligible.
The parameter space of the B 3 − L 2 Z ′ model model was adapted to b → sµ + µ − data and various phenomenological constraints were applied in Refs. [31,51,52,80]. Some more recent flavour data and LHC Z ′ search limits were employed to constrain the model in Ref. [81]. In that reference, similar constraints were also applied to the S 3 leptoquark model (estimated before with earlier data by Refs. [82,83]), before calculating future collider sensitivities of the two models. In all of these prior works, the parameters of each model were fit to central values for C 9 and C 10 obtained from fits to b → sℓ + ℓ − data, assuming that all other SM effective field theory (SMEFT) operators were zero 8 . The additional constraints, including those from meson mixing and LHC searches, were individually applied to obtain 95% CL bounds. The bottom line of Ref. [81] relevant for the present paper is that, in each model, there is a sizeable parameter space that evades current search limits and satisfies other phenomenological constraints, and that simultaneously explains the b → sµ + µ − anomalies. In the present paper we go beyond the prior studies by performing a side-by-side global flavour fit of each model to identical empirical datasets. All dimension-6 SMEFT operators are included -not just those that generate C 9 and C 10 -including 4-quark and 4-lepton operators induced at tree-level in the Z ′ model, and those operator contributions generated by one-loop renormalisation effects. We then perform a goodness-of-fit test on each model, comparing the best-fit points to quantify any statistical preference of the data for either model.
The paper proceeds as follows: §2, we introduce the S 3 leptoquark model and the B 3 − L 2 Z ′ model. In order to make further progress, some assumptions about fermion mixing must be made and we lay these out first; they are identical for each of the two models. In §3 we present the matching to the WCs of dimension-6 operators in the SMEFT resulting from integrating out the new physics state. The fit to each model will be described by 2 effective parameters extra to the SM. For given values of these parameters, the dimension-6 SMEFT operator WCs are then given and can be passed as input to smelli2.3.2, which calculates the observables and quality-of-fit. Fit results are presented in §4, where we shall observe that the leptoquark model has a similar quality-of-fit as the Z ′ model, indicating that the ∆M s observable was less discriminating than one might have expected. We summarise and conclude in §5.
building fully into the UV, being content to describe the important phenomenology such that experimental assessments of the b → sµ + µ − anomalies can be made. We begin by writing down our conventions as regards fermion mixing, since they are common to both the S 3 model and the B 3 − L 2 Z ′ model.

Fermion mixing conventions
In the fermionic fields' gauge eigenbasis, which we indicate via 'primed' symbols, we write along with the SM fermionic electroweak doublets (2. 2) The SM fermions acquire masses after the SM Brout-Englert-Higgs mechanism through where Y u , Y d and Y e are dimensionless complex coupling constants, each written as a 3 by 3 matrix in family space. For both models that we will consider, some of these Yukawa couplings should be interpreted as effective dimension-4 operators that arise in the SMEFT limit, once heavier degrees of freedom are integrated out; see §2.5 for details. Gauge indices have been omitted in (2.3). The matrix M is a 3 by 3 complex symmetric matrix of mass dimension 1, Φ c denotes the charge conjugate of a field Φ andH := (H 0 * , −H − ) T . We may write H = (0, (v + h)/ √ 2) T after electroweak symmetry breaking, where h is the physical Higgs boson field and (2.3) includes the fermion mass terms the see-saw mechanism via a 6 by 6 complex symmetric mass matrix. Since the elements in m ν D are much smaller than those in M , we perform a rotation to obtain a 3 by 3 complex symmetric mass matrix for the three light neutrinos. These approximately coincide with the left-handed weak eigenstates ν ′ L , whereas three heavy neutrinos approximately correspond to the right-handed weak eigenstates ν ′ R . The neutrino mass term of (2.4) becomes, to a good approximation, be diagonal, real and positive for I ∈ {u, d, e}, and V T ν L m ν V ν L to be diagonal, real and positive (all in ascending order of mass from the top left toward the bottom right of the matrix), we can identify the non-primed mass eigenstates 9 We may then find the CKM matrix V and the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U in terms of the fermionic mixing matrices: (2.8)

Fermion mixing ansatz
To make phenomenological progress with our models, we shall need to fix the V P . We make some fairly strong assumptions about these 3 by 3 unitary matrices, picking a simple ansatz which is not immediately ruled out by strong flavour changing neutral current constraints on charged lepton flavour violation or neutral current flavour violation in the first two families of quark. Thus, we pick V e R = V d R = V u R = V e L = I, the 3 by 3 identity matrix. A nonzero (V d L ) 23 matrix element is required for both the S 3 model and the B 3 − L 2 Z ′ model to mediate new physics contributions to b → sµ + µ − transitions. We capture the important quark mixing (i.e. that between s L and b L ) in V d L as where we use the experimentally determined values for the entries of V and U via the central values in the standard parameterisation from Ref. [85].
Having fixed all of the fermionic mixing matrices, we have provided an ansatz that could be perturbed around for a more complete characterisation of the models. We leave such perturbations aside for the present paper. Having set the conventions and an ansatz for the fermionic mixing matrices, we now introduce the S 3 leptoquark model and the B 3 − L 2 Z ′ model in turn.

Benchmark S 3 leptoquark model
For a bottom-up leptoquark model to explain the b → sµ + µ − anomalies we consider a scalar leptoquark rather than a vector leptoquark, since inclusion of the latter would require an extension of the SM gauge symmetry (which necessarily requires more fields beyond the leptoquark). To fit the b → sµ + µ − anomalies alone, the S 3 ∼ (3, 3) 1/3 leptoquark is a good candidate because it couples to left-handed quarks (and left-handed muons), which can broadly agree with b → sµ + µ − data as described in §1.
The Lagrangian density includes the following extra terms due to the S 3 leptoquark: where x, y ∈ {1, 2, 3} are family indices, a, b, d ∈ {1, 2, 3} are SU (2) L adjoint indices (other gauge and spinor indices are all suppressed) and σ a are Pauli matrices for SU (2) L . Here, λ ϵH3 , λ H3 ∈ R. By charging the leptoquark under a lepton-flavoured U (1) X gauge symmetry [37][38][39][40] one can ban the proton-destabilising di-quark operators (κ xy → 0), and furthermore pick out only couplings of the leptoquark that are muonic: We shall assume (2.11). The leptoquark's couplings to fermions are now specified by a unitnormalised complex 3-vector γ in quark-doublet family space. The factor of λ fixes the overall strength of the interaction of the leptoquark with the SM fermions. We shall moreover assume that γ = (0, 0, 1), which could be enforced with further family symmetries, for example by an approximate U (2) Q global symmetry, or by a particular choice of anomaly-free gauged U (1) X that is quark non-universal. The space of anomaly-free solutions U (1) X symmetries with these properties was explored systematically in Ref. [39] (see Sec 2.3), building on the results of Ref. [86]. For concreteness, we here choose to gauge with an S 3 charge of X S 3 = −2. This U (1) X symmetry allows the coupling to Q ′ 3 c L ′ 2 but no other quark-lepton pairs, and moreover bans the coupling to di-quark operators. We shall discuss the Yukawa sector in §2. 5.
We now perform a global rotation in family space on Q ′ i such that the d L i fields are in their mass basis: whereas in a similar fashion, L i already have the e L i fields in their mass basis: and since (V d L ) 23 ̸ = 0 and (V d L ) 33 ̸ = 0 for θ 23 ̸ = 0 in (2.9), the model possesses the correct couplings to mediate b → sµ + µ − transitions as depicted in the left-hand panel of Fig. 1. Table 2. Representations of fields under the SM gauge factors, which are family universal, together with their representations under the family non-universal gauged U (1) X symmetry on which our Z ′ model is based. We use the minimal integer normalisation for the charges under each U (1) factor and we shall specify q θ ∈ {−1, 1}. All fields are Weyl fermions except for the complex scalar Higgs doublet H and the complex scalar flavon θ.

Benchmark
We extend the SM by a U (1) X gauge group and a SM-singlet complex scalar θ, with field charges as in Table 2. This is the B 3 − L 2 Z ′ model of Ref. [31] that provided a simple bottom-up description of the models in Refs. [52,65] (these possess additional fields), which was shown to explain the b → sµ + µ − anomalies. Note that we neglect any effects coming from U (1) X − U (1) Y mixing. This approximation may be motivated (at tree-level) by further model building, for example by embedding each U (1) Lie algebra generator in the Cartan subalgebra of some semi-simple gauge group which subsequently breaks to 10 The mixing would be set to zero at the scale of this breaking and then generated at one-loop order by running down to the mass scale of the Z ′ . The resulting mixing is small unless the two scales are separated by a large hierarchy.
The U (1) X symmetry is broken by ⟨θ⟩ := v X / √ 2 ∼ O(TeV) and so the Z ′ acquires a tree-level mass where g Z ′ is the U (1) X gauge coupling. The couplings of the Z ′ boson are then 15) in the primed gauge eigenbasis.
Re-writing (2.15) in the mass basis, we obtain the Z ′ couplings to the fermionic fields For θ 23 ̸ = 0, (2.17) implies that the diagram in the right-hand panel of Fig. 1 contributes to the b → sµ + µ − observables. We note here that the s L γ µ Z µ b L coupling is proportional to the gauge coupling multiplied by from (2.9) and (2.17).

The origin of quark mixing
The flavoured U (1) X symmetries that we gauge do not allow the complete set of Yukawa couplings at the renormalisable level, for either benchmark model. For both models, the following quark Yukawa textures are populated at dimension-4: consistent with the SU (2) q × SU (2) u × SU (2) d accidental flavour symmetry of the gauge sector (under which the light quarks transform as doublets while the third family transform as singlets [89,90]). For the gauge symmetry (2.12) in our leptoquark model, the charged lepton Yukawa matrix is strictly diagonal at dimension-4, while for the B 3 − L 2 symmetry of §2.4, the dimension-4 charged lepton Yukawa matrix is: Either of these Yukawa structures for the charged leptons is sufficient to reproduce the observed charged lepton masses.
On the other hand, in order to reproduce the non-zero CKM mixing angles between (lefthanded) light quarks and third family quarks, some of the zeroes in (2.20) must be populated by operators that are originally higher-dimension and that originate from integrating out more massive quantum fields. These mixing angles are therefore expected to be small in this framework, in agreement with observations.
For either of our benchmark models, let us set the charge of the U (1) X -breaking scalar field θ to be q θ = 1. The remaining Yukawa couplings can then arise from dimension-5 operators, as discussed explicitly in Ref. [39] (Section 2.3), where in each term the index i ∈ {1, 2} runs over the light families, and where Λ U,D,Q are effective scales associated with the renormalisable new physics responsible for these effective operators. The first two operators on the right-hand-side of (2.22) set the 1-3 and 2-3 rotation angles for left-handed quarks, which must be non-zero to account for the full CKM matrix, while the second two operators would give rotation angles for right-handed fields, which are not required by data to be non-zero. If these right-handed mixing angles are approximately zero, then the resulting Z ′ will couple only to left-handed not right-handedsb currents. (Sizeable Z ′ couplings to right-handedsb currents are less favoured by the b → sµ + µ − data [91]). A simple UV origin for such operators, which can moreover predict the desirable relatioñ C u i ≈C d i ≈ 0, is to integrate out heavy vector-like quarks (VLQs). Since these VLQs also give other contributions to low-energy observables, notably to B s meson mixing which is a particular focus of this work, we expand upon this aspect of the UV set-up in a little detail. 11 Consider adding a pair of VLQs in the representations This permits the following terms in the renormalisable Lagrangian, where M D,U ≫ v X are mass parameters of the VLQs. 12 Without loss of generality, we can take the coefficients κ D and κ U to be real, while the λ i D,U couplings are all complex. Integrating out the VLQs Ψ D and Ψ U at tree level, at scales Λ D := M D and Λ U := M U respectively, we generate the dimension-5 operators written in (2.22) from Feynman diagrams such as the one depicted in Fig. 3. The (dimensionless) WCs are Once U (1) X is broken by θ acquiring its vacuum expectation value, in either of our benchmark models, these operators match onto dimension-4 up-type and down-type Yukawa couplings suppressed by one power of ϵ U : respectively, thus populating the third column of both quark Yukawa matrices (but not the third row) with 11 We do not consider this discussion to be definitive of our model -there could be other ways to provide the desired patterns -but rather as an existence proof. 12 One might worry that, since the components Ψ U,D R have the same quantum numbers as the light righthanded quark fields, we have neglected dimension-3 mass terms coupling the light right-handed quark fields to Ψ U L or Ψ D L . But such terms can be removed by a change of basis in the UV theory which has no other physical effect. In other words, the fields Ψ U,D R are identified with the linear combinations that couple via a dimension-3 mass term to Ψ U,D L , and MU,D denotes the mass eigenvalue. The light quark fields are identified with the zero mass eigenstates, which only acquire their mass after electroweak symmetry breaking. Figure 3. Feynman diagram that gives the dimension-5 effective up-type Yukawa operators in (2.22), which themselves match onto the 1-3 and 2-3 up-type Yukawa couplings. The corresponding diagrams for down quark Yukawa operators can be obtained by trading every sub/superscript 'U ' for 'D', and swappingH for H.
small couplings, where each '×' denotes a dimension-4 renormalisable Yukawa coupling. It is therefore natural that the CKM angles mixing the first two families with the third are small. 13 In particular, the left-handed mixing angle between the second and third family downtype quarks is sin while the corresponding angle for up-type quarks is κ U λ 2 U v X /M U . Working perturbatively in the small angles, the CKM angle V cb is then (2.28) (Note that this admits a complex phase because the Yukawa couplings λ i U,D are complex). Thus, if there is a mild hierarchy between M D and M U or between λ 2 D and λ 2 U then sin θ 23 -which enters the B-anomaly phenomenology of our models -is predicted to be an order of magnitude or so smaller than |V cb |.
We can invert this argument to estimate the rough mass scale of the VLQs. If we assume that the couplings λ 1,2 U,D are order unity, then we expect the lighter of the two VLQ masses to be of order M ∼ v X /( √ 2|V cb |) ≈ 18v X . In the case of the Z ′ model, we can relate this to the mediator mass via (2.14): M ∼ 18M Z ′ /g X ≫ M Z ′ . Indeed, if one looks ahead to the global fits (see Fig. 6), we see that M Z ′ /g Z ′ is larger than 5 TeV in the 95% CL fit region, meaning that we expect both VLQ masses to be M U,D ≳ 90 TeV .
(2.29) 13 Of course, such a model sheds no light on the mass hierarchies of either up-type or down-type quarks.
The VLQs in these simple models decouple from the low-energy phenomenology. (For the leptoquark model, the scale v X of U (1) X breaking is not tied to the b → sµ + µ − phenomenology at all and so v X , and thus M D,U can be higher still.) 3 Tree-level matching to the SMEFT While recent work has derived one-loop SMEFT matching formulae of models such as the ones we consider [92,93], for our purposes tree-level matching [94] is a sufficiently accurate approximation. 14 In Table 3 we record the tree-level dimension-6 SMEFT coefficients for our benchmark scalar leptoquark model resulting from integrating the S 3 field out of the theory, as computed in Ref. [92]. In Table 4 we record the tree-level dimension-6 SMEFT coefficients for our benchmark B 3 − L 2 Z ′ model obtained by integrating out the Z ′ boson. For both the leptoquark and the Z ′ model, the VLQs introduced in §2.5 to account for the CKM mixing with the third family are significantly more massive than the leptoquark or Z ′ , with masses of at least 90 TeV as discussed above. Their effects on the SMEFT matching are therefore sub-leading. We take care to check the size of their effect in B s meson mixing in §3.1, which we find to be negligible.
WC value Table 3. Non-zero tree-level dimension-6 SMEFT WCs predicted by integrating out the S 3 in our benchmark model (2.10), in units of λ 2 , in the Warsaw basis [96]. Here the EFT matching scale is Λ = M S3 .
In order to make contact with the WCs of the WET, we must first match the B 3 − L 2 model to the SMEFT WCs. Strictly speaking, this should be done at the mass of the Z ′ or the S 3 leptoquark. One then renormalises down to the Z 0 boson mass before matching to the WET. While we shall ignore such renormalisation effects in our analytic discussion (because they are small), one-loop renormalisation effects are fully taken into account in our numerical implementation of both models and thus in our fit results.' The most relevant WET operators 14 One potentially significant effect that we have considered is from the one-loop contributions to four-quark operators in the leptoquark model, which are absent at tree-level but which give the leading contributions to Bs − Bs mixing. The relevant WCs generated by one-loop diagrams are [84,92,95], for general flavours, (C ) kj . Specialising to those relevant to Bs − Bs mixing, We include these one-loop contributions in our analysis. Table 4. Non-zero tree-level dimension-6 SMEFT WCs predicted by integrating the Z ′ out of the B 3 − L 2 Z ′ model, in units of g 2 Z ′ , in the Warsaw basis [96]. Here the EFT matching scale is Λ = M Z ′ .

WC value WC value
for our discussion are those in the bsµ + µ − system, including the operators O 9 and O 10 defined in Eqs. 1.1-1.2, but also the scalar operator O S and pseudo-scalar operator O P where s W = sin θ W with θ W the Weinberg angle, and K = 2 √ 2π 2 /(e 2 G F Λ 2 ) where Λ is the EFT matching scale which is identified with the heavy particle mass in either case.
In our conventions, C 9 > 0 can fit the flavour data well [91] 15 . Substituting in the SMEFT WCs from Table 3, we see that, for the S 3 leptoquark model 15 Note that in some conventions a factor of (VtsV * tb ) is included in the definition of C9, as in Ref. [91], which reverses the sign of Re(C9).
On the other hand, for the B 3 − L 2 Z ′ model, substituting the SMEFT WCs in from Table 4 and using (2.19), A crucial difference between the WCs of the two benchmark models is that in the Z ′ model one turns on not only semi-leptonic operators, but also four-quark and four-lepton operators, whereas a leptoquark gives only semi-leptonic operators in isolation. Since fourquark operators can give contributions to processes such as B s − B s meson mixing (via the process depicted in the right-hand panel of Fig. 2) which do not exhibit a significant disagreement between measurements and SM predictions, it is often stated that leptoquark models provide a 'better explanation' of the data. One purpose of the present paper is to quantitatively assess this claim, by comparing the statistical preference of the b → sµ + µ − data including B s − B s mixing for either model. Before describing the fits in detail, we shall review the dependence of certain key observables, namely B s −B s mixing and the B s → µ + µ − branching ratio, on the WCs. These are important observables that, given sufficient precision, could discriminate between the B 3 − L 2 Z ′ model and the S 3 leptoquark model.
where R loop SM = 1.3397×10 −3 and η(Λ) parameterises renormalisation effects between the scale Λ where the SMEFT WCs are set and the bottom quark mass.
The B 3 − L 2 Z ′ model induces (C (1) qq ) 2323 at tree-level (see Table 4), where the cut-off scale Λ is identified with M Z ′ . η(M Z ′ ) varies between 0.79 and 0.74 when M Z ′ ranges between 1 and 10 TeV [98]. Substituting in for (C (1,3) qq ) 2323 from Table 4 (3.7) The SM prediction of ∆M s is higher than the average of experimental measurements (but roughly agrees with it) and so ∆M s restricts the size of g 2 Z ′ sin 2 2θ 23 /M 2 Z ′ to not be too large [98].
Finally, we note that in both models there are one-loop contributions to B s − B s mixing coming from box diagrams involving exchange of the VLQs. These contributions scale as 1/(16π 2 M 2 U,D ), with further suppression factors coming from suppressed flavour-changing interactions, and so are smaller still than the one-loop contribution coming from integrating out the S 3 leptoquark. Indeed, the largest of these contributions comes from a diagram involving W -boson exchange and the Ψ U L,R fields running in the loop, which scales as sin 2 θ 23 where g is the SU (2) L gauge coupling. We therefore drop these tiny contributions to 4-quark operators in the rest of this paper. As discussed in footnote 14, the S 3 model does not generate (C qq ) 2323 at tree-level; the dominant contribution appears at one loop as in (3.1). Substituting this into (3.6), we obtain (3.8) As expected, the non-SM term in (3.8) is suppressed relative to the equivalent term in the Z ′ model by a factor of [coupling] 2 /(16π 2 ). Numerically, we ran the global fits (see §4) both with and without these 1-loop contributions for the S 3 model, and found the differences in the fit quality (and in particular, in the fit to the ∆m s observable) to be negligible for M S 3 = 3 TeV, the canonical value that we shall take in §4. They could become important, however, for much larger values of M S 3 , since the rest of the observables in the fit prefer λ proportional to M .

New physics effects in
Below, we shall update the smelli2.3.2 prediction of the branching ratio of B s and B d mesons decaying to muon/anti-muon pairs with more recent data. In the presence of new physics parameterised in terms of WET WCs, the B s → µ + µ − branching ratio prediction is changed from the SM one by a multiplicative factor 16 [100] where m Bs , m b , m s and m µ are the masses of the B s meson, the bottom quark, the strange quark and the muon, respectively. Eqs. (3.4) and (3.5) imply that C P , C S , C ′ P , C ′ S , C ′ 9 and C ′ 10 are all negligible both in the B 3 − L 2 Z ′ model and in our S 3 model. Thus, (3.9) implies that, in these two models, Since (aside from small loop-level corrections induced by renormalisation group running) On the other hand, C 10 ̸ = 0 in our S 3 leptoquark model and so BR(B s → µ + µ − ) can significantly deviate from its SM limit. This difference between the models, which has a noticeable impact on the fits (see §4), is just a consequence of fitting the anomalies with a pure C 9 new physics contribution in the Z ′ model verses a C 9 = −C 10 new physics effect in the S 3 leptoquark model.

Fits
We now turn to fits of each benchmark model to flavour transition data. Each simple benchmark model has been characterised by three parameters: Z ′ , say, the missing correction is of order 1/(16π 2 ) log(M Z ′ ). To this good approximation then, each fit is over two effective parameters: one is θ 23 and the other is g Z ′ /M Z ′ or λ/M S 3 , depending upon the model. 17 Direct search limits are a different function of the coupling and mass of the new physics state, however. Currently, the most stringent 95% CL LHC experimental lower limits from the production of di-leptoquarks (which each decay to a jet and (anti-)muons or (anti-)neutrinos) is around 1.4 TeV from the ATLAS collaboration [101]. The B 3 − L 2 Z ′ model has a lower M Z ′ limit of around 1 TeV within the 95% CL parameter space region of a previous b → sµ + µ − anomaly fit [80]. Here, we shall illustrate with a reference mass of 3 TeV for either M Z ′ or M S 3 depending on the model, noting that such a mass is allowed by direct searches. We comment further on the constraints coming from the LHC in §4.4.
We use the default smelli2.3.2 development version constraints upon all observables except for R K ( * ) and {BR(B → µ + µ − ), BR(B s → µ + µ − )}, which have a new joint experimental measurement. We now detail our implementation of the new measurement (we also ran smelli2.3.2 to recalculate all of the covariances in the theoretical uncertainties after this change).

Fit to BR(B
In July 2022, CMS released an analysis of 140 fb −1 of LHC data [5], putting constraints jointly upon BR(B s → µ + µ − ) and BR(B → µ + µ − ). We combine this with the most recent similar measurements from ATLAS [9] and LHCb [102]. We approximate the two-dimensional likelihood from each constraint as a Gaussian in order to easily combine them. Following Ref. [26], we approximate each measurement as being independent, which should be a reasonable approximation because each measurement's uncertainty is statistically dominated. The Gaussian approximations to the measurements and our combination of them are depicted in Fig. 4. The combination corresponds to with a correlation coefficient ρ = +0.28. These are jointly displayed in Fig. 4 Table 5. Quality-of-fit at the best-fit points of the S 3 leptoquark model and the TeV. The first column displays the category of observable. n is the number of measurements in each set. The χ 2 values for each best-fit point are also shown and ∆χ 2 := χ 2 SM − χ 2 and s := sign(∆χ 2 ). The (updated-R K ( * ) ) best-fit parameters for the leptoquark model are λ = 1.3, θ 23 = −0.000467. For the Z ′ model, they are g Z ′ = 0.239, θ 23 = −0.0133. Only the combined category includes the two fitting parameter reduction in the number of degrees of freedom when calculating p. The results using the R K ( * ) LHCb measurements [1][2][3] prior to December 2022 (but including the BR(B s → µ + µ − ) update in §4.1) were performed separately and have different best-fit parameter values to these.
with the two dimensional experimental likelihood, we obtain a one-dimensional pull of 1.6σ, if both branching ratios are SM-like. The recent update of the CMS measurement (which analysed significantly more integrated luminosity than it had previously) has reduced this one-dimensional pull from 2.3σ [26] 18 .

Fit results
We display the best-fit points together with the associated χ 2 breakdown and associated p−value in Table 5. From the table, we note that both models can provide a reasonable fit to the data with p−values greater than .1, in contrast to the SM (see Table 1). We also notice that the improvement of the fit to data over that of the SM is considerable and similar in each case: ∆χ 2 = 3.6 for the S 3 model and ∆χ 2 = 3.6 for the Z ′ model. We see that while the December 2022 LHCb reanalysis of R K ( * ) has reduced the improvement of each model with respect to the SM (evidenced by a lower ∆χ 2 ), the p−value, and therefore the actual quality of the fit, has improved. This is because the SM was a poor fit previously and now it fits the 'LFU' observables well (see Table 1) and because, previously, the new physics models were not very well equipped to explain some of the LFU data (in particular, R K * in the low Q 2 bin) which showed a larger deviation than expected in each model. We pick out some b → sµ + µ − anomaly observables of interest in Fig. 5 to display the pull P i of observable i, defined as where T i is the theory prediction, E i is the experimental central value and S i is the experimental uncertainty added to the theoretical uncertainty in quadrature, neglecting correlations with other observables. From the figure, we notice that P ∆Ms , i.e. the pull associated with B s − B s mixing, is negligibly far from the SM prediction in the S 3 model best-fit point, as expected since the contributions from the S 3 are at the one-loop level (as are SM contributions) and furthermore are mass suppressed. While the Z ′ contribution to ∆M s is mass suppressed, it is at tree-level and therefore enhanced compared to the S 3 contribution. Even though the Z ′ prediction for ∆M s does noticeably differ from that of the SM, it is only by a small amount and the effect on quality-of-fit is consequently small. We see from (3.5) that the new physics contribution to b → sµ + µ − observables is proportional to g 2 Z ′ sin 2θ 23 /M 2 Z ′ , whereas the new physics contribution to ∆M s is proportional to g 2 Z ′ sin 2 2θ 23 /M 2 Z ′ , as can be seen from (3.7), i.e. it involves one additional power of sin 2θ 23 . Thus, in the B 3 − L 2 Z ′ model, decreasing θ 23 but increasing g Z ′ /M Z ′ can then provide a reasonable fit both to the b → sµ + µ − anomalies (which require a sizeable new physics contribution) and to ∆M s , which prefers a small new physics contribution.
We also see that the BR(B s → µ + µ − ) prediction agrees with the discussion in §4.1, i.e. it is approximately equal to that of the SM in the Z ′ model, whereas it is quite different (and slightly ameliorated) in the S 3 model. However, the angular distribution P ′ 5 measurements favour the B 3 − L 2 Z ′ model in comparison to the S 3 leptoquark model. Various observables differ in their pulls between the two models, but in the final analysis, when tension in one observable is traded against others, each benchmark model can achieve a similar quality-of-fit, as summarised in Table 5.

Delineating the preferred parameter space regions
We now delineate the preferred regions of parameter space for each benchmark model, to help facilitate future studies. In Fig. 6, we display a scan over parameter space for each model. The black curves enclose the 95% CL combined region of parameter space, defined as χ 2 − χ 2 min < 5.99, where χ 2 min is the minimum value of χ 2 on the parameter plane. We see in the figure that in the S 3 model, the preferred region reaches to large values of −θ 23 and particular values of λ/M S 3 . In fact, the χ 2 minimum valley is almost flat. On the other hand, B s − B s mixing limits −θ 23 to be not too large in the Z ′ model, producing a less degenerate valley in χ 2 (note the different scales between the left-hand and right-hand panels' vertical axis). We note, however, that values |θ 23 | ∼ O(|V cb |) = 0.04 are easily within the good-fit region for each model, meaning that there is no need for any particular flavour alignment. If θ 23 were constrained to be much smaller, for example, V cb would be required (speaking from the point of view of the initial gauge eigenbasis) to come dominantly from the up-quark sector. For each model, a very thin region of allowed parameter space extends to the right-hand side of each plot for small values of −θ 23 which is not visible (but which can be seen when the ordinate is plotted with a logarithmic scaling, for example).

High-energy constraints from the LHC
We remind readers of the M S 3 > 1.4 TeV bound resulting from an ATLAS di-leptoquark search [101]. Given that the leptoquark is coloured, its pair-production proceeds dominantly by QCD production and so this bound is insensitive to the coupling λ to SM fermions. There are also high transverse momentum LHC constraints on our scalar leptoquark model coming from e.g. the Drell-Yan process pp → µ + µ − , which receives a new physics correction due to the S 3 leptoquark being exchanged in the t-channel. These Drell-Yan constraints, unlike the di-leptoquark search, do not scale only with the leptoquark mass but also scale with the size of its coupling λ to quark-lepton pairs, and so they provide complementary information. Using the HighPT package [103,104], we computed the χ-squared statistic as a function of the leptoquark model parameters (λ, θ 23 ) using the CMS di-muon search [105] implemented within, for a benchmark leptoquark mass of M S 3 = 3 TeV (a value which satisfies the dileptoquark production bound mentioned above). For this mass, we find the 95% CL limit from pp → µ + µ − to be weak, essentially giving no further constraint on the best-fit region to flavour data that we plot in Fig. 6 (left panel)  . Shaded regions are those preferred by the global fit at the 95% confidence level (CL). The yellow region displays the domain |θ 23 | ∈ {|V cb |/2, 2|V cb |} to guide the eye. The region below the green curve in the right-hand plot is incompatible with the ∆M s measurement considered on its own at the 95% CL. There is no such domain in the leptoquark case, explaining the lack of a green curve in the left-hand panel. Our datasets are identical (save for BR(B s → µ + µ − ) as detailed in §4.1 and R K ( * ) ) to those defined by smelli2.3.2 and we refer the curious reader to its manual [6], where the observables are enumerated. The fits have been performed with M S3 = 3 TeV and M Z ′ = 3 TeV, respectively; however each fit is approximately independent of the precise value of the new particle mass (see text). In the right-hand panel, the region to the right-hand side of each vertical dashed magenta contour is excluded for M Z ′ /TeV values below x (where x is the labelled value), to 95% CL [81].
shown in that plot, the model is consistent with pp → µ + µ − for λ < 1.3 or so (λ = 1.3 is further to the right of the plotted region).
For the Z ′ model the present LHC constraints are also relevant. Fig. 6 displays 95% CL lower limits upon M Z ′ from CMS measurements of the di-muon mass spectrum [105]. The values of g Z ′ /M Z ′ have been read off from a re-casting of these measurements [81], neglecting the dependence of the collider bounds upon θ 23 . This is a good approximation because the dominant LHC signal process is bb → Z ′ → µ + µ − , whose amplitude is proportional to cos 2 θ 23 ≈ 1 − θ 2 23 and −θ 23 < 0.18 in the well-fit region. Fig. 6 shows that λ × 3 TeV/M S 3 > 0.03 in the 95% CL region in the S 3 leptoquark model, which can be written M S 3 /TeV < 100λ. We then find an upper bound on λ by considering the width-to-mass ratio of the S 3 , which is Γ S 3 /M S 3 = λ 2 /(8π) [106,107]. Requiring that the S 3 is perturbatively coupled, this width-to-mass ratio should not be too large. Requiring Γ S 3 /M S 3 < 1/3 for perturbativity, we obtain λ < 2.9. Substituting this in, we find M S 3 < 290 TeV, from the combination of the fit to flavour data and perturbativity. This is higher than previous estimates in Refs. [31,81] due in part to the LHCb December 2022 reanalysis of R K ( * ) , which made the minimum value on the abscissa smaller. The bound is also subject to large fluctuations from the implementation of the B s − B s mixing bound due to different lattice inputs for the SM prediction. Since the LHC di-leptoquark searches only constrain M 3 > 1.4 TeV, we see that there is plenty of viable parameter space for a perturbatively coupled S 3 that explains the b → sµ + µ − anomalies. The perturbativity and di-leptoquark searches bound combine to imply that λ × 3 TeV/M Z ′ ≤ 6.2 in the S 3 model. Using a similar argument for the Z ′ , we can first infer from Fig. 6 that the 95% CL region requires g Z ′ × 3 TeV/M Z ′ > 0.03. The width of the B 3 − L 2 Z ′ is approximately Γ Z ′ /M Z ′ = 13g 2 Z ′ /(8π) [31]. Taking Γ Z ′ /M Z ′ < 1/3, we find g Z ′ < 0.80, implying the following upper bound on the Z ′ mass coming from perturbativity and the global fit: M Z ′ < 80 TeV, higher than previous estimates in Refs. [31,81]. Given that the LHC constraints only exclude M Z ′ up to about 4.3 TeV, as indicated by the vertical purple dashed lines on the right-hand plot of Fig. 6, we see that there is plenty of viable parameter space for a perturbatively coupled B 3 − L 2 Z ′ that explains the b → sµ + µ − anomalies. We see from Fig. 6 and Ref. [81], that in the region of 95% CL fit to flavour data, M Z ′ ≥ 2.0 TeV. Combining this and perturbativity, we find that g Z ′ × 3 TeV/M Z ′ ≤ 1.2 in the Z ′ model.

Summary
We have contrasted two bottom-up beyond the SM physics models that can significantly ameliorate the b → sµ + µ − anomalies in global fits 19 . Since one is a leptoquark model and one a Z ′ model, a naive expectation is that the Z ′ model is disfavoured by the fact that measurements of B s − B s mixing are broadly compatible with the SM. We have shown that, contrary to the naive expectation, the B s − B s mixing prediction of the best-fit point of each model is similar, being close to the SM prediction; ∆M s is therefore not a key discriminator. That said, as Fig. 6 shows, the B s − B s mixing constraint limits the B 3 − L 2 model's value of θ 23 in the Z ′ model, albeit within a natural range of values that comfortably includes θ 23 ≈ |V cb |, whereas no such preference is strongly evident for the S 3 model. The Z ′ model that we picked differs in other ways to the leptoquark model in that it couples to di-muons through a vector-like coupling 20 , whereas the leptoquark couples through a purely left-handed coupling. This shifts the predictions for some of the other flavour observables around but in the end the fits are of a similar quality to each other; the improvement in quality-of-fit with respect to the SM of ∆χ 2 = 3.6, 3.6 in the S 3 model and in the Z ′ model, respectively.
The LHCb December 2022 reanalysis of R K ( * ) pushes the fits to smaller new physics couplings to fermions for a given b L − s L mixing angle. One unexpected consequence of including the reanalysed measurements was that the fit to each new physics model improves, as Table 5 displays. The improvement with respect to the SM decreases as expected, however each new physics model (with two fitted parameters) still improves upon the SM χ 2 by some 12 units. This realisation should be tempered with the risk that the remaining discrepant b → sµ + µ − observables could potentially receive unaccounted-for long-distance contributions from charm penguins 21 . In Ref. [108], for example, the discrepancy with SM predictions is ascribed completely to such contributions and is used to parameterise them.
As previously mentioned, both in the B 3 − L 2 Z ′ model and in the S 3 leptoquark model, B−physics fits are insensitive to scaling the coupling and the mass by the same factor. However, direct searches for the new hypothesised particles do not exhibit this scaling symmetry. The S 3 model di-leptoquark search constraints only depend sensitively on the mass of the leptoquark [83], since the production is governed by QCD and thus via the strong coupling constant, which is of course known from experimental measurement. Single leptoquark production in the S 3 model does depend sensitively both upon the mass of the leptoquark and upon its coupling [82]. Direct B 3 − L 2 Z ′ search constraints also depend upon the coupling and the Z ′ mass [31]. However, the scaling symmetry of the fits allows us in §4.3 to present complete fit constraints upon the three parameters of each model in a two-dimensional plane. We hope that this will facilitate future direct searches for either explanation of the b → sµ + µ − anomalies.