The Case for Future Hadron Colliders From $B \to K^{(*)} \mu^+ \mu^-$ Decays

Recent measurements in $B \to K^{(*)} \mu^+ \mu^-$ decays are somewhat discrepant with Standard Model predictions. They may be harbingers of new physics at an energy scale potentially accessible to direct discovery. We estimate the sensitivity of future hadron colliders to the possible new particles that may be responsible for the anomalies: leptoquarks or $Z^\prime$s. We consider luminosity upgrades for a 14 TeV LHC, a 33 TeV LHC, and a 100 TeV $pp$ collider such as the FCC-hh. Coverage of $Z^\prime$ models is excellent: for narrow particles, with perturbative couplings that may explain the $b$-decay results for $Z^\prime$ masses up to 20 TeV, a 33 TeV 1 ab$^{-1}$ LHC is expected to cover most of the parameter space up to 8 TeV in mass, whereas the 100 TeV FCC-hh with 10 ab$^{-1}$ will cover all of it. A smaller portion of the leptoquark parameter space is covered by future colliders: for example, in a $\mu^+\mu^-jj$ di-leptoquark search, a 100 TeV 10 ab$^{-1}$ collider has a projected sensitivity up to leptoquark masses of 12 TeV (extendable to 21 TeV with a strong coupling for single leptoquark production), whereas leptoquark masses up to 41 TeV may in principle explain the anomalies.


