Sensitivity of Future Hadron Colliders to Leptoquark Pair Production in the Di-Muon Di-Jets Channel

We estimate the future sensitivity of the high luminosity (HL-) and high energy (HE-) modes of the Large Hadron Collider (LHC) and of a 100 TeV future circular collider (FCC-hh) to leptoquark (LQ) pair production in the muon-plus-jet decay mode of each LQ. Such LQs are motivated by the fact that they provide an explanation for the neu-tral current B -anomalies. For each future collider, Standard Model (SM) backgrounds and detector effects are simulated. From these, sensitivities of each collider are found. Our mea-sures of sensitivity are based upon a Run II ATLAS search, which we also use for validation. We illustrate with a narrow scalar (‘ S 3 ’) LQ and ﬁnd that, in our channel, the HL-LHC has exclusion sensitivity to LQ masses up to 1.8 TeV, the HE-LHC up to 4.8 TeV and the FCC-hh up to 13.5 TeV.


Introduction
There has been much attention recently on various measurements in rare B-meson decays, since several of them are in apparent disagreement with Standard Model (SM) predictions. The measurements that we focus on here involve processes with muon pairs, bottom quarks and strange quarks. Most of the disagreements are at the 2 − 3σ level and so do not warrant particular consternation of and by themselves. However, if one takes them collectively, it seems as if a pattern may be emerging. Of particular interest are the ratios of branching ratios both predicted by the SM to be 1.00 with high precision in the di-lepton invariant mass squared bin m 2 ll ∈ [1.1, 6] GeV 2 . LHCb measurements [1,2] imply R K = 0.846 +0.060 −0.054 +0.016 −0.014 a e-mail: mum20@cam.ac.uk (corresponding author) and R K * = 0.69 +0.11 −0.07 ± 0.05 in this bin. Other observables disagree with their SM prediction, despite larger theoretical uncertainties. For example, the branching ratio of B s → μ + μ − [3][4][5][6] is also measured to be lower than the SM prediction. Discrepancies with SM predictions [7,8] include some of the angular distributions in B → K ( * ) μ + μ − decays [9][10][11][12]. Collectively, we refer to the discrepancies between these measurements and SM predictions as the neutral current B-anomalies (NCBAs). 1 We shall take the hypothesis that the NCBAs are harbingers of new particles with flavour dependent interactions. We suppose that the new particles are much heavier than B mesons, such that they will not have necessarily already been produced and discovered in previous high energy colliders. In SM effective theory, global fits find that the two most relevant beyond the SM (BSM) Lagrangian operators for describing such new particles are [13][14][15][16][17][18] The (1/36) 2 TeV normalisation makes the C L ,R dimensionless, and for C L ,R ∼ O(1) they can fit the NCBA data. Bearing Eq. 2 in mind, the only tree-level solutions explaining these deviations are those of Z s and leptoquarks. Pioneering projections for the direct discovery of such new particles at the LHC and future colliders, in a simplified framework, were made in Refs. [19,20]. More accurate estimates involving simulations of collisions and detector response were made in Refs. [21,22] for the Z case. Heavy LQs with couplings to b quarks, s quarks and muons predict the BSM operators in Eq. 2 [23][24][25][26][27][28][29][30]. Here we shall consider the easier example of a scalar LQ. Vector leptoquarks can also successfully explain the NCBAs, but to obtain sensible high energy behaviour, require a full ultraviolet model to be devised, introducing further model dependence.
A complex SU (2) L triplet scalar S 3 , with quantum numbers (3, 3, 1 3 ) under the standard model gauge group SU (3) c × SU (2) L × U (1) Y , can produce the BSM operator with coefficient C L in Eq. 2. In fact, it is the only single scalar LQ progenitor. There are two alternative vector boson LQ progenitors: U 1 and U 3 with quantum numbers (3, 1, 2 3 ) and (3, 3, 2 3 ), respectively. Many of our results (for example, bounds on production cross section times branching ratio) shall apply equally to other LQs such as U 1 and U 3 . However, we shall illustrate some other estimates that depend on calculation of signal cross section and branching ratio (for example, sensitivity to m LQ ) solely for the S 3 case, vector LQ simulations being beyond the scope of the present paper, since they add model dependence concomitant with the necessity of providing a more complete ultraviolet model. 2 The Yukawa couplings of the LQ to the i th -family SM quark (Q i ) and lepton (L i ) SU (2) L doublets are given by (in the primed weak eigenbasis) [31,32]: where we have suppressed QCD gauge indices, i, j ∈ {1, 2, 3} are family indices (repeated indices have an implicit summation convention), a, b, c ∈ {1, 2} are fundamental SU (2) L indices, k ∈ {1, 2, 3} is an adjoint SU (2) L index, the superscript C denotes a charge conjugated fermion, ab is the Levi Civita symbol, the τ k ab are the Pauli matrices and Y L and Y Q are 3 by 3 matrices of complex dimensionless Yukawa couplings.
In order to avoid proton instability we assume that baryon number is conserved, setting (Y Q ) i j to zero in consequence. After electroweak symmetry breaking, S 3 becomes (S −2/3 , S +1/3 , S +4/3 ) where the superscript denotes electric charge. The left-handed quarks and leptons mix according to where 3 P ∈ {u L , d L , e L , ν L }, bold face denotes a 3vector in family space and the unprimed basis is the mass eigenbasis. V P are then unitary dimensionless 3 by 3 matrices, being experimentally constrained via the Cabbibo-Kobayashi-Maskawa combination 2 Vector leptoquarks' couplings to gluons depend on a free parameter in contrast to scalar leptoquarks whose interactions are fixed by the QCD coupling. 3 Here, the transpose denotes a column vector.
In the mass eigenbasis, where In order to describe the NCBAs, we require (Y de ) 32 = 0 and (Y de ) 22 = 0. Then, Fig. 1 demonstrates how S 3 contributes to the NCBAs via tree-level exchange. For LQ masses much larger than B meson masses (m LQ m B ), the effective field theory matches C L = −1.06 ± 0.16 fits combined NCBA data [16]. We shall typically use the central value from the fit in order to fix (Y de ) 32 (Y * de ) 22 for a given value of m LQ . Although in general Y de are complex, we shall here take real values for simplicity and because we are not considering C P-violating observables.
Our LQs contribute to B s −B s mixing at one-loop order, as shown in Fig. 2. A recent determination of the SM prediction for B s −B s mixing is broadly in agreement with the experimental measurement and so upper bounds can be placed on the LQ contribution. However, a recent determination [33] finds that this is not very constraining for S 3 leptoquarks that fit the NCBAs. Perturbative unitarity provides the stronger constraint that m LQ < 68 TeV.
At hadron colliders, the dominant mechanism for pair production of the LQ is by gluon fusion, as shown in Fig. 3. Since by definition, a LQ couples to a quark and a lepton, it Fig. 3 Example Feynman diagrams of LQ pair production at a hadron collider followed by subsequent decay of each LQ into μj must carry colour to preserve SU (3) and therefore, by SU (3) gauge symmetry, must couple to gluons via the QCD coupling constant. Run II of 13 TeV running at the LHC has produced some 139 fb −1 of integrated luminosity each for ATLAS and CMS. A search for NCBA-solving LQs using all of this data is eagerly awaited, having not appeared yet. The next phase of LHC running will be in the 14 TeV high-luminosity phase (HL-LHC), with a design integrated luminosity of 3 ab −1 . This phase will provide much new information on the NCBAs [34] concurrently with Belle II [35,36]. At the same time, direct searches for new particles [37] may include S 3 . One potential future LHC upgrade would be to insert 16 Tesla magnets into the current LHC ring, resulting in the high energy LHC (HE-LHC), running at a nominal energy of 27 TeV with a design integrated luminosity of 15 ab −1 [38]. Ultimately, such magnets could be placed within a much larger tunnel, resulting in the Future Circular Collider, which could collide protons (FCC-hh) at a centre of mass energy √ s = 100 TeV with a design luminosity of 20 ab −1 [39,40]. We arrive at the question central to this paper, which is:

