Large Hadron Collider Constraints on Some Simple $Z'$ Models for $b\to s\mu^+\mu^-$ Anomalies

We examine current Large Hadron Collider constraints on some simple $Z'$ models that significantly improve on Standard Model fits to $b\to s\mu^+\mu^-$ transition data. The models that we consider are the 'third family baryon number minus second family lepton number' $(B_3-L_2)$ model and the 'third family hypercharge' model and variants. The constraints are applied on parameter regions of each model that fit the $b\to s\mu^+\mu^-$ transition data and come from high-mass Drell-Yan di-muons and measurements of Standard Model processes. This latter set of observables place particularly strong bounds upon the parameter space of the $B_3-L_2$ model when the mass of the $Z'$ boson is less than 300 GeV.


Introduction
Certain experimental measurements of particular B hadron decays are currently in tension with Standard Model (SM) predictions.
The ratios of branching ratios BR(B → K ( * ) µ + µ − )/BR(B → K ( * ) e + e − ) [1][2][3], BR(B s → µ + µ − ) [4][5][6][7][8], angular distributions in B → K * µ + µ − decays [9][10][11][12][13][14] and the branching ratio BR(B → φµ + µ − ) [15,16] are some examples, which we dub collectively as 'b → sµ + µ − anomalies'. Some of these observables have small theoretical uncertainties in their SM predictions, whereas others have more sizeable theoretical uncertainties. No individual measurement is yet in sufficient tension to claim unambiguous 5σ evidence of new physics. However, collectively, the tensions point to the same conclusion even when theoretical uncertainties are taken into account: that a beyond the SM (BSM) contribution to a process connecting a left-handed bottom quark b L , a left-handed strange quark s L , a muon µ − and an anti-muon µ + is preferred (together with the anti-particle copy of the process). 1 One recent estimate puts the combined global significance at 3.9 standard deviations [18] once the look elsewhere effect and theoretical uncertainties are taken into account. Recent fits to single BSM effective field theory operators broadly agree with 1 Although their tensions with SM predictions are mild, recent experimental determinations of the ratios of branching ratios BR(B 0 → K 0 S µ + µ − )/BR(B 0 → K 0 S e + e − ) and BR(B + → K ( * )+ µ + µ − )/BR(B + → K ( * )+ e + e − ) bolster this conclusion further [17]. each other [19][20][21]: a correction to the coefficient of the effective field theory operator (b L γ µ s L )(μγ µ P X µ) in the Lagrangian density can greatly ameliorate the tension between the predictions and the measurements. Here, P X is a helicity projection operator. According to the fits, both P X = P L (a coupling to left-handed muons) as well as P X = 1 (i.e. a vector-like coupling to muons) work approximately as well as each other, but P X = P R is disfavoured. The fits are typically done in the approximation that the energy scale (here given by the bottom meson mass m B ) relevant to the measured observables is much smaller than the scale of new physics producing the operator. A BSM contribution to such effective field theory operators can be generated by the tree-level exchange of a massive electrically-neutral gauge boson (dubbed a Z ), as depicted in Fig. 1, if it has family dependent couplings. In particular, to generate the effective field theory operator required, it should have a coupling both to muon/antimuon fields and to b L s L + s L b L . In the fits of a BSM effective field theory operator to the measurements, terms of order m 2 B /M 2 Z are implicitly neglected in the effective field theory expansion, where M Z is the mass of the Z boson.
In order to explain the b → sµ + µ − anomalies with such a Z , the product of the Z -coupling to b L s L + s L b L and its coupling to µ + µ − divided by M 2 Z has a range of values that does not include zero. It is then of interest to ask if such a Z boson could be produced directly at high energy proton-proton (pp) colliders and detected so as to provide a smoking-gun signal of the model. The most obvious decay channel is into µ + µ − : the b L s L +s L b L coupling is strongly bounded from above by B s − B s mixing constraints, where the data and SM prediction are consistent with each other, but where tree-level Z contributions from the diagram in Fig. 2 are predicted. It was shown in Refs. [22,23] that one expects a 100 TeV future circular hadron-hadron collider to be able to cover much of the available parameter space of generic toy Z models in the µ + µ − channel in parameter regions consistent with the b → sµ + µ − anomalies. Such searches have been carried out by the Large Hadron Collider (LHC) general purpose experiments ATLAS and CMS, but so far no significant signal has been found. Constraints from a 139 fb −1 13 TeV pp ATLAS high-mass Drell-Yan di-lepton search were placed upon the toy models in Ref. [24]. The parameter spaces of toy models which fit the b → sµ + µ − data were only weakly constrained by the search. However, a more complete model, the Third Family Hypercharge (Y 3 ) Model [25], was found to be more strongly constrained by the search; in particular, it was found that M Z > 1.2 TeV, calculated in the tree-level limit with an on-shell Z (and not including associated production with a jet). Here, the fit to the b → sµ + µ − data was rather crude: dominant tree-level effects from a SM effective field theory operator were included, and all renormalisation effects were ignored.
In the Y 3 model as well as other typical explicit models, massive Z gauge bosons originate from spontaneously broken U (1) gauge symmetries. Several other models with additional U (1) gauge symmetries 2 and family dependent charges have been proposed to obtain a Z with the correct properties to explain the b → sµ + µ − anomalies . A variant of the Y 3 model, the Deformed Third Family Hypercharge (DY 3 ) was introduced in Ref. [61] where limits from the ATLAS high-mass Drell-Yan di-lepton search were placed upon it. However, the global fits to the Y 3 model and DY 3 model have subsequently changed significantly from the inclusion of precision electroweak observables, which tend to pull the fit more toward the SM limit [62]. A new variant (DY 3 ) was introduced in Ref. [62] and a fit of the model to electroweak and b → sµ + µ − data was performed. The aforementioned ATLAS highmass Drell-Yan di-lepton search has also constrained the parameter space of the baryon number minus second family lepton number (B 3 − L 2 ) model [45,46]. In the B 3 − L 2 model, in contrast to the third family hypercharge type models, a region of parameter space with M Z < 300 GeV was found to simultaneously explain the b → sµ + µ − anomalies and pass all of the other relevant experimental constraints [46]. Such low values of M Z are unfeasible in third family hypercharge type models because the fit to electroweak data would become too poor.
Our aim here is to update LHC constraints on the parameter space regions of the Y 3 , DY 3 , DY 3 and B 3 − L 2 models that fit the b → sµ + µ − anomalies. For the first three models in this list, the good-fit parameter space has changed significantly due to the inclusion of electroweak precision observables in the fit of Ref. [62] (pertinent new data from LHCb were included and the theory predictions were improved to include a proper matching to the SM effective field theory calculation and renormalisation effects, but these had a less dramatic effect). We will check the constraints coming from measurements of SM-predicted processes on the DY 3 , DY 3 and B 3 − L 2 models for the first time. Notably, we shall show that the M Z < 300 GeV region of the B 3 − L 2 model is strongly constrained by such measurements. The calculation of high-mass Drell-Yan di-lepton search constraints upon the DY 3 model is also new. The calculation of limits in the higher mass régimes of the models is more accurate than previous determinations in the literature because it includes associated production of a jet. 3 The paper proceeds as follows: the models are defined in §2 along with a characterisation for each model of the parameter region that fits the b → sµ + µ − anomalies. In §3, we then go on to introduce the various measurements used to constrain the models, which are calculated by the Contur2.2.0 computer program [63,64]. We present the resulting constraints in §4 before summarising and discussing them in §5.