Introduction
Perhaps the most convincing anomalies 1 observed in the LHC data thus far are those seen in ratios of branching ratios of semi-leptonic B-to -K or -K * decays in LHCb [1,2]. Though they involve sensitive measurements of rare processes, they are theoretically clean [3] and apparently a clear signal of violation of lepton universality, a principle that is sacrosanct in the gauge interactions of the Standard Model (SM). Moreover, the fact that such processes arise only at loop level in the SM means that, even though the observed deviations are large compared to the SM contribution, they could plausibly be explained by tree-level exchange of new particles at the TeV scale, with couplings of comparable size to those present in the SM. 2 To put the measurement of these ratios in context, we summarise some of the related anomalies that preceded them: the first sign of a discrepancy appeared in the P 5 observable [11] of angular distributions in B → K * µ + µ − decays [12][13][14][15][16], designed in such a way that hadronic uncertainties cancel out and are under control. LHCb found a 3.4σ anomaly [14], supported somewhat at the 2σ level by a later BELLE measurement [17]. These were also consistent with a 3.2σ tension in B s → φµ + µ − [18]. Indeed, various global fits including LHCb, Belle, BaBar, CMS, and ATLAS data to a variety of b → sµ + µ − kinematic observables indicated a non-zero value for a particular Wilson coefficient parameterising new physics coupling to left-handed quarks and muons, with a statistical pull 4σ [19][20][21][22][23][24][25][26][27][28]. However, these observables could still have been heavily affected by residual theoretical uncertainties in the SM prediction. It was therefore notable that subsequent measurements of the more theoretically clean ratios R K = BR(B → Kµ + µ − )/BR(B → Ke + e − ) and R K * = BR(B → K * µ + µ − )/BR(B → K * e + e − ) both observed deviations at around the 2.5σ level each [1,2]. Moreover, fits to these two clean observables alone demonstrate a pull away from the SM at more than around 4σ on the same Wilson coefficient as the one from the global fit to the other observables [23][24][25][26][27][28][29]. This non-trivial consistency of the various anomalies goes some way towards explaining the level of interest in them, despite the significance of each individual measurement being low.
Even if the anomalies really are signatures of physics beyond the SM (and further data or a better understanding of the SM predictions may well indicate that they are not), we face the problem that the effects we see in B decays arise indirectly, via exchange of virtual states that are far from being on-mass-shell. To confirm the presence of new physics, and to begin the long, but tremendously exciting, programme of exploring Nature's next layer, we will need to produce the new particles directly on-shell, at a current or future collider. But in trying to plan for this, we must overcome a serious obstacle: the size of the effects being seen currently fixes neither the identity, nor the mass, nor the couplings of the new particles. So, at least without further consideration, not only do we not know what energy threshold a collider would need to reach to produce the new states, but also, even if we did know what energy were needed, we do not know what sort of detector, triggering, or cuts might be needed to make the discovery, nor which backgrounds we should strive to better control, nor how much luminosity might be required, and so on.
At least naïvely, we can make some progress on these issues by appealing to the arguments of perturbative unitarity: we know that the loop expansion of quantum field theory, and hence its predictability, breaks down when couplings approach values of 4π or so, and imposing this as an upper bound imposes an upper limit of O(100) TeV or so on the possible masses of new particles [24,30]. A more refined analysis of partial wave unitarity shows that the scale of unitarity violation is actually ∼ 80 TeV [30], and can be even lower in more specific model-dependent cases. Such unitarity arguments successfully predicted the appearance of a Higgs boson at the LHC [31,32], but in the case of the physics inferred from b−decays, the cut-off scale is too high to form a similar no-lose theorem for the next generation of colliders.
Here, we attempt to carry out a rather more detailed analysis of the prospects for discovery of the new physics underlying the b → sµ + µ − anomalies, at current colliders or at proposed future facilities while making as few assumptions about the models as possible. Thus, in the most pessimistic possible scenario, we shall not assume universal couplings to different generations or minimal flavour violating couplings; we only include the minimal new physics that explains the b → sµ + µ − anomalies whilst refraining from adding more model-specific structure (that would typically only lower the scale of new physics or make it more easily discoverable). It turns out that (at least if one is prepared to accept a few simple assumptions along the way) one can make rather detailed and quantitative statements. This is possible for a variety of reasons, which we now describe in turn.
One reason is that, on the theory side, the possible underlying new physics models are rather limited, at least if one assumes that the new physics results in an effective lowenergy operator coupling a left-handed quark current to a left-handed leptonic current.
This assumption is reasonable not only because doing so results in a very good fit to the data (in fact, the best fit to the data, as discussed above), but also because it is highly plausible theoretically, given that the basic objects in the SM are the distinct left-and righthanded fermion multiplets. With this operator, one is limited at tree-level to models with either vector leptoquarks (LQs), scalar LQs or models with new neutral vector particles (Z s) coupling to left-handed currents.
Within this limited range of possible models, there is still a great deal of room to manoeuvre in terms of choosing couplings. But again we can make headway by adopting a conservative approach, leading to predictions that are as pessimistic as possible. For LQ models, for example, it is perfectly consistent mathematically (although highly unlikely in practice), that the LQ has only the Yukawa couplings needed to explain the anomalies, namely to the left-handed lepton doublet containing the muon mass eigenstate 3 and to the two left-handed quark doublets containing the b-and s-quark mass eigenstates. The presence of any other Yukawa couplings (especially those to electrons or light quarks) is likely only to increase the discoverability of the LQ, by providing additional channels for production at a hadron colliders and additional final states that are relatively easy to observe. 4 For Z models, things are a little more complicated, because it is not possible to switch on a coupling to b and s quarks alone: any assignment of charges under the corresponding U (1) gauge symmetry to the three quark doublets in the electroweak basis will lead to other couplings being present in the mass basis. 5 So we consider two different conservative models featuring Z states. In the first model, we allow only a bs coupling in the mass basis. Though mathematically inconsistent, strictly speaking, no inconsistencies arise in the collider phenomenology that we consider here. In the second model, we assume that there is only a coupling to a single generation of quark and lepton doublets (which are those that are mostly b and µ, respectively) and assume that all of the CKM rotation takes place in the down quark sector. Again, for both models our expectation is that any couplings that are additionally present are likely to increase discoverability.
There are also reasons on the experimental side for why a more detailed analysis of the prospects for discovery at a current or future collider is possible. Most importantly, we can extrapolate based on the performance of current colliders, making the conservative assumption that the detector performance will remain roughly the same. This extrapolation is simplified by the fact that the discovery potential of a given machine is largely fixed by our understanding of the backgrounds. In the particular case of searching for a narrow resonance in a given channel at a given centre of mass energy, for example, what is needed is an understanding of the different background contributions (and their uncertainties) in that channel at that centre-of-mass (CM) energy. These backgrounds come, of course, 3 The anomalies in the ratios RK and RK * could, a priori be due to physics in either muonic or electronic operators. But the presence of additional anomalies in purely muonic processes [12][13][14]18], together with the difficulty of accommodating large deviations in flavour physics processes involving electrons, both lead us to assume that the new physics states couple to muons rather than electrons. 4 There is a danger, e.g. by adding charm/tau couplings, of diluting the LQ decays to clean final states, but to study this fully would require an analysis at a level of detail that seems overly premature. 5 Consistency of the theory also requires additional particles for the U (1) gauge symmetry to be anomalyfree [33].
from a combination of the underlying SM physics, which we understand well, together with its manifestation in the detector, which we assume remains similar to current detector performance at the extrapolated energies. This extrapolation is further helped by the fact that the SM is essentially scaleless at the multi-TeV energies that we consider, so that extrapolation amounts to a simple re-scaling, using a procedure outlined and validated in Ref. [34]: in a nutshell, the idea is that the equivalent CM energy at a future collider that gives the same number of background events as a given CM energy in a current search will also yield the same upper limit on a putative signal cross section at that equivalent CM energy.
Proceeding in this way, we are able to obtain a number of simple results, that we believe to be robust within our reasonable assumptions. We find that a 33 TeV high energy upgrade to the LHC 6 should be able to cover most of the Z parameter space that is under perturbative control in our first model with only bs couplings, while it can cover all of the parameter space for our second model with CKM-induced couplings to the first two generations of down-type quarks. A 100 TeV hadron collider has complete coverage for both models; it can therefore discover or exclude any perturbative Z explanation of the anomalies (where the Z width does not exceed 10% of its mass). On the LQ side, considering only pair production via QCD interactions, we find that masses up to 12 TeV can be ruled out in the scalar case. Limits from single production are more modeldependent but become important for O(1) couplings, with sensitivity to LQ masses up to 21 TeV for coupling values up to 4π.
All of this assumes, of course, that the anomalies currently observed are really due to new physics. If it turns out that they are not, the exercise that we have carried out becomes much more academic. But even so, we think that it gives a useful illustration of the complementarity between indirect and direct searches and how one can use anomalies that may plausibly arise in the future, wherever they might occur, to build a concrete strategy for future colliders and particle physics in general 7 .
The paper is organised as follows: in Section 2 we summarise the effective field theory description of the possible new physics parameterising the anomalies, justifying our choice of operator, then describing the possible models that may explain the discrepancy with the SM. In Section 3 we describe the extrapolation method that we adopt for our study, and present our results. We conclude with a summary and outlook in Section 4.

Effective field theory description
Processes involving b → sl + l − transitions can be described by a low-energy effective Lagrangian below the weak scale with the W ± boson, Z boson, Higgs boson and top quark 6 Studies to date have assumed a 33 TeV centre of mass energy, which we choose as a benchmark, but in the future we shall also consider the reduced energy of 27 TeV that can be attained using the 16 T beam magnets currently being designed for FCC-hh. 7 For some reviews of physics at a 100 TeV hadron collider, see for example Refs. [35][36][37][38][39]. The indirect sensitivity of future lepton colliders has been explored in e.g. Refs. [40][41][42][43][44][45].
integrated out. 8 The relevant indirect effects of new physics (and SM weak interactions) are encapsulated by the following four-fermion operators, 9 In the second line we defined dimensionless Wilson coefficientsc l ij normalised by a conventional factor involving elements of the CKM matrix V and ratio of the EFT cut-off scale Λ to the weak scale v 174 GeV such that If new particles with couplings to leptons and quarks of size g NP are integrated out at tree-level, then c l ij ∼ O(g 2 NP ) and since, according to our criterion, the limit of validity of perturbative unitarity is reached when g NP ∼ 4π, this sets an approximate upper limit on the cut-off scale 10 For example, with |c l ij | 1.33, as found in certain best fit values [28], we have Λ max 390 TeV. A more detailed analysis of partial wave unitarity yields a 80 TeV bound [30]. However, other experimental and theoretical bounds will lead to a more restrictive upper limit on the scale of new physics, as we discuss below.
Many global fits to the flavour anomalies have been performed e.g. [19][20][21][22][23][24][25][26][27][28]. Ref. [28], for example, finds that an individual fit to one operator at a time in the muonic sector favoursc µ LL −1.33 at > 4σ significance. A similar conclusion holds for a global fit allowing several operators to vary simultaneously, which then allows an additional subdominant contribution fromc µ LR (thoughc µ LR alone cannot explain the anomalies since it predicts the pattern R K * > 1 when R K < 1 or vice versa). The coefficientsc µ RR andc µ RL , whose contributions must be large to have an effect since their SM interference terms are suppressed, are disfavoured by the relative directions of their pulls on R K and R K * .
In individual fits to R K ( * ) for electronic operators, the anomalies are also well described by eitherc e LL ,c e LR , orc e RR (though the latter two require larger coefficient values due to their suppressed SM interference). Nevertheless the significance decreases substantially 8 If the new physics responsible for the B anomalies is not at low energies [46][47][48][49][50] then the low-energy effective theory can be matched to the SM effective field theory (EFT) [29,51,52]. 9 The relation to coefficients of the O9,10 operators in another commonly used basis is given by c9,10 = ±(cLL ± cLR)/2 [28]. 10 The perturbativity condition is sometimes also taken to be g 2 NP ∼ 4π [24], in which case the cut-off is Figure 1. Feynman diagrams of the two tree-level possibilities for mediating an effective operator that explains discrepancies in B → K ( * ) µ + µ − decays as compared to SM predictions. The diagram on the left hand side shows mediation by a scalar, whereas the right-hand side shows mediation by a flavour dependent Z .
in a global fit including other observables, which shows a clear preference for non SM contributions in decays to muons rather than in decays to electrons 11 . We shall therefore assume new physics to reside solely in the muonic sector and inc µ LL in particular. This restricts the type of heavy particles that can be integrated out to givec µ LL in the EFT, as we discuss next.