For LQs which fit the NCBAs, what is the LQ mass sensitivity of future hadron colliders?
We wish to estimate the sensitivity for the Run II, HL-LHC, HE-LHC and FCC-hh options. Table 1 summarises the centre of mass energies and integrated luminosities that will be used for each collider in our estimates. We hope that this will help inform the European Strategy for particle physics, which is currently deliberating on various scientific priorities.
A previous estimate of future collider sensitivity to S 3 LQs consistent with the NCBAs was made in Ref. [19], which projected current sensitivity to higher centre of mass energies and luminosities. However, the sensitivity estimate had two highly dubious approximations. The first was that experimental efficiency and acceptance did not change with centre of mass energy. In fact, at large m LQ and at high energies (particularly at FCC-hh), the decay products from LQs will be highly boosted. This has two effects: the muons will be pushed closer to the jets, meaning that more of them will fail isolation criteria. Also, at higher energies, the muon momentum resolution is likely to be very poor, since such hard muons can only be bent to a limited extent by the magnets. This will also affect the signal efficiency from peak broadening. The second dubious approximation was that the LQs are produced exactly at threshold. This is likely to introduce large uncertainties. We shall rectify these approximations in our paper by performing a fast simulation of the signal and backgrounds, as well as including detector response. The first of these approximations has already been found to have non-trivial effects upon the predicted future hadron collider sensitivity of Z explanations of the NCBAs [22,41]. The estimate in this paper should be much more accurate than the previous pioneering determination in Ref. [19].
Searches for LQ pair production with subsequent decays of each into a muon and a jet have already been performed at the 13 TeV LHC. The ATLAS Collaboration set a 95% confidence level lower limit on m LQ of 1.05 TeV from 3.2 fb −1 of pp collisions [42]. This is a simple cut-based analysis, which we adopt for estimating future hadron collider sensitivity. More recent experimental analyses were made more sophisticated in order to squeeze more sensitivity out of them. The CMS Collaboration maximise their sensitivity using a multi-dimensional optimisation of the final selection for each m LQ in 36 fb −1 of delivered beam at the LHC [43], finding a 95% CL lower bound of m LQ > 1.28 TeV. The ATLAS collaboration has also performed a search in 36 fb −1 of 13 TeV pp collisions for LQs decaying to muons and jets. They utilise differential cross-section measurements and boosted decision trees to obtain a lower bound of m LQ > 1.23 TeV. However, such a level of sophistication is unnecessary for our purposes, where the uncertainties involved in estimating future collider sensitivities (for example because we do not yet know the experimental design) are much larger than the gain in sensitivity. Thus, following the much simpler methodology in Ref. [42] is sufficient for our purposes. The NCBAs predict that there should be couplings between S 3 andb, μ from the first term in Eq. 5. Thus we expect a decay channel S 3 →bμ to be open. In the experimental analysis we choose, the bottom quark remains untagged and is counted merely as a light jet. We note that the second term in Eq. 5 may simultaneously predict S 3 decays to top quarks and muons. This mode is more complicated than the one we choose for analysis, and we leave it for future work. For a discussion of potential analysis strategies for this decay mode, see Ref. [44].
Collider sensitivity to LQ pair production is limited by SM background rates. Therefore the estimation of such back-ground rates are of vital importance to the estimate of the sensitivity to LQ pair production. The paper proceeds as follows: in Sect. 2 we describe the SM backgrounds, how they are simulated. We validate our estimate of the backgrounds and resulting LQ limits against the ATLAS determination for √ s = 13 TeV. Then we estimate the backgrounds at future hadron colliders. Next in Sect. 3 we present the sensitivity estimates for future hadron colliders, before summarising in Sect. 4.