Models
In this section we shall introduce the four models under study in the present paper. The chiral fermion content of each model is that of the SM augmented by three righthanded neutrinos (SM+3ν R ); these are added in order to obtain neutrino masses. Each model extends the SM gauge group by an additional U (1) gauge group factor and each is anomaly-free. In order to distinguish it from U (1) Y , we shall call the additional is broken by the vacuum expectation value θ of a SMsinglet complex scalar θ, the flavon, which has a non-zero U (1) X charge X θ . Then, the U (1) X gauge boson X µ acquires a mass via the Brout-Englert-Higgs mechanism, where g X is the U (1) X gauge coupling. The family dependent charges of the other fields are given in Table 1 for the four models under study. We use the following notation for the representation of the chiral fermion fields under (SU (3), and L Li := (ν Li , e Li ) T ∼ (1, 2, −1/2) and the L and R suffix refers to left and right-handed chiral fermions, respectively. The complex scalar Higgs doublet field transforms as H ∼ (1, 2, 1/2). In what follows, we denote 3-component column vectors in family space with bold font, for example u L := (u L1 , u L2 , u L3 ) T . Each of our models then has X µ -fermion couplings in the Lagrangian density as dictated by the charges of the fermions. For each fermion species we have a Lagrangian density term where we have defined the three-by-three hermitian matrices The V F matrices are 3 by 3 unitary matrices that transform F from the (primed) weak eigenbasis to the (unprimed) mass eigenbasis, i.e. F := V F F . The CKM matrix V and the PMNS matrix U are then predicted to be respectively.
In order to specify the models further for phenomenological investigation, one must make some assumptions about the V F matrices. They are taken to be consistent with (4) once empirical inputs are taken for the central values of the entries of U and V . For simplicity, we take the 3 by 3 identity matrix. V ν L = U may be fixed by using empirical inputs for U (however, neutrinos shall play no further role in our study). We shall require V d L = I 3 , since we require a coupling between d L2 and d L3 and the X µ boson in order to explain the b → sµ + µ − anomalies. A 'standard parameterisation' [65]  where s ij := sin θ ij and c ij := cos θ ij , for angles θ ij , δ ∈ R. θ 23 was allowed to vary, since it is this coupling that controls the X µ coupling to b L s L +s L b L . The assumptions on the angles and phase in V d L in the fits we use are different for the third family hypercharge type models [62] and the B 3 − L 2 model [54], so we consider each in turn.

Third family hypercharge models
In this subsection, we shall characterise recent global fits of the three third family hypercharge type models under investigation (Y 3 , DY 3 and DY 3 ) to b → sµ + µ − data [62]. Aside from the effects of varying θ 23 , we set V d L equal to the CKM matrix. In more detail, the other angles and phases in V d L are set equal to the central values of the parameters of the CKM matrix inferred by experimental measurements, i.e. s 12 = 0.22650, s 13 = 0.00361 and δ = 1.196 [65].
X H = 0 in third family hypercharge type models (as Table 1 attests) implying that they predict 'Z 0 -Z mixing' (i.e. mixing between the X µ boson, the hypercharge gauge boson and the electrically neutral SU (2) L boson) and thus they affect the prediction of precision electroweak observables. Ref. [62] went on to fit the third family hypercharge type models to 219 electroweak and b → sµ + µ − data. The effective field theory calculation of observables in the electroweak sector implicitly misses relative correc- and so is only valid if M X M Z 0 , the mass of the Z 0 gauge boson. After electroweak symmetry breaking, the mass of the gauge boson is corrected to follow All three models (Y 3 , DY 3 , DY 3 ) provided a much better fit than that of the SM: improving ∆χ 2 by 33-43 Fig. 3: Global fits of third family hypercharge models from Ref. [62] for M X = 3 TeV. The shaded regions are the 95% CL fit regions for the model, as according to the legend. The inner contours within each shaded region enclose the 68% CL region. The parameter point of best-fit is labelled by a dot in each case. To a good approximation, the fits are independent of M X , provided that it is much larger than M Z 0 and provided that g X is scaled proportional to M X , as implied by the abscissa [62]. We show our characterisation of the good-fit region by the dashed line in each case, detailed in (6) and Table 2. units depending upon the model, for two fitted parameters (θ 23 and g X /M X ). The 95% confidence level (CL) region of each model is shown in Fig. 3 as a shaded region, as according to the legend. We wish to reduce the three independent parameters (g X , M X and θ 23 ) down to two in order to display search constraints in terms of twodimensional plots. To this end, we characterise each fit in Fig. 3 by a dashed curve. The equation of each curve is parameterised by where and a and b were 'fit by eye' for each model. The values taken, along with the domain of good fit, are displayed in Table 2 for each model. It will suit us below to scan in x and M X ≈M Z , constraining θ 23 to satisfy (6) in order to stay within the region of good fit for the domain of x given in Table 2. We shall refer to this region, where flavour data (and electroweak data for the third-family hypercharge type models) are within the 95% CL, as the 'favoured region' of each model.   Fig. 4: Tree-level Feynman diagram of a Z contribution to the neutrino trident process. N represents a nucleon.

B 3 − L 2 model
In the B 3 − L 2 model [45,46], since X H = 0, there is no Z 0 − Z mixing at tree-level and so M Z = M X . Electroweak precision observables then follow the SM predictions, to a good approximation. We may thus entertain lower values of M Z , since neither the theoretical consistency constraint nor the need to avoid large corrections to SM predictions of electroweak observables (both of which imply M X M Z 0 in third family hypercharge type models) apply to the B 3 − L 2 model. Ref. [46] showed that, as well as possessing viable M Z > 1 TeV parameter space, the model has a region of parameter space within the domain 20 GeV < M X < 300 GeV which simultaneously evades other constraints whilst providing much improved fits to b → sµ + µ − data.
The B 3 − L 2 model was matched to fits of b → sµ + µ − data in Ref. [54]. The assumptions on V d L were equivalent to taking s 12 = s 13 = δ = 0 in (6) and it was found, after matching to the b → sµ + µ − fit, that θ 23 satisfies (6) with the values shown in Table 2. The upper bound on the domain of x comes from measurements of the trident process, which roughly agree with SM predictions and so cannot receive large corrections from the Zmediated process shown in Fig. 4. The lower bound on the domain of x comes from measurements of B s − B s mixing, which bounds the contribution coming from the process in Fig. 2. We note in passing that the fit we match to in the B 3 − L 2 model taken is less sophisticated than the global fits used for the third family hypercharge type models; the fit to the B 3 − L 2 model was at the tree-level and did not include renormalisation group effects. In contrast to the third family hypercharge type models, the predictions of the B 3 − L 2 model for electroweak observables are identical to those of the SM and so they were not included in the fit.

LHC Constraints from Contur
The Z -fermion interactions of the four models that we use have been encoded into UFO format [66] 5 by using Feyn-Rules [67]. Currently, the flavon is neglected in the files, however Z 0 − Z mixing effects [25] have been included to leading order in (M Z 0 /M Z ) 2 for the third-family hypercharge type models.
The UFO files are then passed to the Herwig7.2.2 [68,69] event generator, which calculates the total width and branching fractions of the Z , and generates full final-state simulated pp collision events 6 . These events are passed through the Rivet3.1.5 [70] library of analyses. This constitutes a 'signal injection' of the putative BSM contribution to several hundred differential cross-sections measured at the LHC and stored in HepData [71]. Contur then evaluates whether this BSM contribution would have been visible given the experimental uncertainties, and if it should have been 'seen already', derives an exclusion probability. The approach is described in more detail in [64].

Results
We find that high-mass Drell-Yan LHC searches into a dimuon final state provide the most constraining bound at large values of M Z , irrespective of the model. The dominant partonic Z production processes are shown in Fig. 5, where it is emphasised that one requires a bb partonic initial state from the LHC proton pairs pp. The calculated Z production cross-section is thus dependent upon which parton distribution function (PDF) set is used, since the b(b) content of the proton differs from PDF-set to PDFset at relatively high values of the ratio of partonic centre of mass energy to pp centre of mass energy 7 For the results discussed below, the default Herwig 7.2.2 choice of CT14 [72] was used.
The exclusion limit for the Y 3 model is shown in Fig. 6. The sensitivity is dominated by the di-muon channel in the high-mass ATLAS Drell-Yan di-lepton search, with the CMS measurement [73] (which uses only 3.2 fb −1 of integrated luminosity) also contributing. The data barely impinge on the allowed parameter space, with only a small region at low mass and large coupling disfavoured. Here, we see that the favoured region g X < 0.2 is not constrained at the 95% CL in the parameter space considered. The colours indicate the bound giving the dominant sensitivity as in the key below. The region above the white solid line is excluded at the 95% CL and the region above the white dashed line is excluded at the 68% CL. The region favoured by the fits is below the blue dashed line.

ATLAS and CMS high-mass Drell-Yan
For the DY 3 model, Fig. 7 shows that the sensitivity is somewhat greater, extending to M Z ≈ 1.2 TeV for the highest values of g X × 1 TeVM Z considered, with the same two datasets contributing. At low M Z , the favoured region 0.1 ≤ g X × 1 TeV/M Z ≤ 0.2 is constrained by the ATLAS high-mass Drell-Yan ll search at the 95% CL.
The exclusion is stronger still for the DY 3 model, shown in Fig. 8. At high coupling, masses up to M Z ≈ 1.3 TeV are excluded, although in the region favoured by the fits, only the low M Z region is impacted.
The comparative strength of the bounds in the three third family hypercharge models can be understood in terms of the size of the absolute additional U (1) charges of the left-handed and right-handed muons. For a given point in parameter space, the larger the absolute value of the charge, the larger is BR(Z → µ + µ − ) and so the bound from high mass Drell-Yan di-lepton searches is concomitantly stronger 8 . As a careful reading of Table 1 allows, whilst the quark charges are identical for the three models in question, the muon charges are largest for the DY 3 model, next largest for the DY 3 model and smallest The colours indicate the bound giving the dominant sensitivity as in the key below. The region above the white solid line is excluded at the 95% CL and the region above the white dashed line is excluded at the 68% CL. The region favoured by the fits is below the blue dashed line.

ATLAS and CMS high-mass Drell-Yan
For the B 3 − L 2 model, the range 200 < M Z < 1000 GeV is excluded for all allowed couplings, as shown in Fig. 9. In this case di-lepton-plus-photon final states [74-76] also contribute. At high couplings the limit extends up to around 3.7 TeV. However, as previously mentioned, in this model an open parameter window at low masses M Z < 300 GeV also exists.  Fig. 9: Exclusion in the high mass region of the B 3 − L 2 model. The colours indicate the bound giving the dominant sensitivity as in the key below. The region above the white solid line is excluded at the 95% CL and the region above the white dashed line is excluded at the 68% CL. In the region above the black dashed line, the width of the Z is more than M Z /3, and so the perturbative cross-section calculation becomes unreliable. The region favoured by the fits is between the blue dashed lines.

ATLAS and
CMS high-mass Drell-Yan ATLAS γ In Fig. 10, we examine this low-mass region. The picture in terms of contributing analyses is more complex. Dilepton measurements at the Z 0 pole, particularly those in Refs. [77,78], exclude much of the region for g X × 1 TeV/M Z > 1. Measurements targeted at W bosons decaying leptonically [79-81] play the dominant role around M Z = 300 GeV for high couplings. Photon-plus-di-lepton measurements again contribute, especially at higher masses and lower couplings. Lower mass di-lepton measurements [82] contribute for M Z < M Z 0 . In the region where the sensitivity runs out, the inclusive four-lepton measurement [83] also contributes, due to a clean (but very low cross-section) (Z Z 0 )-production contribution. With more integrated luminosity this observable would become more sensitive.
The fact that the cross-section for this model is large enough to significantly distort the expected distributions, even in the presence of a large SM cross-section, means that the model is disfavoured over the majority of the previously open parameter window, leaving only a small region at low mass and low values of g X × 1 TeV/M Z still allowed. We note that towards the left-hand side of the plot, the accuracy of the previous fit of the model to flavour data is called into question since unaccounted-for relative corrections O(m 2 B /M 2 Z ) become sizeable. For an example point in parameter space for each model, we display some relevant cross-sections in Table 3. The region above the white solid line is excluded at the 95% CL and the region above the white dashed line is excluded at the 68% CL. In the region above the black dashed line, the width of the Z is more than M Z /3, and so the perturbative cross-section calculation becomes unreliable. The region favoured by the fits is between the blue dashed lines.  Table 3: Example leading order matrix element crosssections calculated by Herwig for some parameter points for the four models considered. The cross-sections quoted are pp cross-sections at 13 TeV centre-of-mass energy, in fb. σ Z +q,g is the cross-section for associated Z -quark production plus the cross-section for associated Z -gluon production, with a minimum transverse momentum of 20 GeV for the quark or gluon.

Summary and Discussion
We have calculated LHC bounds upon four models that have been fitted to b → sµ + µ − anomalies. Each of the models includes an electrically neutral, massive Z gauge boson which has family dependent couplings to SM fermions, the most important for our discussion being the couplings to µ + µ − and tobs + bs. 9 Three third family hypercharge type models (the Y 3 , DY 3 and DY 3 models) were recently fit [62] to b → sµ + µ − data, which is in tension with SM predictions. The fits included electroweak data, since the third family hypercharge models alter the SM predictions of the electroweak observables. In the present paper, we have calculated LHC constraints upon the parameter space of these three models that fit the b → sµ + µ − data. The LHC constraints upon the three third family hypercharge models are much weaker than those calculated previously in Refs. [24,61]. This is due to the fact both the electroweak data and recent b → sµ + µ − data pushed the fit towards the SM limit, preferring smaller values of g X × 1 TeV/M Z as compared to previous preferred values. Our calculation is at a higher level of precision compared to previous estimates, since it includes the effects of associated production and renormalisation of the SM effective field theory in the fit. Neither the DY 3 model nor the DY 3 model had been checked previously against measurements of SM-predicted processes. However, for the parameter regions of interest where M Z > 300 GeV, these are not as constraining as the high-mass Drell-Yan di-lepton searches, as Figs. 6-8 show. In general even the high-mass Drell-Yan limits upon the third family hypercharge models are not very constraining; one cannot quote a lower bound on M Z independent of the coupling for any of the three models.
In contrast, the B 3 − L 2 model is more tightly constrained by the high-mass Drell-Yan di-lepton searches at the LHC, as Fig. 9 shows: M Z >1 TeV in the favoured region. The measurements of SM-predicted quantities (calculated here for the first time) have a large impact on the low M Z window of the B 3 − L 2 model, as shown in Fig. 10. Differential cross-sections in an ATLAS µµ+jet analysis play a particularly important role at low M Z . The constraints in the low M Z window are significantly stronger than those calculated previously in Refs. [46,54].
Generally, measuring τ leptons is more difficult than measuring muons in LHC experiments. However, at higher transverse momenta, the hadronic τ energy resolution is expected to improve [86], since it is dominated by calorimetery, whereas the muon resolution degrades, due to the lower curvature in the trackers [87,88]. Higher values of M Z lead to final state particles at higher transverse momenta on average, so τ leptons may become relatively more important. We note that while no relevant measurements of τ final states are currently available in Rivet, for the Y 3 model, the branching fraction Z → τ + τ − is around 10 0.3, compared to Z → µ + µ − ≈ 0.075. For DY 3 However, in order to satisfy other experimental constraints, simple Z models such as those deployed in the present paper are forced into a parameter space where the beyond-the-SM contribution to aµ from the Z is too small to explain the discrepancy and so further model building involving additional fields and/or additional Z couplings [85] would be required to explain the measured value of aµ. 10 We quote representative values of the branching fractions in the régime where M Z is much greater than twice the top quark mass. The predicted branching fractions we quote here all increase for lower values of M Z . the corresponding branching fractions are 0.4 and 0.1, and for DY 3 0.35 and 0.22, suggesting that τ measurements could make an important contribution in future, despite the additional experimental challenges involved.
For each of the four models that we analyse, an appreciable portion of parameter space remains where future LHC analyses may search for (and hopefully find) a signal.