Z and LQ models to explain the discrepancy
At tree level there are only a few candidates to consider for mediating the interactions responsible for the B anomalies. These are so-called LQs, that can be either scalar or vector, and Z vector bosons. We shall assume that in each scenario, the new fields are unique representations of the Lorentz group and the SM, i.e. we are not considering multiple identical fields. Feynman diagrams for the relevant interactions are shown in Fig. 1. When the mass of the LQ or Z is much larger than the mass of the decaying B meson, matching to the effective field theory in Eq. 2.1 should provide an accurate approximation to order m B /Λ, where Λ is the mass of the LQ or Z .
Other explanations for the anomalies arise at the loop level. In this case, in order to explain the required size of the non-standard contributions to B → K ( * ) µ + µ − decays, the new particles mediating the interaction must be relatively light and so are more easily discoverable; we therefore restrict our attention to the more conservative case of heavier tree-level induced new physics.
The preference of fits for the O µ LL operator picks out particular combinations of quantum numbers allowed for the LQs [23,25,28] . For the scalar case this is the triplet LQ S 3 , with quantum numbers (3, 3, whose Yukawa couplings to the quark and lepton doublets Q and L are of the form The term proportional to y q induces proton decay and is typically set to zero by imposing baryon number conservation. For the vector case, the O LL operator may be generated by integrating out a singlet V 1 or a triplet V 3 with quantum numbers (3, 1, 2 3 ) and (3, 3, 2 3 ), respectively. The possible couplings are We focus on the couplings generating our operator of interest, O µ LL . Integrating out the LQs with mass M and coupling y gives the Wilson coefficient [25] where κ = 1, −1, −1 and y = y 3 , y 1 , y 3 for S 3 , V 1 , V 3 , respectively. For Z vector bosons, the minimal Lagrangian containing the couplings responsible for generating O µ LL at low energy is given by [28,55]. (2.9) Couplings to some other SM fermions are required by SU (2) L invariance and some additional couplings to other flavours of quark are necessarily generated by CKM rotations when going from the weak to the mass eigenbasis. However, given that these additional interactions are more model-dependent than the ones we write above, we shall take the Lagrangian of Eq. 2.8 as our minimal model (which we call the naïve Z model). Although strictly, the model is incomplete without the additional couplings, the naïve Z model is the most conservative possible case to study; additional couplings will only raise the Z production cross-section, by including couplings to the first two quark generations, and increase the total decay width which is in tension with other constraints. Hence, if a future collider covers some portion of the viable parameter space of the naïve Z model, then we know that a more realistic and complete model will also be covered there (and then some).
To illustrate the size of such effects in a more complete model we shall also consider the case where the Z couples only to third generation left-handed quarks and left-handed muons and neutrinos in the weak basis. The couplings to the first two generations of quarks then arise from CKM rotations, which we assume to be entirely in the down sector. Additionally, if we assume that in the weak eigenbasis all left-handed lepton mixing resides in the neutrino sector, we have a logically consistent model which contains only a coupling to left-handed muons and some family mixture of neutrinos. The precise family mixture of neutrinos is immaterial for collider experiments, since each neutrino is essentially massless and leaves an identical missing momentum signature in detectors. The relevant interaction terms in the Lagrangian for this '33µµ' model are given by where U denotes the PMNS matrix involved in lepton mixing.
With these LQ and Z models in hand, we now turn to their discovery prospects. Some previous studies have examined the 13 TeV LHC's ability to discover other effects caused by new physics involved in the errant b−decays assuming any mediator is not too heavy. In Ref. [55], LHC bounds on Z models that explain the b−anomalies from di-muon resonances were placed, assuming a universal Z coupling to the first two generations of quarks. Ref. [56] also examined current LHC constraints on LQs and performed naïve rescaling to estimate the sensitivity at higher luminosities in models which explain both the b → sµ + µ − anomalies that we consider and additional ones inferred in b → cτν τ decays 12 . Ref. [57] also examined the LHC's ability to detect scalar LQs of 1 TeV mass of the type that we shall examine. Ref. [58] examines the di-lepton final state for effective field theory operators caused by LQs or Z s. Some sensitivity is found under the assumption of minimal flavour violation for light enough Z s.
In our study we look towards future colliders at higher luminosity and energy. In the next Section we shall estimate the projected limits on the Z and LQ masses in our conservative models by extrapolating from direct searches at the LHC.

Limit extrapolation method
We follow the approach of Ref. [34] to extrapolate the limits on direct searches for new resonances at the LHC to higher energy and luminosity. The method assumes that such a limit is entirely driven by the number of background events, so that finding the equivalent mass at a future collider that gives the same background as a given mass in a current search will also yield the same upper limit on a putative signal cross section at that equivalent mass.
Concretely, the background cross-section at a resonance mass M and centre of mass collision energy √ s is whereσ ij (ŝ) is the partonic cross section for production of the resonance by partons i and j evaluated at a partonic centre of mass energy √ŝ and the parton luminosity function dL ij /dŝ for the initial state parton pair labelled by i and j is given by We set the factorisation scale µ = √ŝ . We assume that the resonance is sufficiently narrow, ∆ŝ M 2 , such that the partonic luminosity is approximately constant in the integration 12 In the present paper, we do not consider physics due to these charged current decays because their SM predictions are subject to larger theoretical errors. Moreover, the size of those effects requires a low mass scale that would make the new physics responsible more easily discoverable than the source of the neutral current b-anomalies.
region. Since the background consists of SM processes at energies far above the weak scale, the partonic cross-section should scale likeσ ij ∝ 1/ŝ. The total background cross-section then simplifies to where C ij =ŝσ ij is approximately constant. The number of background events at a given luminosity L is N B = L · σ B (M, s). If a 95% confidence level (CL) limit on a signal crosssection is set for a given resonance mass M 0 at a present collider (with energy √ s 0 and luminosity L 0 ), then we find the equivalent mass M for which the limit applies at a future collider (with energy √ s and luminosity L ) by the assumption that the same limit is applicable when The fixed relative width ∆ŝ/M 2 and other prefactors have cancelled out, leaving a straightforward equation to solve for M . The constants C ij can be normalised such that they represent the relative weights of the contributions from each parton pair. This method introduces some arbitrariness in the starting point of the extrapolated exclusion curve, since it depends on a re-scaling by the luminosity ratio L 0 /L . If L = L 0 then the smallest mass M 0min at the lower end of the current collider sensitivity will be extrapolated to the starting point M min of the exclusion curve at the future collider. On the other hand if L > L 0 then the starting point will be at a higher mass point, while L < L 0 would reach lower masses. A conservative procedure to account for this artificial effect is to smoothly vary the future collider luminosity up to L during the extrapolation and take the strongest limit for each mass point, which only affects the limit for masses below M min , and in any case is more conservative than a realistic limit [34].
This extrapolation method has been validated against a cut-and-count-based analysis for di-lepton searches in Ref. [34], where agreement is found up to a factor of two for a width of ∆ŝ/M 2 = 10%. Results from the approximate method outlined here can then be trusted in so far as a more complete analysis does not give limits too far off from a cut-and-count-based one. While more realistic experimental analyses will certainly use more refined methods that go beyond our assumptions, the approximation is sufficient for a rough estimate of future collider sensitivity and should help motivate a more detailed study.

Z sensitivity
We extrapolate limits from the ATLAS 13 TeV search in the di-muon final state at √ s = 13 TeV and 3.2 fb −1 [59] 13 . The dominant backgrounds come from Drell-Yan, tt and di-boson production. Using the procedure described in Section 3.1, we obtain the projected limits displayed in Fig. 2. The solid black line in the left plot is the current 95% CL limit from the ATLAS 13 TeV analysis. In dashed black is the projected limit for HL-LHC at 14 TeV with 3 ab −1 , while the solid and dashed lines in cyan are for the HE-LHC at 33 TeV with 1 and 10 ab −1 , respectively. The plot on the right shows the corresponding FCC-hh 100 TeV limits in solid (dashed) red for 1 (10) ab −1 . The shaded regions on the curves indicate the point at which the extrapolation method underestimates the actual limit at low masses, as explained in Section 3.1.
One may note various features in Fig. 2 that might seem surprising prima facie: for example, it appears that the 14 TeV 3 ab −1 HL-LHC can reach lower in σ ×BR than the 10 ab −1 33 TeV HE-LHC for M < 6 TeV. This is caused by the behaviour of regions dominated by high backgrounds at lower masses: if one increases the centre of mass energy from the LHC to higher collider energies then this background-dominated region will correspondingly move to higher masses. On the other end we see that at the highest values of M the HE-LHC is the most sensitive, as expected. While these sensitivity limits are purely a function of the background, the actual limit set for a given Z mass and coupling also depends on the signal cross-section, which is larger at higher collider energies. Therefore a lower-energy collider whose limit curve reaches further down than that of a higher-energy collider does not necessarily translate to better sensitivity in a model's parameter space.
The actual Z mass that can be excluded for a B-anomaly-compatible model depends on the specific couplings of the Z and its total decay width. We calculated the Drell-Yan cross-section for pp → Z → µ + µ − as a function of these couplings using the following expression in the narrow width approximation, cross-checked with MadGraph [60], where S i and C i are the number of spin and colour degrees of freedom of parton i respectively, and the parton luminosity function is The decay rate for Z into fermions with coupling g f , assuming m f M Z , is given by For the parton distribution functions f (x, µ 2 ) we use the 5-flavour NNPDF2.3LO (α s (M Z ) = 0.119) set [61] with LHAPDF [62] and fix the factorisation scale to be M 2 Z . We consider b quarks to be in the initial PDFs of the proton, thus re-summing large logarithms on the initial b quark line [63]. The Feynman diagram for hadron collider production is therefore identical to the right-hand plot of Fig. 1.
Using these expressions and the extrapolated limits of Fig. 2, the resulting parameter space for the naïve model is shown in Fig. 3. As discussed in Section 2.2, we take the naïve model of Eq. 2.8 defined by only a Z coupling tobs +sb and µ + µ − , and nothing else, as the most conservative possible case. While other couplings should necessarily be present, Figure 2. Projected sensitivity of di-muon resonance searches of some future hadron colliders to Z models that may explain anomalous B → K ( * ) µ + µ − decay results for the luminosities and centre of mass energies given in the legend. Shaded parts of the curve indicate the conservative extrapolation method at low masses that underestimates the actual limit. the naïve Z serves as a useful scenario to assess the sensitivity of a future collider since any model that seeks to explain the B-anomalies must have at least these couplings, while other interactions are more model-dependent.
The line and colour coding for Fig. 3 is as follows: the blue-shaded region covering the area vertically towards the right corresponds to the extrapolated 95% CL limit for the highest luminosity at the collider energy shown in the plot title; the grey-shaded region excluding the area horizontally towards the top is where the Z width exceeds 10%; the vertical red region is excluded by too large a contribution to B s −B s mixing which constrains |ḡ sb L | √ 2M Z /(210 TeV) [28]; the green region is compatible with the B anomaly within 1σ of the best fit value of Ref. [28]; the blue (mostly) vertical dashed lines are the contours of cross-section in units of femtobarns; and the horizontal grey dashed line is where [55] , (3.8) indicating that the Z couplings will hit a Landau pole before the Planck scale; the region above this dashed line is therefore theoretically disfavoured. This last condition is model dependent as the Feynman diagram contributing to the decay width, given by the imaginary part of the Z propagator, will also contribute to the renormalisation group running from the real part of the propagator. While this perturbativity condition is weakened by new vector bosons contributing to the running, it is strengthened by the addition of scalars or fermions.
The top left plot in Fig. 3 indicates that the HL-LHC at 14 TeV and 3 ab −1 is barely sensitive to a naïve Z when its mass is 1.5 TeV. This may seem low but we recall that in the naïve model the only production mechanism for Z in Drell-Yan is through b and s initial state partons 14 . On the other extreme end of collider reach is the FCC-hh at 100 14 For a study of other possible production mechanisms with these couplings, see Ref. [64]. TeV, shown on the top right for M Z = 15 TeV. We see that a 15 TeV Z is at the limit of being anomaly-compatible and evading the constraints from both B s −B s mixing and Landau poles. Nevertheless, the blue region corresponding to FCC-hh with 10 ab −1 can easily cover all of the parameter space of interest. For lower luminosities the sensitivity can be read off from the cross-section contours and the corresponding limits in Fig. 2.
Between the CM energies of HL-LHC and FCC-hh is the HE-LHC at 33 TeV CM energy, displayed in the bottom row of Fig. 3 for a Z mass of 2 TeV on the left and 5 TeV on the right. The sensitivity drops off such that the HE-LHC no longer covers any non-excluded parameter space for M Z 7 TeV.
To illustrate the possible sensitivity to a more realistic model, in Fig. 4 we show the reach for the 33µµ Z model defined by the Lagrangian of Eq. 2.10. There, the couplings to third-generation left-handed quarks induce a coupling to the first two generations of quarks through the CKM matrix. This raises the production cross-section through the additional initial state partonic channels, and also increases the total decay width. In the top left-hand plot we see that a 1.5 TeV Z is now accessible to the HL-LHC in all of its favoured parameter space, with the top right-hand plot indicating that the new limit of sensitivity of the HL-LHC for this more realistic model is raised to M Z 4 TeV. From the bottom two plots, with M Z = 5 (10) TeV on the left (right), we conclude that the 33 TeV HE-LHC at its highest luminosity can cover all the parameter space of interest for all favoured masses. Indeed, we see that for M Z 10 TeV the anomaly-compatible region lies entirely within the grey and red areas and yet is still covered by the blue-shaded area. The FCC-hh with even more energy will therefore also be sensitive to the entire mass range, so we omit its plot.
To summarise the projected reach, we now study the behaviour of the bounds and future collider coverage of Z models shown in Figs. 3 and 4 for a continuously varying M Z , shown on the abscissa. We scan along the central green line in those figures, corresponding to the central inferred value ofc µ LL = −1.33 [28], and plot the value of g µµ L − g sb L along this line on the ordinate. We see from the right-hand side plots in Fig. 5 that the naïve model is not covered much at all by di-muon resonance searches at the LHC, even at high luminosity, but that a 100 TeV 10 ab −1 collider can cover all of the viable parameter space where the Z is narrow (we note that the sensitivity at low masses is underestimated by our limit extrapolation technique, as explained in Section 3.1). However, the naïve model is a limiting case that underestimates both the potential sensitivity and the current constraints for a more realistic model. We see in the left-hand plots that in a more complete 33µµ model, a 14 TeV 1 ab −1 LHC can cover a decent portion of the viable parameter space and a 33 TeV LHC collider is sensitive to all of it.  Figure 6. Example Feynman diagrams of LQ production at a hadron collider followed by subsequent decay of each into µj.

LQ sensitivity
There are many dedicated experimental studies of LQs. For some recent examples, CMS have searched for first and second generation LQs in pair production [65][66][67] and single production [68] at 8 TeV centre of mass energy, while ATLAS set limits on the pair production of third generation LQs using 7 TeV data [69] and first and second generation LQs with 13 TeV [70]. A summary of LQ searches by ATLAS and CMS can be found in Ref. [71]. LQs were recently reviewed in Refs. [72,73].
As the basis for our extrapolation, we take the 95 % CL limits from the CMS 8 TeV search for a pair of second generation scalar LQs with 19.6 fb −1 of integrated luminosity [65], focusing on the µµjj channel in particular, as shown in Fig. 9. The current limits exclude masses up to 1070 GeV, assuming a 100% branching fraction into a charged lepton and quark. We note here that pair production proceeds through the strong interaction and so limits coming from the experimental search may be phrased as only depending on the LQ mass, once the assumption about its branching fraction is made.
Following the extrapolation procedure detailed in Section 3.1, we obtain the weighted sum of parton pair luminosities for the dominant contributions to the background processes, in this case Z/γ * +jets and tt, then find the equivalent mass at a future collider that gives the same number of background events. The results for the projected limits are shown in Fig. 7. In the left hand plot, the exclusion curve in solid black is the current CMS 8 TeV exclusion curve, while the dashed black line shows that the LHC reach can be extended for 14 TeV at high luminosity (HL-LHC) with 3 ab −1 . The cyan-coloured limits are for a potential highenergy upgrade to the LHC (HE-LHC) that could reach up to 33 TeV centre of mass energy. The solid and dashed lines represent 1 and 10 ab −1 of integrated luminosities, respectively. It appears that at low masses, the CMS 8 TeV analysis is more sensitive (when phrased in terms of σ × BR) than when the energy is upgraded to 14 TeV at the HL-LHC. This is an artefact of the arbitrariness in the starting point of the extrapolated exclusion curve, as explained in Section 3.1, where below this point lower luminosities can set limits at lower masses, though this conservative procedure underestimates the actual limit. The regions below the extrapolated starting point are shaded on top of their respective curves. On the right-hand side of Fig. 7 we display the limits for a 100 TeV proton-proton future circular collider, the FCC-hh, at 1 (10) ab −1 in solid (dashed) red.
The dotted lines superimposed on both plots are theoretical calculations at next-toleading order for the LQ pair production process, using the code of Ref. [74]. Up to O(1) uncertainties, we see that HL-LHC can exclude LQ masses up to 2 TeV, while HE-LHC Figure 7. Projected sensitivity of future colliders to di-LQ production (where each decays to a muon and a jet) for the luminosities and centre of mass energies given in the legend. We also show the scalar LQ cross-section times branching ratio predicted for some future collider scenarios by the curves labelled σ N LO × BR. Shaded parts of the curve indicate the conservative extrapolation method at low masses that underestimates the actual limit. can roughly double that to 4 (5) TeV with 1 (10) ab −1 . At FCC-hh the limits are improved by an order of magnitude with respect to current searches, reaching exclusions up to 10 and 12 TeV for 1 and 10 ab −1 , respectively.
These projected bounds on the LQ mass are to be compared with the upper limit allowed by B s −B s mixing. The relevant four-fermion operator of the effective Lagrangian for this process can be written as  Figure 9. Example Feynman diagrams of single LQ production at a hadron collider followed by its subsequent decay into µj. Note that the LQ production cross-section depends upon its coupling to fermions, in contrast to the pair production cross-section depicted in Fig. . The Wilson coefficient gets a contribution from the coupling combination |y bµ y * sµ | that is given by [25] c bb LL = k (3.10) where y = y 3 , y 1 , y 3 and k = 5, 4, 20 for the S 3 , V 1 , V 3 LQs, respectively. Using this expression, together with Eq. 2.7 and the experimental limit from B s −B s mixing that constrains c bb LL 1/(210 TeV) 2 [28], we obtain the parameter space shown in Fig. 8. The couplings as a function of mass are displayed on a log-log scale, and the green strip represents the parameter space compatible with the B-anomalies at 1σ. The different shades of red are excluded by B s −B s mixing for the S 3 , V 1 , V 3 LQs up to the solid red, dotted brown, and dashed pink lines, respectively. We see that the maximal values of the LQ masses allowed by B s −B s mixing correspond to M LQ = 37, 41, 18 TeV for S 3 , V 1 , V 3 , respectively. The blue region shows the 95% CL limits for scalar LQs at a 100 TeV collider with 10 ab −1 , such as the FCC-hh. The pair production process for vector LQs is more model-dependent (unlike scalar LQs, whose gluon interactions are fixed by the SU (3) c gauge couplings) but is typically stronger than the scalar case [75][76][77].
The direct search sensitivity may also be extended to heavier LQs by considering single LQ production [78], as in Fig. 9. For large enough couplings the limits may be be stronger than those obtained in pair production [68], but the product of the bµ and sµ couplings must remain within the stringent bounds from B s −B s mixing.
We extrapolate the current limits from a direct search by CMS for a single scalar LQ produced at 8 TeV with 19.6 fb −1 [68]. CMS places a bound of M LQ 660 GeV for a second generation LQ with coupling to sµ of order unity. For our signal cross-section we also include a bµ coupling since we take the b quarks to be in the 5-flavour initial parton distribution function NNPDF2.3LO (α s (M Z ) = 0.119) [61]. This re-sums the large logarithms of the initial state b-quark line. We integrate the partonic cross-section with the parton distribution functions using LHAPDF [62]. The partonic cross-section at leading order for a scalar LQ φ is given by [79] σ(qg → φl) = y 2 α S 96ŝ 1 + 6r − 7r 2 + 4r(r + 1) ln r , (3.11) where r = M 2 LQ /ŝ and we set y sµ = y bµ = y for simplicity. This expression has been crosschecked with Fig. 8b of Ref. [68] and found to agree within partonic uncertainties. The extrapolated limits and production cross-sections for a coupling and branching ratio set to Figure 10. Projected sensitivity of future colliders to single LQ production that decays to a muon and a jet, for the luminosities and centre of mass energies given in the legend. We also show the cross-section times branching ratio for some future collider scenarios by the curves labelled σ y ×BR, where y is the scalar LQ coupling to bµ, set equal to the coupling to sµ. Shaded parts of the curve indicate the conservative extrapolation method at low masses that underestimates the actual limit.
1 are displayed in Fig. 10, with the same colour coding as Fig. 7. The signal cross-sections at 14 and 33 TeV are shown as dotted lines for y = 1 in black and cyan respectively, on the left plot. On the right we have the signal cross-section for 100 TeV with y = 1, 2, and 12 in red, black, and green dotted lines respectively. We see that for couplings y = 1, the limits are comparable to pair production but slightly lower. On the other hand for y = 2 the limits at 100 TeV go up to 15 TeV for 10 ab −1 , extending to 21 TeV for y ∼ 4π. The reach can be further extended for a model with additional quark couplings. Note however that in a realistic model the product of y bµ and y sµ must still be anomaly-compatible within the B s −B s mixing bounds shown in Fig. 8, so that these limits only apply when one coupling is taken large with the other small.

Conclusion
Some measurements of B → K ( * ) l + l − decays disagree with SM predictions: using only theoretically clean quantities, the discrepancy on a Wilson coefficient with respect to the SM value is at around the 4σ level [28]. More specifically, the ratio of decays to muon pairs and electron pairs is predicted to be 1.0 in the SM, but is measured to be lower than this value both for K and K * in the final state, each in two different bins of transferred 4-momentum (squared). Moreover, the 4σ pull from these clean observables on a Wilson coefficient parameterising new physics is not only statistically significant, but also in the same direction as another independent 4σ pull due to other (less clean) observables. The combined significance in a global fit is then significantly larger than 4σ.
Many authors have constructed bottom-up models containing new particles in order to change the apparent predictions and explain the discrepancies. In particular, it appears that lepton flavour universality should be broken by the new particles, which should have chiral interactions. At tree level, there are only two classes of new particle which explain the discrepancies: flavourful Z s and LQs. We choose these two cases to examine future hadron collider sensitivities: there are other possibilities from particles which affect the decays at the loop level, but because of the loop suppression, these particles should be a factor of roughly 4π lighter than the tree-level cases, and should therefore be easier to detect directly by production in a collider. Studying the tree-level possibilities is then conservative: if one shows that one can discover these, it should be easier to find the lighter particles that are predicted by the loop effects.
We found that for the Z models, a 100 TeV future collider will essentially cover all of the parameter space that can explain the B → K ( * ) l + l − decay data in a resonant dimuon search. Examining more complete models than the naïve Z model such as the 33µµ model, we see that even a 33 TeV run of the LHC may cover all of the relevant viable parameter space. The more complete models contain more model dependence, but have stronger bounds and more coverage than the naïve model. One caveat to our analysis is that we only consider a narrow Z with width less than a tenth or so of the mass. This will not necessarily be the case, but wide Z s invalidate the procedure we use to extrapolate current bounds from scaling the background detailed in Section 3.1, and so require a more detailed simulation of backgrounds and signal at high energies. Nevertheless, a wide Z is theoretically disfavoured by the large couplings required: they run into Landau poles. The Z case is summarised in Fig. 5, where expected sensitivities, bounds and the validity limit of our analysis are plotted for various different future hadron collider assumptions.
Coverage of the relevant LQ models is significant, but somewhat less complete than the Z models: whilst LQs of varying kinds up to masses of 41 TeV may explain the B → K ( * ) µ + µ − decay data, we show that the expected sensitivity of di-LQ production into a µµjj final state reaches up to 12 TeV for a scalar LQ. In more model-dependent cases the reach can be higher, as for example in pair production for vector LQs and single production for both vector and scalar LQs, which depend on a choice of couplings. We estimated the sensitivity of single scalar LQ production and found that O(1) couplings can reach a sensitivity up to LQ masses of around 21 TeV at the strong coupling limit. Whilst our extrapolation of current LHC limits is rather rough (one may expect an uncertainty of a factor of two in the cross-section times branching ratio for the limit due to PDF uncertainties and different detector effects etc.), we estimate that this only results in an uncertainty on the sensitivity of around ±1 TeV when expressed in terms of the mass of leptoquarks or Z particles (as evidenced by the steep model prediction curves in Figs. 7 and 10, for example).
A potential loop-hole in our analysis would occur if one assumed the existence of multiple (N ) mediators of the c µ LL operator. Under the assumption that each of the N mediators (Z s or LQs) has an equivalent mass and identical couplings, one obtains a contribution to c µ LL that is proportional to either N |y| 2 /M 2 in the LQ case, or N g sb L g µµ L /M 2 Z in the Z case. For an identical effect on the measured b decays as in the unique mediator case, each of the N mediators could therefore be heavier by a factor √ N or more weakly coupled. The LHC mediator production cross-section falls with a power of the mass that is significantly higher than two, resulting in weaker collider sensitivity despite a factor of N from the production of more new mediators.
Of course, there is always the possibility that the current discrepancy with SM predictions is due to a fluke. In this case, our paper still serves a purpose, estimating the reach of future colliders into particular flavourful Z or LQ models and demonstrating the interplay between indirect and direct searches for new physics. In any case, new empirical data on the B → K ( * ) l + l − decays are expected from Belle II and LHCb in the next few years. Ref. [80] points out that by 2020, the number of bb pairs produced inside the acceptance of LHCb should increase by a factor of 3.7 as compared to those produced before and during 2012. For example, this would result in a discovery of a non-SM effect beyond the 5σ level in R K , and close to a 5σ level effect in R K * from LHCb data alone if the central values were not to change from their current values [80].
We have shown that there is significant coverage of all beyond the SM explanations of the current anomalies in proposed future hadron colliders. Thus, if the signal significance of non-SM effects in B → K ( * ) µ + µ − decays increases, so does this particular motivation for higher energy future colliders 15 .