Standard model backgrounds
Consider the pair production of LQs and their decay to a μμj j final state. Following previous searches for leptoquark pair production, we define the parameter m min (μ, j) from the kinematics of these four final state particles by finding the configuration of muon-jet pairings which minimises the difference in invariant masses |m(μ 1 , where j 1 and j 2 are the hardest two jets in an event. In an on-shell LQ pair production event this parameter will approximate the LQ mass m LQ . To estimate the sensitivity of future colliders to LQ pair production in the μμj j channel we simulate the distribution of the SM background in m min (μ, j). We select events containing exactly two muons with no charge requirement and at least two jets with no flavour requirement. Our background simulations for 13 TeV are validated against the results presented in [42]. More details can be found in Sect. 2.2. We place limits on σ ×BR and determine the maximum LQ mass m LQ which could be excluded at 95% CL by each collider, assuming the observed data is consistent with the SM background. Alternatively, assuming a LQ exists at mass m LQ , we estimate the discovery potential by finding the integrated luminosity required for a 5σ significance.

Methodology
We generate the parton level SM background events at leading order in Madgraph5 [45]. These events are then passed to Pythia8 [46] for the simulation of initial state radiation, parton showering and hadronisation. Finally, the hadronlevel events are passed to Delphes3 [47] for detector simulation. We use the 5-flavour NNPDF2.3LO [48] parton distribution function via LHAPDF6 [49] for all background simulations except for di-boson production, for which the 4flavour NNPDF2.3LO parton distribution function is used. This choice is made to remove interference in di-boson production, as outlined in more detail later in this section.
There are four significant contributions to the SM background in the μμj j channel. These are Drell-Yan (Z /γ * → μ + μ − ), top pair production (tt), single top production in 6 Single top production association with a W boson (W t) and di-boson production (W + W − ), where top quarks decay leptonically to muons. An example of the production of each component of the background is shown in Figs. 4, 5, 6 and 7. Other sources of background include misidentified muons from W +jets, single top production in the s and t channel or multi-jet events. These form a negligible component of the background in comparison and therefore we treat Drell-Yan, top pair production, single top and di-boson production as the only sources of background.
To contribute to the μμj j signature, Drell-Yan and diboson production require the addition of at least two jets from initial and final state QCD radiation. Similarly at least one extra jet must be added to single top production. To account Di-boson production for this we generate events from processes of a range of different jet multiplicities according to the following definitions: We include processes with less than two final state jets at parton level to account for the possibility that sufficiently hard jets may be produced by the parton shower algorithm. We use MLM matching [50] to match the final state partons generated from matrix elements in Madgraph5 to those produced by parton showering in Pythia8. This removes overcounting between multi-jet final states and accounts for the fact that while soft and collinear jets are well described by the parton shower, the matrix elements are more suited to simulating hard and well-separated partons. MLM matching is implemented in Madgraph5 by specifying a nonzero value of the jet cut-off xqcut to be approximately 1/3 times a hard scale in the process for each component of the SM background and for each collider. We confirm our choice of xqcut value in each case by checking that the differential jet rate distributions are smooth and that observables such as the total cross section are insensitive to changes in xqcut about the chosen value.
Interference arises between W W + 2 jets, W t + 1 jet and tt production: all three processes may produce a W W bb final state via tt production. Ideally we would simulate all contributions to the W W bb final state simultaneously, but this was found to be very computationally difficult and impractical. Instead, we generate events from each process separately and then combine them to produce the SM background. This means we must define each process in our simulations in such a way that any overcounting is removed.
A number of methods have been suggested to define W t+1 jet production such that large contributions from tt diagrams are removed [51]. We use the diagram removal method as it is the most straightforward to implement in Madgraph5 [52]. Let us denote the amplitude for W t+1 jet by A W t . We can write this as A W t = A 1 + A 2 where A 2 is the amplitude for all diagrams containing tt production. Double counting arises from the appearance of . Diagram removal is implemented by setting A 2 = 0 in our definition of the W t+1 jet process, which removes the double counting. Although this method also neglects the interference term 2Re(A † 1 A 2 ), it has been shown that the effect of this on observables is moderate and that this method approximates W t production well. We will validate this choice by comparing our simulations to data in Sect. 2.2. The violation of gauge invariance in the diagram removal method is found to have no observable effect [51].
The production of a W + W − j j final state in the di-boson channel features overcounting as a result of interference with both tt and W t+1 jet production. This happens only when the two jets originate from b quarks. In our simulations we remove this interference by treating the b quarks as massive and neglecting them from the definitions of the proton and jets i.e. by working in a 4-flavour scheme. This is the method used by ATLAS in their simulations at 13 TeV [53]. 4-flavour parton distribution functions are used. This removes all tt and W t+1 jet production from the di-boson channel, but also neglects processes with initial and final state b quarks which do contribute to di-boson production. A study of how well this 4-flavour scheme approximates the full di-boson production cross section at centre of mass energies √ s = 14,100 TeV was undertaken in Ref. [54] by comparing the leading order cross sections of di-boson production in the 4 and 5-flavour schemes, where it was found that the difference is negligible at 14 TeV and ∼ 5% at 100 TeV.
To produce an accurate simulation of the SM background at each future collider, we use Delphes3 to simulate the response of the detectors and the decay of short-lived particles. Jets are clustered using the anti-k T clustering algorithm [55] with jet radius R = 0.4. This value is chosen from the ATLAS analysis at 13 TeV [42] to reproduce their analysis as closely as possible. To mimic the response of different detectors at each future collider we specify detector configurations as follows. The ATLAS configuration is used in all simulations at 13 TeV. At 14 TeV and 27 TeV we use the Delphes3 HL-LHC configuration designed to reproduce the average response of the ATLAS and CMS detectors at high energies and luminosities. Similarly in our simulations of the 100 TeV FCC-hh we use the FCC-hh configuration provided by Delphes3. We maintain the default settings in our simulations except in the case of muon isolation. Muon isolation is defined by finding the sum of the transverse momentum p T of all objects within a cone of radius R max around a muon, excluding the p T of the muon itself. If the sum satisfies p sum T < p max T for fixed p max T , the muon is considered isolated. At 13 TeV and 14 TeV we select only isolated muons with p max T = 0.2 GeV and R max = 0.2, choosing these parameters to reproduce the 13 TeV ATLAS analysis. At 27 TeV and 100 TeV we make no selection on the muon isolation criteria, following the same reasoning as in [41]. This choice is made because the overall normalisation of the SM background is found to be very dependent on the muon isolation criteria and the specific selection made will likely vary in different future analyses. Relative to our simulations, any selection on muon isolation at future experiments will only reduce the SM background producing a better sensitivity to the LQ signal.
We are interested in the search for TeV-scale LQs, which are expected to manifest as a resonance at high m min (μ, j). Producing a large number of events in the tail of the m min (μ, j) distribution is therefore necessary to achieve good statistics in this region. We find that binning the generation of events in m min (μ, j) at parton-level or in parameters such as the dimuon invariant mass M μμ and H T = p is inefficient for producing a sufficient number of tail events. Instead we reweight the generation of each event x by applying a bias b(x) ∝ s(x) 5 . For each SM background process s(x) is defined at parton-level as the invariant mass of the final state muons and jets, where we only include the minimum number of jets in the multi-jet process definitions of Eq. 7, accounting for jets originating from top quarks. For example, for Drell-Yan we define s(x) as the invariant mass of the di-muon final state. All physical observables and distributions shown in this paper have been obtained by unweighting the events after parton showering and detector simulation, in order to remove the effect of this bias.

Validation
We first validate our methods by simulating the SM background at √ s = 13 TeV for an integrated luminosity L = 3.2 fb −1 and comparing with the ATLAS search for second generation LQs at the same centre of mass energy and integrated luminosity [42]. We compare our simulations to the ATLAS data in two regions of phase space: the preselection region and the signal region. Both are defined by cuts on p T , |η| and R designed to increase the significance of a LQ signal above the SM background and are summarised in Table 2. All jet cuts are placed on the two hardest jets in the event denoted by j 1 , j 2 . The signal region is subject to further cuts on S T = p T and M μμ . These significantly reduce features of the SM background due to soft jets and W and Z boson resonances. In both the signal and preselection regions we also reject muons falling in the range 1.01 < |η| < 1.1 as specified by the ATLAS analysis to avoid potential p T mismeasurement in this range. A pre- liminary selection on muon isolation is made at the level of detector simulation as outlined in Sect. 2.1.
To efficiently simulate events in these regions of phase space, we generate events subject to a subset of the phase space cuts. These are applied at parton-level in the Madgraph5 run card and summarised in Table 3. We will refer to this subset as generator cuts. The jet cut off xqcut required for MLM matching is found for each process in the presence of the generator cuts, as outlined in Sect. 2.1. Note however that we set the parameter auto_ptj_mjj = True for DY and di-boson production, allowing the jet matching procedure to automatically set the cuts on p j T and M j 1 j 2 equal to the chosen value of xqcut. We use xqcut = 30, 60, 60, 30 GeV for DY, top pair, single top and di-boson production, respectively. Figure 8 shows the distribution of preselection events in the parameter M μμ . The Monte Carlo error on each bin is shown in grey and is computed from the Monte Carlo event weights w i by Err i = w 2 i . Systematic uncertainties are not included. Our simulations are not in perfect agreement with the ATLAS data (simulations) shown in black (green). This is expected because we generate all events at leading order, and the dominant process in this region is tt production which has large NLO corrections. However, our simulations provide a good estimation of the order of magnitude of the SM background in each bin.
Our methods are further validated in Fig. 9 which shows the distribution of signal region events in the parameter m min (μ, j). As in the preselection region, we underestimate the SM background slightly by working only at leading order. However, compared to the preselection region, this provides a less fair comparison as normalisation factors have been applied to rescale to the ATLAS simulations in the signal region. Overall we take this comparison as a validation of our methods for simulating the SM background.

Future collider backgrounds
We generate the SM background at the LHC with the full Run II integrated luminosity (13 TeV, 140 fb −1 ) and at the three future colliders in Table 1. We define the signal region at each collider by a set of phase space cuts based on the 13 TeV ATLAS analysis as follows. All angular separations R are kept unchanged. For the LHC and HL-LHC we keep the same cuts on |η|, while at the HE-LHC and FCC-hh these are increased to |η| < 4 to allow for a difference in detec-  Table 6 Values of the jet cut-off parameter xqcut in GeV used in MLM matching of multi-jet events at the 13 TeV LHC and future colliders. All jet matching parameters are found in the presence of the generator cuts summarised in Table 5 Collider  Table 4. As in our 13 TeV simulations, we generate events subject to generator cuts applied at parton-level in the Madgraph5 run card, summarised in Table 5. Table 6 specifies the values of xqcut used in MLM matching for each component of the SM background. Figures 10 and 11 show the resulting distributions of the SM background at the LHC, HL-LHC, HE-LHC and FCChh respectively. As before the Monte Carlo error is shown in grey and systematic uncertainties are not included.

Signal simulations
To find the significance of a LQ at mass m LQ relative to these backgrounds, we simulate the distribution of a LQ resonance in m min (μ, j). We simulate LQ pair production and decay into a μμj j final state at leading order 4 in Madgraph5. We work in a 5-flavour scheme using the 5-flavour NNPDF2.3LO [48] parton distribution function, and correct for parton showering and detector effects using Pythia8 and Delphes3 as before. We simulate LQs from the S 3 LQ model provided by 5 [61]. We specify LQ couplings as follows. For each m LQ the product of couplings |(Y de ) 32 (Y * de ) 22 | is fixed by fits to the NCBAs as given in Eq. 6. We choose ( 32 and set all other (Y de ) i j to zero. This couples the LQ to b L μ L and s L μ L pairs as required by the NCBAs. Following the conventions of Ref. [31], all CKM and PMNS mixing occurs within the up and neutrino sector respectively i.e. we set 4 The leptoquark pair production signal at hadron colliders was calculated some time ago at leading order [56,57]. Next-to-leading order effects [58,59] and parton shower effects [60] on the signal have also been studied. 5 A link to the UFO model files can be found within this reference. V d L = V e L = I in Eq. 5. Our choice of couplings then corresponds to setting (Y L ) 22 = (Y L ) 32 = 0, inducing further couplings of the LQ to u L μ L , b L ν L , s L ν L and u L ν L pairs where u L and ν L denote the vectors of left-handed up-type quarks and neutrinos respectively. These include CKM suppressed couplings to u L μ L and c L μ L which will contribute to the muon-jet decay channel of the LQ, increasing the number of events in the μμj j signal. Similarly, we account for these additional CKM and PMNS suppressed couplings in calculating the theory predictions for σ × BR, taking the central values of V CKM and U PMNS from [62] assuming normal ordering of neutrino masses.
As outlined in Sects. 2.2 and 2.3 we generate events subject to the generator cuts in Table 5 at parton-level, applying the full set of signal region phase space cuts from Table 4 in the analyses. Examples of the predicted distributions of LQ events at each future collider are included in Figs. 10 and 11.

Statistics
Each experiment consists of measurements of events in N bins of a histogram, denoted by n i where i = 1, . . . , N . We find the expected number of background events b i and signal events s i in bin i from our Monte Carlo simulations, and parametrise the signal present in our data sample by the signal strength μ ∈ [0, 1]. The likelihood is defined by taking the product of Poisson probabilities in each bin where θ denotes all nuisance parameters. Defining the profile likelihood ratio as λ(μ) = L(μ,θ)/L(μ,θ) whereμ andθ are the maximum likelihood estimates of μ and θ , andθ is found by maximising the likelihood with fixed μ. Suppose the measured data n i shows no fluctuations above the SM background b i . To set exclusion limits on σ × BR we test the b + μs hypothesis and find the maximum value of μ compatible with the data. We quantify compatibility by computing the p value from the modified frequentist CL s method [63] and the test statistic q μ defined by The upper limit at 95% CL on μ is then given by the value of μ at which CL s = 0.05. We compute the CL s values using pyhf [64], a Python implementation of HistFactory [65]. By comparison with the theoretical predictions, any (σ × BR) theory > (σ × BR) lim can then be excluded. This determines the mass sensitivity i.e. the maximum m LQ that could be excluded at 95% CL at each future collider. Alternatively, suppose an excess of events is seen in the data n i . To find the significance of such an observation we test the compatibility of the background-only hypothesis μ = 0 with the data. The test statistic q 0 is defined by The discovery reach of each future collider is found by determining, for each m LQ of interest, the integrated luminosity L required for a p-value of CL s = 2.9 × 10 −7 or equivalently a statistical significance of 5σ . In determining the discovery and exclusion sensitivities, we work in the large sample approximation and use the Asimov data set to calculate the median CL s [66].
We validate this method by using the signal region data generated at √ s = 13 TeV, L = 3.2 fb −1 in Fig. 9 to place limits on LQs in the range m LQ ∈ [500, 1400] GeV. The resulting limits on σ × BR as a function of m LQ are shown in Fig. 12, excluding LQ masses up to approximately 1.15 TeV. This limit is compared to the exclusion limits found by ATLAS, shown by the black dashed curve, indicating sensitivity to LQ masses up to 1.05 TeV. Note that for the purposes of this comparison only we generate events and compute the Allanach,Corbett and Madigan 2020 √ s = 13 TeV, L = 3.2 fb −1 limits at 95% CL σNLO × BR ATLAS σNLO × BR Expected limit ATLAS expected limit Expected ±1σ Expected ±2σ Fig. 12 Validation plot comparing our expected limits at 95% CL on σ × BR for LQ pair production and decay into a μμj j final state at √ s = 13 TeV, L = 3.2 fb −1 to the expected limits obtained by ATLAS σ ×BR from a model of second generation LQs decaying into a μ − μ + cc final state with coupling y μc = √ 0.01 × 4πα em from the minimal Buchmüller-Rückl-Wyler model [67], following the ATLAS 13 TeV analysis. All other LQ events and values of σ × BR in this paper are found as outlined in Sects.

and 3.4.
This shows that our methods have slightly overestimated the sensitivity to LQ. This is to be expected from the fact that we have underestimated the SM background and do not include systematic uncertainties in setting limits. However, as an estimate of the sensitivity this is a good approximation, and so we take this comparison as a validation of our methods and proceed by using the same methods for future colliders.

Future colliders
The resulting limits on σ × BR as a function of m LQ are shown in Figs. 13 and 14 for the LHC, HL-LHC, HE-LHC and FCC-hh at design integrated luminosities. We compare our limits with theory predictions for σ × BR, shown by the blue curves, and determine the mass to which each collider is sensitive from the point of intersection. We see that while the sensitivity will be increased up to m LQ = 1.75 TeV and eventually m LQ = 2.5 TeV by the LHC Run II and HL-LHC respectively, the HE-LHC and FCC-hh have the potential to explore a much larger range of LQ parameter space, excluding masses up to m LQ = 4.8 TeV and 13.5 TeV respectively.
To further investigate the potential of future colliders to exclude high-mass LQs, we scan over a range of integrated luminosities up to L = 3, 15 and 20 ab −1 for the HL-LHC, HE-LHC and FCC-hh respectively and determine the mass sensitivity at 95% CL for each. Similarly, we perform a scan over integrated luminosities and determine the discovery reach of each future collider. These results are shown in Fig. 17. In both plots the points correspond to the design integrated luminosities of L = 3, 15 and 20 ab −1 . The highest m LQ that can be observed with a 5σ significance Fig. 13 Expected limits at 95% CL on σ × BR for the pair production of LQs and decay into a μμj j final state at the LHC with full Run II integrated luminosity (a) and the HL-LHC (b). Theory curves σ NLO × BR are calculated for narrow width LQs with couplings chosen to fit the NCBAs  is m LQ = 9.5 TeV: we predict that narrow width scalar LQs could be discovered at this mass assuming the FCC-hh operates at the full L = 20 ab −1 . Similarly, the HE-LHC and HL-LHC have the potential to observe narrow width scalar LQs of masses up to m LQ = 3.6 TeV and 1.9 TeV respectively. Finally, we compute the discovery reach of the LHC Run II with √ s = 13 TeV, L = 140 fb −1 to be m LQ = 1.2 TeV, right on the edge of the 95% exclusion limits already found by the 13 TeV LHC as discussed in Sect. 1. Table 7 summarises the maximum 5σ discovery reach and mass exclusion at 95% CL of each future collider.

Wide resonances
If we depart from the NCBA limit, the LQs may acquire an appreciable and non-negligible width. In this subsection, we wish to estimate how big the effect might be on the resulting sensitivity. Our SM background simulations and statistical methods can be applied to determine the approximate change in sensitivity to wider LQ resonances, where some |y lq | may be large, 6 with one important caveat: we do not include interference between signal and background. Interference can only occur between the LQ signal and DY, as all other SM background processes include neutrinos in the final state. Our final estimate of sensitivity will be an overestimate because we may expect the signal to be broadened further by signal-background interference effects. Our purpose however, is just to see the approximate shift in sensitivity rather than to provide a true and accurate calculation of the sensitivity itself. We shall see that the sensitivity is not drastically changed by including large width effects, and we expect that this qualitative conclusion holds once signal-background interference effects have been included. The partial decay width of a LQ into a lepton l and quark q is related to the mass m LQ and coupling y lq by [68] = |y lq | 2 m LQ 16π .
Given the choice of couplings in our signal simulations as outlined in Sect. 3.1, we have so far only considered narrow LQ resonances satisfying /m LQ < 0.01. Any narrow width LQ will still produce a wide resonance in the distribution of m min (μ, j) as shown in Figs. 10 and 11. This is an effect of changes in the kinematics of the final state particles due to parton showering and detector resolution, as well as the ambiguity in defining m min (μ, j), and determines the experimental resolution. By fitting a normal distribution to these resonances and approximating the resolution res by twice the standard deviation, we estimate the resolution to be res /m LQ = 0.1. To investigate the effects of wide resonances we then simulate LQ events with decay width ≥ res . We do this by switching on the same couplings (Y de ) 22 = (Y de ) 32 = 0 as before, determining their values from Eq. 11 for = 0.1, 0.2 and 0.5.  Figure  15b shows the corresponding expected limits on σ × BR at 95% CL for the HE-LHC at 27 TeV, 15 ab −1 . To provide a sample estimate of sensitivity we compare our limits to the values of σ LO × BR calculated for a wide LQ signal with nonzero couplings (Y de ) 22 and (Y de ) 32 as outlined above. Figure 15 shows that the distribution of signal events spreads out in m min (μ, j) with increasing . We expect that the sensitivity to LQs is decreased as a result of the signal events spreading out in this way rather than being peaked around a few bins. This effect is seen in the increase in the upper limits on σ × BR with increasing in Fig. 15b. We can see from the intersection of the /m LQ = 0.1 theory curve (purple, dashed) with each set of expected limits (solid curves) that if the theory predictions for σ × BR were independent of LQ couplings (Y de ) 22 and (Y de ) 32 , an increase in LQ width from /m LQ = 0.1 to /m LQ = 0.5 would result in a loss of sensitivity from approximately m LQ = 4.8 TeV to m LQ = 4 TeV. However, this effect is mitigated by the fact that at such large couplings, pair production is no longer dominated by gluon-gluon interactions. Instead, pair production via quark-lepton interactions has a significant contribution to the total cross section. As a result, the theory prediction for   32 . This can be seen by the overall increase in the number of signal events with in Fig. 15a, and by the large increase in the values of σ × BR with in Fig. 15b. Overall this leads to an increase in sensitivity to LQs with increasing . A similar effect is seen in our predictions for wide LQs at the HL-LHC and FCC-hh, as shown in Fig. 16. Were we to include signal-background interference effects in the calculation, the sensitivity would be degraded. This leads us to conclude that the overall change in sensitivity from the larger widths is not dramatic and would remain small were we to include signal-background interference effects.

Conclusions
We have estimated the exclusion and discovery sensitivities of future hadron colliders to LQ pair production for the case that each LQ decays to a muon and a jet. Such a decay channel is motivated in part by the LQ solution to the NCBAs. It is also motivated by the fact (regardless of the NCBAs) that muons are empirically robust objects, which are good for tagging and beating down irreducible backgrounds. By concentrating on LQ pair production (rather than single LQ production, for example) we cover a large volume of model parameter space where LQs, being perturbatively coupled, are narrow and the pair production cross-section varies only with the LQ mass m LQ . For such LQs, their production is dominated by production from glue-glue interactions, their interactions with initial state quarks being negligible. This is typically true for LQs that have a coupling-mass relation consistent with the NCBAs, but we emphasise that our sensitivities extend beyond this coupling-mass relation more generally, as discussed below.
The previous estimate of the exclusion sensitivity in Ref. [19] extrapolated LHC search limits using two highly dubious approximations. The first being that experimental efficiency and acceptance would not change with centre of mass energy, and the second that LQs are produced exactly at threshold. With respect to the first point, at large m LQ and at high energies (particularly at FCC-hh), the decay products from LQs will be highly boosted resulting in muons collinear to the jets resulting in more muons failing isolation criteria. The muon momentum resolution is also likely to be very poor at higher energies, since such hard muons will only be bent to a limited extent by the magnets. This also affects signal efficiency due to peak broadening. Secondly, the assumption LQs are produced exactly at threshold is likely to introduce large uncertainties. 7 We rectify these two bad approximations in our present paper by performing a fast simulation of the signal and detector response. We summarise our expected discovery and exclusion sensitivities in Table 7. Reference [19] estimated that the HL-LHC could exclude 2.2 TeV at 95% CL, to be compared with 1.8 TeV. The HE-LHC was estimated to cover up to 4.1 TeV, but this was for a higher centre of mass energy (33 TeV) and a different luminosity (15 ab −1 ), precluding a direct comparison. The FCC-hh exclusion sensitivity was calculated at an integrated luminosity of 10 ab −1 to be 12.0 TeV, to be compared with 12.5 TeV from our estimate (see Fig. 17). It is somewhat surprising that the comparable estimates are so similar, since as we have argued, the old ones were based on self-admitted bad approximations. The results in Table 7 are on a much firmer footing. This is the first time that the 5σ S 3 discovery sensitivities for future colliders have appeared. It is also the first time that S 3 sensitivity estimates for varying luminosities have been calculated as in Fig. 17.
The sensitivities phrased in terms of LQ mass have a dependence on the LQ-fermion couplings assumed in the model, since these may affect the BR of the muon-jet decay rate. However, all limits in the narrow LQ limit on σ × B R also apply to models with different (but still small) LQ couplings to fermions. Only when one or more of the LQ couplings approaches the non-perturbative régime does the LQ width become comparable to the experimental resolution, potentially affecting sensitivity. To cover this case, we con-sidered a wider LQ: see Sect. 3.4. Of and by itself, the width does not change the sensitivity much. Increasing the width divided by mass of the LQ from 0.1 to 0.5 but keeping the cross-section times branching ratio constant only results in a 10% degradation or so in FCC-hh mass reach, as the righthand panel of Fig. 15 shows.
We hope that the results of our study will be useful for the current European Strategy in Particle Physics [69] and provide a part of the physics case for future hadron colliders [37][38][39][40]. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .