keywords

1 Introduction

The nature of dark matter (DM) remains one of the outstanding mysteries in the particle physics today. Among the many possible particle candidates known [1], weakly interacting massive particles (WIMP) are arguably the most theoretically motivated and well-studied scenarios; for a recent review, see [2]. WIMPs with masses in the range of 10 GeV–20 TeV [3] can be stable on the age of the universe, and once they are in thermal equilibrium in the early universe remain so, even after decoupling occurs. Moreover, such stable particles are naturally present in many beyond the Standard Model (BSM) frameworks. For example: In weak-scale supersymmetric (SUSY) theories [4, 5], if one assumes R-parity conservation, then lightest SUSY particle is stable, and hence is a potential component of DM. In the universal extra dimension (UED) models [6], the lightest Kaluza-Klein excitations of neutral electroweak bosons can be viable DM candidates.

Notably, a number of these new physics scenarios also involve additional heavy colored particles, \(\mathcal {Q}\), that couple to DM candidate(s). Hence, if \(\mathcal {Q}\) and DM are kinematically accessible at the Large Hadron Collider (LHC), or its potential successors, such as the High Energy (HE)-LHC, then it may be possible to study DM in a laboratory setting. In particular, once produced and if allowed, \(\mathcal {Q}\) can decay into DM and SM particles leading to a plethora of interesting signatures at the LHC.

At hadron colliders, search strategies for these hypothetical colored particles usually involve investigating jets and leptons produced in association with final-state DM candidates manifesting as large missing transverse energy (\(E_{\mathrm {T}}^{\mathrm {miss}}\)). In the context of simplified SUSY models, such signatures are now strongly constrained by LHC data if \(\mathcal {Q}\) has a mass around 1 TeV and below [7,8,9,10]. Such constraints on \(\mathcal {Q}\), however, can be circumvented. One of the most celebrated examples of this is the compressed spectrum scenario [11,12,13]. In this situation, the DM and new colored particles have a small mass splitting. Consequently, the visible decay products in the \(\mathcal {Q}\rightarrow \)DM+SM process do not have sufficient momenta to be readily distinguished from SM backgrounds. In other words, the compression of the mass spectrum constrains the visible decay products of \(\mathcal {Q}\) to possess such low momenta that they fail experimental selection criteria. This leads to significantly smaller selection and acceptance efficiencies, and hence significantly weaker bounds on heavy particle masses.

Even though a compressed mass spectrum represents a special corner of a typical BSM parameter space, the most attractive feature of this situation is that it allows relatively light, colored DM partners in light of present-day LHC data. This is particularly true for stops (\({\tilde{t}}\)) and gluinos (\({\tilde{g}}\)) in SUSY [14,15,16,17]. If a compressed scenario is realized in nature, then one can experimentally resolve the soft, i.e., low \(p_{T}\), visible decays of \(\mathcal {Q}\) by recoiling against a relatively hard, i.e., high \(p_{T}\), electroweak or QCD radiation that, in its own right, is sufficiently energetic to satisfy trigger criteria. One such process, shown diagrammatically in Fig. 1 and the focus of this study, is the inclusive monojet plus \(E_{\mathrm {T}}^{\mathrm {miss}}\) collider signature.Footnote 1

Fig. 1
figure 1

Diagrammatic depiction of \(\mathcal {Q}\overline{\mathcal {Q}}\) pair production with an extra hard QCD radiation in pp collisions

Were evidence for the new particle \(\mathcal {Q}\) established at the LHC, or a successor experiment, it would be crucial to determine \(\mathcal {Q}\)’s properties, especially its mass, spin, and color representation. Generically, such a program would involve investigating various collider observables that can discriminate against possible candidates for \(\mathcal {Q}\). For example: cross sections are highly sensitive to the aforementioned spacetime and internal quantum numbers. Consider the cases of a scalar top \({\tilde{t}}\) (color triplet), a spin-1/2 top partner \(T_p\) (color triplet), and spin-1/2 gluino \({\tilde{g}}\) (color octet). For a fixed mass, i.e., \(m_{{\tilde{t}}} =m_{T_p} =m_{{\tilde{g}}}\), the pair production cross sections for these particles exhibit the hierarchy

$$\begin{aligned} \sigma ({\tilde{t}}{\tilde{t}}^*) \ll \sigma (T_p\bar{T}_p) \ll \sigma ({\tilde{g}}{\tilde{g}}). \end{aligned}$$
(1)

Conversely, for a fixed cross section, i.e., \(\sigma ({\tilde{t}}{\tilde{t}}^*)=\sigma (T_p\bar{T}_p)=\sigma ({\tilde{g}}{\tilde{g}})\), one finds that

$$\begin{aligned} m_{{\tilde{g}}}> m_{T_p} > m_{{\tilde{t}}}. \end{aligned}$$
(2)

This implies, however, that were a monojet cross section measured, the result could be replicated by different \(\mathcal {Q}\) scenarios by a simple tuning of mass \(m_{\mathcal {Q}}\). In another way, one cannot constrain the mass of \(\mathcal {Q}\) from cross section measurements alone without first asserting its color representation and spin. Quantitatively, this is more nuanced due to fact that leading order (LO) calculations are poor approximations for QCD processes, even when using sensible scale choices. Using Ref. [19], one can easily verify that, like the top quark [20, 21], QCD corrections at next-to-leading order (NLO) increase the production cross section \(\sigma (T_p\bar{T}_p)\) by \(\mathcal {O}(50)\%\) for TeV-scale \(T_p\). This is the case at both \(\sqrt{s} = 14\) and 27 TeV, and despite scale uncertainties at LO and NLO spanning \(\mathcal {O}\)(20–30)% and \(\mathcal {O}(10)\%\), respectively, for \(\sigma (T_p\bar{T}_p)\). Moreover, it is well-known that next-to-next-to-leading order (NNLO) corrections are non-negligible for SM top production [22]. It is also known that such large theoretical uncertainties can greatly limit the interpretation of the experimental results, particularly in searches for so-called top-philic dark matter [23].

The situation, however, is more hopeful following the advent of general-purpose, precision Monte Carlo (MC) event generators. With software suites such as Herwig [24], MadGraph5_aMC@NLO+Pythia8 (MG5_aMC@NLO +PY8) [19, 25], and Sherpa [26], automated event generation at NLO in QCD with parton shower (PS) matching is now possible for both SM and BSM [27] processes. Not only can one now readily include potentially important \(\mathcal {O}(\alpha _s)\) corrections to cross sections normalizations, but parton showers augment fixed order predictions with resummed corrections to at least the leading logarithmic (LL) level. As a consequence, such observables like the associated jet multiplicity in the monojet process, an exclusive observable that is critical to search strategies, is automatically modeled at LO+LL accuracy. This is the lowest order at which the quantity is qualitatively correct. In light of the availability of such sophisticated technology, one is now in position to systematically investigate the impact of QCD corrections on the inclusive monojet process.

In this report, we perform such a dedicated precision study on the inclusive monojet signature in the context of a compressed mass spectrum. As mentioned, observables associated with this process are highly sensitive to the mass, spin, and color representation of the mediating states. Hence, we consider benchmark models with representative mass, spin and color configurations for \(\mathcal {Q}\), with \(\mathcal {Q}\in \{{\tilde{t}},T_p,{\tilde{g}}\}\). Our study is aimed at the HL-LHC, assuming \(\mathcal {L} = 3\mathrm{~ab^{-1}}\) at \(\sqrt{s}=14\) TeV, and the HE-LHC, assuming \(\mathcal {L} = 15\mathrm{~ab^{-1}}\) at \(\sqrt{s}=27\) TeV. We show how the monojet search strategy can be used to identify the nature of \(\mathcal {Q}\) from the \(p_{T}\) dependence of the leading jet in \(\mathcal {Q}\overline{\mathcal {Q}}\) pair production. We also estimate the precision required to distinguish these new physics scenarios. We quantitatively discuss various sources of theoretical uncertainty, including event generator dependence. Although we do not discuss the second jet distribution intensively in this paper, azimuthal angle correlation of the first and second jet contains information of the spin of \(\mathcal {Q}\), as discussed in [28, 29].

The remainder of this study continues in the following manner: In Sect. 2 we provide in-dept detail of our computational setup. In Sect. 3, we discuss observed and expected sensitivity of monojet searches at present and hypothetical future facilities, and address various theoretical uncertainties in Sect. 4. A brief outlook on the impact of this work is discussed in Sect. 5, and we conclude in Sect. 6.

2 Computational and theoretical setup

Systematic studies of QCD radiation in the production of hypothetical, TeV-scale colored particles are now possible due to the availability of precision, general-purpose MC event generators. In practice, this nontrivial task is handled by using several individually published formalisms and software packages that have largely been integrated into a single framework or well-specified tool chain sequence. In this section, we describe our computational and theoretical setups for modeling \(\mathcal {Q}\overline{\mathcal {Q}}\) production with various associated jet multiplicities at LO+PS and NLO+PS in MG5_aMC@NLO +PY8 and Sherpa. For numerical results, readers can go directly to Sect. 3.

The section continues as follows: In Sect. 2.1, we briefly summarize the \(\mathcal {Q}\) we consider from representative BSM scenarios. In Sect. 2.2, we enumerate the several methods for incorporating additional QCD radiation into \(\mathcal {Q}\overline{\mathcal {Q}}\) production that we employ, and briefly note their main features and formal accuracies. We describe our setup for MG5_aMC@NLO +PY8 and Sherpa, respectively, in Sects. 2.3 and 2.4, and our detector simulation in Sect. 2.5. In Sect. 2.6, we summarize the SM inputs.

2.1 Framework for new heavy colored particles

In this analysis, we consider three benchmark BSM candidates for \(\mathcal {Q}\): a stop squark \({\tilde{t}}\), a gluino \({\tilde{g}}\), and a fermionic top partner \(T_p\). To model these states in pp collisions at our desired accuracy, we use the NLO in QCD-accurate Universal FeynRules Object (UFO) [30, 31] model libraries available from the FeynRules model database [32]. The \(\mathcal {O}(\alpha _s)\) counterterms required for NLO computations and contained in these libraries are generated with FeynRules [30], using NLOCT [27] and FeynArts [33]. For illustrative purpose, we choose three mass values for \(\mathcal {Q}\), namely \(m_{\mathcal {Q}}\) = 400, 600, and 800 GeV. We note that as the spacetime and \(\mathrm {SU}(3)_c\) quantum numbers for \(T_p\) are identical to those of the SM top quark, several publicly available calculations can be adapted in straightforward ways for \(T_p\). This includes total cross section predictions for inclusive \(pp\rightarrow T_p\bar{T}_p\) production at NNLO in QCD, which we obtain using the HATHOR package [34]. In Sect. 3.3, where we study the dependence of the selected observables on the color and spin representations of \(\mathcal {Q}\), we also briefly consider the well-studied [35,36,37] case of a scalar gluon \((\sigma )\) at \(m_\sigma = 600\mathrm{~GeV}\).

Table 1 Summary of signal particles, their SU\((3)_c\) and Lorentz representations (Rep.), and decay mode to stable DM candidate \((X)\)

Furthermore, in order to focus on the ISR from \(\mathcal {Q}\overline{\mathcal {Q}}\) production and also to use traditional analysis techniques, we assume that the decay of \(\mathcal {Q}\) is prompt, with its characteristic lifetime (\(\tau _{\mathcal {Q}}\)) satisfying \(d_0 = \beta c\tau _{\mathcal {Q}} \ll 100~\mu m\), where \(\beta \) denotes the relative velocity of \(\mathcal {Q}\). As a result, \(\mathcal {Q}\)’s total width must respect the boundary \( \varGamma _{\mathcal {Q}} = \hbar / \tau _{\mathcal {Q}} \gg 2~\mathrm{meV}\). In realizations of the \(\mathcal {Q}\) we consider, e.g., Refs. [38,39,40], this stipulation on \(\varGamma _{\mathcal {Q}}\) is readily satisfied by large regions of the models’ phenomenologically relevant parameter spaces. We also assume that the \(\mathcal {Q}\rightarrow X+\mathrm{SM}\) decay is \(\mathcal {Q}\)’s dominant decay mode. By virtue of DM being weakly coupled, this implies that \({\mathcal {Q}}\)’s total width scales as \(\varGamma _{\mathcal {Q}}\sim g_{\mathcal {Q}}^2 m_{\mathcal {Q}}\), with an effective \(\mathcal {Q}\)-DM coupling \(g_{\mathcal {Q}}\ll 1\). This ensure that the width-to-mass adheres to the inequality

$$\begin{aligned} 1 \gg \frac{\varGamma _{\mathcal {Q}}}{m_{\mathcal {Q}}} \gg \left( \frac{2~\mathrm{meV}}{m_{\mathcal {Q}}}\right) \approx 3\times 10^{-15} \times \left( \frac{600\mathrm{~GeV}}{m_{\mathcal {Q}}}\right) . \end{aligned}$$
(3)

Subsequently, we decouple the production of \(\mathcal {Q}\) from its decay into an invisible particle \(X\) using the narrow width approximation (NWA), as implemented in MadSpin [41]. For technical simplicity, we neglect spin correlation between the production and decay of \(\mathcal {Q}\). For \({\tilde{t}}\) and \(\sigma \) decays, this has no impact due to their scalar nature. For \(T_p\) and \({\tilde{g}}\) decays, this can impact the \(p_T\) and \(\eta \) distributions of X and subleading jets associated with \(\mathcal {Q}\)’s decay, and hence the true \(E_{\mathrm {T}}^{\mathrm {miss}}\) distribution. However, we have checked that this approximation has a relatively small impact (up to 5%) on our specific study since the \(E_{\mathrm {T}}^{\mathrm {miss}}\) observable we use is built from visible objects, and is dominated by contributions from the hard radiation.

More specifically, for each \(\mathcal {Q}\), the particle nature of the DM particle \(X\) is chosen in accordance with its underlying theory: For \(T_p\), \(X\) is a neutral scalar and \(T_p\) are decayed to the two-body final state \(X+q\), where \(q\) is a light QCD quark. For \({\tilde{t}}\), \(X\) is a neutral fermion and \({\tilde{t}}\) are decayed via an off-shell top quark to the four-body final state, \({\tilde{t}}\rightarrow t^* X\rightarrow X+bq\bar{q}'\). For \({\tilde{g}}\), \(X\) is a neutral fermion and \({\tilde{g}}\) are decayed to the three-body final state \(X+q\bar{q}'\). We enforce a compressed mass spectrum by fixing a small mass gap between \(\mathcal {Q}\) and \(X\) to be \(\varDelta m = m_{\mathcal {Q}} - m_{X} = 20\; \mathrm {GeV}\). Since the mass gap is (relatively) small, the SM decay products of \(\mathcal {Q}\) are forced to be (relatively) soft. Hence, the SM decay products of \(\mathcal {Q}\) fail the criteria needed to be identified as the leading jet, further justifying our neglect of spin-correlation in the decays of \(T_p\) and \({\tilde{g}}\). Qualitatively and quantitatively, our results are expected to be robust against varying \(\varDelta m\) factors of 2 as this does not change substantially the likelihood of \(\mathcal {Q}\) decay promptly, nor the likelihood of its decay products being tagged as subleading jets in the event. For pathologically smaller \(\varDelta m_{\mathcal {Q}}\), the lifetime of \(\mathcal {Q}\) can be extended beyond \( 100~\mu m\), but does not change the prompt monojet collider signature, nor our conclusions. For further discussion on displaced (monojet) signatures, see [42, 43]. For hierarchically larger \(\varDelta m_{\mathcal {Q}}\), the SM decay products of \(\mathcal {Q}\) become more energetic and the process transitions to a multi-jet+\(E_{\mathrm {T}}^{\mathrm {miss}}\) collider signature. The signature associated with the latter scenario is outside the scope of our study and hence will not be discussed further.

The choices of \(\mathcal {Q}\), their relevant quantum numbers and decay path, and the corresponding UFO library references are summarized in Table 1.

2.2 Multi-leg matching and merging prescriptions

The collider signature considered in this work is characterized by the presence of a high-\(p_{T}\) jet recoiling against the \(\mathcal {Q}\overline{\mathcal {Q}}\) system. It is thus necessary to include at least one QCD radiation at the matrix element level beyond what is modeled in Born-level \(\mathcal {Q}\overline{\mathcal {Q}}\) production. With presently available MC technology, this can be achieved in different ways and at various formal accuracies. We now briefly describe the several prescription used in this study.

  • Leading Order Multijet Merging: The LO multijet/multileg merging techniques [44,45,46] outline how parton shower emissions can be augmented with full matrix elements. The emissions are classified according to their hardness, i.e., \(p_{T}\), and in terms of a dimensionful variable \(Q_{\mathrm {cut}}\). Emissions above a hardness threshold \(Q_{\mathrm {cut}}\) are described at LO accuracy using the appropriate matrix elements while preserving the all-orders resummation accuracy of the parton shower below \(Q_{\mathrm {cut}}\). In this work we use the MLM scheme [45] as implemented in MadGraph5_aMC@NLO 2.6.0 (MG5_aMC@NLO) [19]. We take into account matrix element corrections to \(\mathcal {Q}\overline{\mathcal {Q}}\) pair production in association with up to two QCD partons. While genuine \(\mathcal {O}(\alpha _s)\) (and higher) corrections are included via this procedure, the calculation remains formally LO accurate (LO+LL after parton showering) due to missing virtual corrections. Some of the earlier studies on monojet spectra with LO multijet merging technique using then available Madgraph/MadEvent v4 [47], we refer [48,49,50,51,52].

  • \(\mathcal {Q}\overline{\mathcal {Q}}\) Production at Next-to-Leading Order with Parton Shower Matching: The accuracy of fixed order (FO) matrix element calculations at NLO can be combined with resummed parton showers at LL by means of NLO+PS matching techniques [53,54,55], such as the MC@NLO formalism [53]. In this approach, the LO matrix element for an extra hard, wide-angle QCD emission in the final state is naturally included as part of the \(\mathcal {O}(\alpha _s)\) FO correction. Extra soft and/or collinear emissions enter through the parton shower. Potential double counting of \(\mathcal {O}(\alpha _s)\) soft/collinear contributions is avoided by the use of additional counter terms. It is worth noting that leading-jet observables in this calculation are at most LO+LL accurate, but are nonetheless well-defined at all \(p_{T}\).

  • \(\mathcal {Q}\overline{\mathcal {Q}}+j\) Production at Next-to-Leading Order with Parton Shower Matching: In order to achieve NLO+LL accuracy for leading-jet observables, the above NLO+PS matching technique must be applied to the \(\mathcal {Q}\overline{\mathcal {Q}}+j\) process. This requires explicitly introducing a regularizing \(p_{T}\) selection on the leading jet at the matrix element level. Since no Sudakov form factor is present for this jet, the \(p_{T}\) selection must also be well above the Sudakov shoulder of the inclusive \(\mathcal {Q}\overline{\mathcal {Q}}\)-system to ensure that the FO is perturbatively valid.

  • Next-to-Leading Order Multijet Merging: The LO multijet merging technique described above can be extended to describe jet observables at NLO precision for jets above a hardness scale \(Q_\mathrm {cut}\). Analogous to the LO case, matrix element merging with parton shower matching at next-to-leading order, known colloquially as MEPS@NLO, is achievable by introducing additional all-orders/resummed Sudakov form factors for each NLO-accurate matrix element in consideration. We use an extension of the Catani-Krauss-Kuhn-Webber (CKKW) [56,57,58] merging formalism as implemented in Sherpa by Refs. [59, 60]. In the following, we employ MEPS@NLO multijet merging with up to one or two jets, meaning that samples will contain up to two or three real radiations, respectively, beyond the lowest order process before parton showering. As in the LO(+LL) case, while \(\mathcal {O}(\alpha _s^2)\) corrections are present in this calculation, the final result remains formally NLO(+LL) accurate.

2.3 Event generation in MadGraph5_aMC@NLO + Pythia8

Cross section calculations and event generation with accuracy up to NLO in QCD are handled using MG5_aMC@NLO. For signal processes, we use  the NLO-accurate UFO libraries described in Sect. 2.1 and listed in Table. 1 as inputs to MG5_aMC@NLO. Within the MG5_aMC@NLO framework, one-loop virtual corrections are evaluated numerically via MadLoop [61] and matched with real emissions using the Frixione-Kunszt-Signer (FKS) subtraction formalism [62], as implement by Ref. [63]. Decays of \(\mathcal {Q}\) are then handled at LO accuracy with MadSpin [41]. The central value \(Q_0\) for the renormalization scale \(Q_R\) and factorization scale \(Q_F\) is set to

$$\begin{aligned} Q_F,~Q_R = Q_0 \equiv \frac{H_T}{2}, \;\text {with}\; H_T \equiv \sum _{k\in \{g,\overset{(-)}{q},\overset{(-)}{\mathcal {Q}}\}} \sqrt{m_k^2 + \vert \mathbf {p}_{T}^k \vert ^2}, \end{aligned}$$
(4)

where \(H_T\) is the scalar sum of the transverse energy of the final state partons and \(\mathcal {Q}\).

In Sect. 3, where we discuss experimental searches and sensitivity to \(\mathcal {Q}\) in monojet searches, we impose an analysis-level selection on the leading jet \((j_1)\) in an event. More specifically, we set the ptj variable in MG5_aMC@NLO’s run_card.dat file to ptj\(=p_{T,\mathrm {cut}}^{j_1}-100\;\mathrm {GeV}\) (and jet radius \(R=0.6\)). Then, to enhance yields at relatively high \(p_{T}\), we generate events by binning the phase space in \(p_{T}\). After preparing event samples for a particular \(p_{T,\mathrm {cut}}^{j_1}\), we apply higher \(p_{T,\mathrm {cut}}^{j_1}\) until the statistical uncertainty of the MC samples is no longer negligible. At this point, we prepare another event sample based on \(p_{T,\mathrm {cut}}^{j_1}\) and apply the procedure iteratively. For the samples with the highest \(p_{T,\mathrm {cut}}^{j_1}\), we apply exponential biasing on \(p_{T}^{j_1}\) at the event-generation level to enhance the tail of jet \(p_{T}\) distributions.

After the \(\mathcal {Q}\overline{\mathcal {Q}}\) pair have been decayed, events are passed to Pythia 8.2.26 [25] for parton showering and hadronization. We choose a shower starting scale \(Q_S\) small enough such that light, colored final-state partons in the matrix element remain the hardest emissions in the full process, if the parton exists. Namely, the parton shower is restricted to operate below the scale [19]

$$\begin{aligned} Q_S = Q_{S*} \approx \min \left[ \frac{H_T}{2}, \sqrt{d_*} \right] . \end{aligned}$$
(5)

Here, \(d_*\) is the minimum \(d_{i,j}\) distance measure as calculated with the \(k_T\) algorithm [64, 65] (for \(R=1\)) over all momentum recombinations of light, colored partons during the clustering phase in fixed order event generation. The events are essentially categorized by whether or not a hard \(\mathcal {O}(\alpha _s)\) emission is present.

To quantify and estimate the size of missing, higher order QCD corrections, we compute the three-point scale-variation envelope. This is obtained in the usual fashion, i.e., by varying discretely and jointly the factorization and renormalization scales \(Q_R\) and \(Q_F\) over the range,

$$\begin{aligned} 0.5\times Q_0 \le Q_F,~Q_R \le 2.0\times Q_0. \end{aligned}$$
(6)

Where necessary, we also consider the uncertainty associated with the parton shower starting scale. We quantify this by discretely and independently computing the scale variation over the range,

$$\begin{aligned} \min \left[ 0.5 \times \frac{H_T}{2}, \sqrt{d_*} \right] \le Q_S \le \min \left[ 2.0 \times \frac{H_T}{2}, \sqrt{d_*} \right] . \end{aligned}$$
(7)

2.4 Event generation in Sherpa

For cross-validation of the fermionic top quark partner \(T_p\), we employ Sherpa 2.2.4 [26]. At LO in QCD, arbitrary BSM models can be simulated through Sherpa’s generic UFO model [31] interface [66]. At NLO, processes involving \(T_p\) can be simulated by slightly modifying the default SM model file, and setting the top quark mass to the mass of \(T_p\). For the decay into a scalar dark matter particle, an additional decay vertex is added using the methodology of Ref. [66]. Tree-level matrix elements of the calculation are provided by Sherpa’s in-house matrix element generators AMEGIC [67] and COMIX [68]. One-loop amplitudes are treated by interfacing with OpenLoops [69]. Parton shower-matching is performed according to the MC@NLO formalism [53, 70, 71], using Sherpa’s Catani-Seymour subtraction-based shower procedure [72, 73].

To study potential improvements to modeling fermionic top quark partners, we also employ multijet merging at NLO in QCD with Sherpa. To account for additional high-\(p_{T}\) QCD emissions at the matrix-element level, beyond what is already present in inclusive \(T_p\bar{T}_p\) production at NLO, we include the NLO-accurate matrix elements for the processes

$$\begin{aligned} pp\rightarrow T_p\bar{T}_p+j\quad \text {and}\quad pp\rightarrow T_p\bar{T}_p+2j. \end{aligned}$$
(8)

We merge these samples with the fully inclusive \(pp\rightarrow T_p\bar{T}_p\) sample following the MEPS@NLO prescription [59, 60]. In this scheme, the nominal values for the factorization, renormalization, and parton shower scale are determined through a backward clustering procedure that maps higher multiplicity configurations to a \(2\rightarrow 2\) configuration. We set all scales to the invariant mass of the top partner pair. As a nominal value for the merging scale \(Q_{\mathrm {cut}}\), we set \(Q_{\mathrm {cut}}= 120\;\mathrm {GeV}\). For parton showering we employ one of Sherpa’s dipole showers, which is published in [73].

In addition to MEPS@NLO merged samples we also simulate \(pp\rightarrow T_p\bar{T}_p+j\) at NLO+PS with Sherpa. We do not add matrix element corrections to these samples beyond what is already present at NLO+PS. For the generation of the \(pp\rightarrow T_p\bar{T}_p+j\) process, we use the scale schemes,

$$\begin{aligned} Q_F = Q_R = m(T_p\bar{T}_p)\, \quad \text {and}\quad Q_S = p_{T}(T_p\bar{T}_p), \end{aligned}$$
(9)

where \(m(T_p\bar{T}_p)\) and \(p_T(T_p\bar{T}_p)\) denote the \((T_p\bar{T}_p)\)-system’s invariant mass and transverse momentum respectively.

2.5 Detector simulation and object reconstruction

For fast detector simulation, we use Delphes 3.3.3 [74] with the default ATLAS card. Jets are constructed from calorimeter tower elements using Fastjet 3.2.1 [75], according to the anti-\(k_T\) jet clustering algorithm [76] with jet radius \(R=0.4\). Analysis-quality jets are required to satisfy the fiducial and kinematic criteria,

$$\begin{aligned} p_{T}^j > 20\mathrm{~GeV}\quad \text {and}\quad \vert \eta ^j\vert < 4.5. \end{aligned}$$
(10)

Events are then accepted or rejected based the monojet selection criteria discussed in Sect. 3.

2.6 Standard model inputs

We assume \(n_f=5\) active/massless quark flavors and a diagonal Cabbibo-Kobayashi-Maskawa (CKM) quark mixing matrix with unit entries. Relevant SM inputs used in our study include,

$$\begin{aligned} m_t(\mathrm pole) = 173.3\mathrm{~GeV}\quad \text {and}\quad \alpha ^\mathrm{\overline{MS}}(M_{Z})= 0.118. \end{aligned}$$
(11)

The evolution of parton distribution functions (PDFs) and the strong coupling constant \(\alpha _s(\mu _R)\) are extracted using the LHAPDF 6.1.6 [77] libraries. As discussed in Sect. 2.2, LO (NLO) multijet merging with two (one) additional partons accounts for new kinematic channels and configurations that first arise at NNLO. (In principle, all one would need to achieve NNLO accuracy are the missing two-loop virtual corrections.) However, such \(\mathcal {O}(\alpha _s^2)\) contributions are already accounted for in the normalizations of NLO PDFs. Hence, to minimize potential double counting of initial-state contributions, we use the NNPDF 3.0 NNLO PDF set (lhaid=261000) [78] for all signal process calculations. For the LO SM \(pp\rightarrow Z + j\) calculation in Sect. 3, we use the NNPDF 3.0 LO PDF.

3 Monojet searches at the HL- and HE-LHC

At hadron colliders, the term “monojets” represents a broad class of sensitive collider signatures and search strategies that assume varying degrees of particle multiplicity and inclusiveness. In this section, we consider specifically the inclusive monojet signature, as implemented by ATLAS during Run II of the LHC’s operations after collecting \(\mathcal {L}=\)36.1 fb\(^{-1}\) of integrated luminosity at \(\sqrt{s}=13\) TeV [18]. After discussing various sources of experimental and theoretical uncertainties, we report the observed and expected sensitivity of the channel at current and proposed pp colliders using the \(\mathrm {CL}_{\mathrm {s}}\) modified frequentist approach [79]. The model-independent limits derived in Sect. 3.1 are then applied in Sect. 3.2 to the heavy colored particles \(\mathcal {Q}\) described in Sect. 2.1. The ability to determine and distinguish principle properties of \(\mathcal {Q}\) is then discussed in Sect. 3.3.

3.1 Monojet searches, LHC data, and model-independent limits

Table 2 The expected number of SM background events and associated errors for the inclusive mode signal regions IM1-IM10 as defined by the ATLAS experiment in Ref. [18]

The ATLAS and CMS collaborations have both reported on their search of early Run II data for anomalous events with significant transverse momentum imbalance and at least one energetic jet [18, 80]. For the case under consideration [18], the ATLAS collaboration has investigated two overlapping signal regions, categorized as exclusive modes (EM) and inclusive modes (IM), based on various \(E_{\mathrm {T}}^{\mathrm {miss}}\) thresholds spanning \(E_{\mathrm {T}}^{\mathrm {miss}}=250\) GeV to 1 TeV. The EM signal regions are defined in terms of \(E_{\mathrm {T}}^{\mathrm {miss}}\) binning. For example: signal region EM1 (EM5) selects for events with \(E_{\mathrm {T}}^{\mathrm {miss}}\) satisfying 250 GeV \(< E_{\mathrm {T}}^{\mathrm {miss}}< \) 300 GeV (500 GeV \(< E_{\mathrm {T}}^{\mathrm {miss}}< \) 600 GeV). The IM signal regions are defined in terms of minimum \(E_{\mathrm {T}}^{\mathrm {miss}}\) selections. For example: signal region IM1 (IM5) selects for events with \(E_{\mathrm {T}}^{\mathrm {miss}}\) satisfying \(E_{\mathrm {T}}^{\mathrm {miss}}>\) 250 (500) GeV. Additionally, events are required to satisfy the following selection criteria:

  • At least one analysis-level jet with \(p_T>250\) GeV and \(\vert \eta \vert <2.4\).

  • A maximum of four analysis-level jets satisfying \(p_T>30\) GeV and \(\vert \eta \vert <2.8\).

  • An azimuthal separation of \(\varDelta \phi (jet, E_{\mathrm {T}}^{\mathrm {miss}}) >0.4\) for each analysis-level jet and the \(E_{\mathrm {T}}^{\mathrm {miss}}\) vector.

For the remainder of this study, we focus on the IM monojet signal regions.

In Table 2, we display the expected number of SM (background) events passing all selection criteria with uncertainties (statistical and total) in each of the inclusive mode signal regions (IM1-IM10), as reported by ATLAS [18]. Non-statistical errors include both experimental and theoretical uncertainties. Sources of systematic uncertainty estimated in Ref. [18] include: dependencies on parton shower and PDF modeling, which span \(\pm \,\) 0.7 to \(\pm \,\) 0.8%; uncertainties in jet energy and \(E_{\mathrm {T}}^{\mathrm {miss}}\) scales, which range \(\pm \,\)0.5% (IM1) to \(\pm \,\)5.3% (IM10); jet quality and pileup descriptions additionally provide uncertainties ranging \(\pm \,\) 0.8 to \(\pm \,\) 1.8%; for more detailed discussions, see Ref. [18]. We note that the total errors for IM1-5 are nearly flat, with 2.2–2.6%, indicating that these signal regions’ uncertainties are systematics dominated. Statistical and systematic uncertainties are much larger in the higher \(E_{\mathrm {T}}^{\mathrm {miss}}\) regions. However, with the HL phase of LHC, one expects to collect significantly more data that will correspondingly reduce statistical errors for the high-\(E_{\mathrm {T}}^{\mathrm {miss}}\) regions. Additionally, the analysis’ control sample will also increase during the HL run, therefore also reducing systematic uncertainties. Thus, one anticipates that total uncertainties will shrink for future inclusive monojet searches at LHC.

Table 3 Model-independent 95% CL upper limit on the visible cross section for each inclusive mode (IM) signal region, after \(\mathcal {L}=\)36.1 fb\(^{-1}\) of data at \(\sqrt{s} = 13\) TeV, as reported by Ref. [18], and the estimated \(\mathrm {CL}_{\mathrm {s}}\) limit assuming \(\mathcal {L}=\)3 ab\(^{-1}\)

The non-observation of data deviating significantly from SM predictions enables one to set model-independent upper limits on the production cross section of new particles. In the second and third columns of Table 3, we tabulate the expected and observed limits at the 95% confidence level (CL) on the visible cross section,Footnote 2 respectively, for each IM signal region, as reported in [18] with \(\mathcal {L}=\)36.1 fb\(^{-1}\) of data. To address the prospect of the HL-LHC with \(\mathcal {L} = 3\) ab\(^{-1}\) at \(\sqrt{s}\) = 13 TeV, we calculate the expected upper limits by scaling the number of events at the 13 TeV run of LHC. We choose two values of the total systematic uncertainty, namely \(\delta _\mathrm{Sys.}=\)2.5% and 1%. The former is a pessimistic assertion that systematic uncertainties will not be reduced beyond present levels (see Table 2), even after 15–20 years of LHC operations. The latter is an optimistic, but benchmark, assumption. The likelihoods of background only and signal-plus-background hypotheses are set as Gaussian, with a standard deviation set to the total uncertainty. We have checked that our likelihoods are in agreement with those reported by Ref. [18] for \(\mathcal {L} = 36.1\) fb\(^{-1}\). We report the scaled limits in the last two columns of Table 3. We now discuss the impact of these limits on the production of \(\mathcal {Q}\overline{\mathcal {Q}}\) pairs in pp collisions.

3.2 HL- and HE-LHC sensitivity to heavy colored particles \(\mathcal {Q}\)

Fig. 2
figure 2

The \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section as a function of minimum \(E_{\mathrm {T}}^{\mathrm {miss}}\) after the experimental selection criteria at 13 TeV, for \(\mathcal {Q}=T_p,~{\tilde{g}}\), and \({\tilde{t}}\), with current 95% CL\(_s\) limits after \(\mathcal {L}=36.1\) fb\(^{-1}\) of data at the 13 TeV LHC. Also shown is the estimated sensitivity with \(\mathcal {L}=3\) ab\(^{-1}\), assuming \(\delta _\mathrm{Syst.}=\)2.5 and 1% systematical errors                                             

We now compare and apply the model-independent upper limits on the cross sections derived in the previous section to the production of \(\mathcal {Q}\overline{\mathcal {Q}}\) pairs at the LHC. In Fig. 2, we show the model-independent 95% CL upper limits along with NLO+PS-accurate cross section for \(\mathcal {Q}\) produced in association with a hard jet at the matrix-element level, viz. \(\mathcal {Q}\overline{\mathcal {Q}}+j\) at NLO. (For details of our computational setup, see Sect. 2.) We overlay FO scale uncertainty, computed according to Eq. (6). For all cases of \(\mathcal {Q}\), we find that the scale uncertainty at NLO is a dominant source of uncertainty and hence take it as a representative measure of the total uncertainty. A dedicated and in-depth discussion of this uncertainty is given in Sect. 4.

From inclusive monojet searches at \(\sqrt{s} = 13\) TeV with \(\mathcal {L}=36.1 \mathrm ~fb^{-1}\) of data, we find that the lower limits on \(\mathcal {Q}\) masses stand at around \(m_{T_p}=\) 400 GeV for the fermionic top partner and \(m_{{\tilde{g}}}=\) 600 GeV for the gluino, while no constraint on stop masses is found within the range under consideration. We observe that in the high-\(p_T\) bins both the systematic and statistical experimental uncertainties play a crucial role. As argued in the previous section, one expects sensitivity to improve at the HL-LHC due to a much larger dataset, leading to better control on both uncertainties. From the scaled limits, we find that fermionic top partners with masses \(m_{T_p}\lesssim \) 800 GeV, gluinos with \(m_{{\tilde{g}}}\lesssim \) 1000 GeV, and stops with masses \(m_{{\tilde{t}}}\lesssim \) 600 GeV, in a compressed spectrum scenario, can be excluded at 13 TeV with \(\mathcal {L}=3\) ab\(^{-1}\), using the inclusive monojet signature.

Along with more data, a possibility that can greatly push the sensitivity to heavy colored particles is increasing the beam energy of the LHC itself. Presently, community discussions are underway on upgrading the LHC’s magnet system to handle a center-of-mass energy up to \(\sqrt{s}=27\) TeV [81]. In light of this prospect, we briefly investigate the impact of a higher beam energy on the production of \(\mathcal {Q}\) and the SM \(Z+\mathrm {jets}\) background, i.e., the dominant background of the monojet signature [18], and estimate the experimental reach of such a collider.

Table 4 The LO \(Z+j\) cross section (pb) for representative \(p_{T}^{\mathrm {gen}}\) (GeV) at \(\sqrt{s}=13\), 14, and 27 TeV

Before this, however, we shortly digress to describe our modeling of the SM \(Z\,+\,\)jets background. In particular, we note that while the signal processes are consistently determined at NLO in QCD with parton shower matching, we consider the \(Z+\)jets background only at LO with an experimentally determined normalization factor. For our purposes, we believe this provides a sufficiently reliable description of the SM backgrounds after selection cuts.

The motivation comes precisely from the fact that proper modeling of the SM EW boson+jets background for monojet searches is highly nontrivial [82], particularly in comparison to simulating the inclusive \(W/Z\,+\,\)jets process. The technical difficulty is due, in part, to strong phase space restrictions (cuts) on the final-state jets (see above Sect. 3.1 for the list of cuts). Present implementations of the MC@NLO formalism into general purpose event generators require that one integrates over the entire phase space of additional real radiation at NLO to ensure infrared pole cancellation. In the present case, this renders event generation at NLO inefficient. For example: the \(pp\rightarrow Z+j\) process with \(p_T^j \gtrsim 250-1000\) GeV is known [83, 84] to have giant QCD corrections stemming from the opening of new kinematic configurations. In this instance, large corrections are driven by the high-\(p_T\) dijet process with a relatively soft Z emission off a final-state quark, a configuration that would otherwise fail the \(\varDelta \phi (jet,E_{\mathrm {T}}^{\mathrm {miss}})\) and minimum \(E_{\mathrm {T}}^{\mathrm {miss}}\) selection criteria. Moreover, after careful consideration, we find that the EW boson + jets background, that survives the selection analysis, is dominated by Born-like configurations. For example: after cuts, partonic channels such as \(qg \rightarrow qZ\) contribute much more to the \(pp\rightarrow Z+\)jets background than channels like \(qq\rightarrow qqZ\).

Once selection criteria have been applied, the difference then between our background modeling and a much more precise determination, e.g., Ref. [82], is largely an overall normalization. This does not necessarily hold true when taking into account EW corrections at NLO and beyond. However, such corrections are beyond the claimed accuracy of our work and we refer readers to Refs. [83, 84] for thorough discussions. For representative minimum \(E_{\mathrm {T}}^{\mathrm {miss}}\) choices, we have checked at 14 TeV that the ratio \((K^\mathrm{NLO})\) of the \(pp\rightarrow Z+j\) rates at NLO+PS and LO+PS after selection cuts are applied is roughly a constant \(K^\mathrm{NLO} \sim 1.2\). This is consistent with the size of finite virtual corrections at NLO in QCD to the DY process [85], and supports our arguments that the residual EW+jets background exhibits Born-like kinematics. Now, to achieve a reliable normalization of the EW+jets background, we scale the generator-level \(pp\rightarrow Z+j\) cross section at LO by (approximately) a factor 1/20 so that the ratio of the curve to the signal cross section is normalized with respect to the post-event selection limit in Fig. 2. For \(\sqrt{s}=14\) and 27 TeV, this is additionally scaled by the cross section ratios \(\sigma (\sqrt{s}=14\mathrm{~TeV})/\sigma (\sqrt{s}=13\mathrm{~TeV})\) and \(\sigma (\sqrt{s}=27\mathrm{~TeV})/\sigma (\sqrt{s}=13\mathrm{~TeV})\), respectively, to account for the increase in parton luminosity. In Table 4, we list the cross sections for the \(Z+j\) process at \(\sqrt{s}=13,~14,\) and 27 TeV for representative \(p_{T}^{\mathrm {gen}}\) selections. We stress that this procedure is only an estimation of the SM background; we do not advocate that this is a suitable replacement for full NLO+PS (or more accurate) computations in experimental searches.

Fig. 3
figure 3

Upper: Cross sections of \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) process as a function of \(p_{T}^{\mathrm {gen}}\) for various \(\mathcal {Q}\) at \(\sqrt{s}=14\)  and 27 TeV. Lower: The ratio of the cross sections at 27 and 14 TeV

For representative masses, \(m_{\mathcal {Q}}\), we show in Fig. 3 the cross sections for the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) and (normalized) \(Z+j\) processes as a function of the leading jet \(p_T\). More specifically, the cross sections are calculated as a function of a generator-level \(p_T\) threshold \((p_{T}^{\mathrm {gen}})\) on the light jet. In the lower panel, we show the 27 TeV-to-14 TeV cross section ratios. There, one sees that the production cross section of \(\mathcal {Q}\) increases faster than the SM background with increasing center-of-mass energy. The enhancement follows from the well-documented [86,87,88] growth in PDF luminosities for fixed partonic mass scales but increasing collider beam energy. Quantitatively, for \(p_{T}^{\mathrm {gen}}=1\) TeV, the \(Z+j\) cross section increases by a factor of 10 with respect to the change of \(\sqrt{s}\), while the fermionic top partner cross section for \(m_{T_p} = 600~(800)\) GeV increases by approximately \(24\times ~(28\times )\), the gluino rate for \(m_{{\tilde{g}}}=900\) GeV by \(40\times \), and the stop rate for \(m_{{\tilde{t}}}=400\) GeV by \(20\times \). Although other sources of SM backgrounds for the monojet signature have been presently neglected, the signal over background ratio (S / N) still increases significantly at higher collider energies due to the larger luminosity enhancement for the signal process than dominant SM backgrounds. Subsequently, the HE-LHC enables ones to investigate parameter regions that are not accessible at the LHC. We do emphasize, however, that S / N ratios can change drastically if additional information is provided to enhance the separation of the signal events from the backgrounds. For example: proposals exist on how to utilize soft leptons, jets, and displaced vertices associated with decays of \(\mathcal {Q}\) that can further reduce SM background rates [89,90,91].

As shown in Fig. 4, once the scaled limits for the SM backgrounds are determined, we can compare the predicted cross sections for \(p p \rightarrow \mathcal {Q}\overline{\mathcal {Q}}+ j\) NLO process and estimate the expected reach at the 27 TeV HE-LHC. We find that with \(\mathcal {L}=3-15\) ab\(^{-1}\), one is sensitive to compressed spectra scenarios featuring fermionic top partners with masses \(m_{T_p}\lesssim 1100\) GeV, gluinos with masses \(m_{{\tilde{g}}}\lesssim 1800\) GeV, and stops with masses \(m_{{\tilde{t}}}\lesssim 600\) GeV.

We end this discussion by providing an estimate of the anticipated sensitivity of new colored particles \(\mathcal {Q}\) at the 27 TeV HE-LHC. Here, we assume that the SM background is dominated by the \(Z+j\) process and simply scale the model-independent 95% CL upper limit at \(\sqrt{s}=13\) TeV by the 27-to-13 TeV production cross section ratio. In other words, the SM background cross section at a given \(\sqrt{s}\) and \(E_{\mathrm {T},\mathrm {cut}}^{\mathrm {miss}}\), denoted as \(\sigma ({\sqrt{s}}; {E_{\mathrm {T},\mathrm {cut}}^{\mathrm {miss}}}) \), is obtained from the relation

$$\begin{aligned}&\sigma ({\sqrt{s}}; {E_{\mathrm {T},\mathrm {cut}}^{\mathrm {miss}}}') = \sigma ({13\;\mathrm {TeV}}; E_{\mathrm {T},\mathrm {cut}}^{\mathrm {miss}}) \nonumber \\&\quad \cdot \frac{\sigma _{pp \rightarrow Zj}(\sqrt{s};p_{T}^{\mathrm {gen}}= {E_{\mathrm {T},\mathrm {cut}}^{\mathrm {miss}}}' )}{\sigma _{pp\rightarrow Zj}({13\;\mathrm {TeV}};p_{T}^{\mathrm {gen}}= E_{\mathrm {T},\mathrm {cut}}^{\mathrm {miss}})}. \end{aligned}$$
(12)

In the above, \(\sigma _{pp \rightarrow Zj}(\sqrt{s};p_{T}^{\mathrm {gen}})\) is the LO \(Z+j\) cross section with \(p_{T}^{\mathrm {gen}}\) at a collider energy \(\sqrt{s}\). We further assume that the detector acceptance and efficiencies are the same at 13 and 27 TeV. This assumption is not as strong as one may anticipate in more general circumstances. The HE-LHC project proposes to refit, replace, and/or upgrade the current LHC magnet system and detector experiments. As the detector experiment caverns themselves cannot physically grow, one is forced to adopt a detector fiducial volume at 27 TeV that is largely unchanged from 13 TeV. Similarly, we also assume systematic uncertainties of 2.5% and 1%, the same considered in Sect. 3.1.

Fig. 4
figure 4

Same as Fig. 2 but scaled limits for \(\sqrt{s}=27\) TeV assuming a \(\mathcal {L}=3\mathrm ~ab^{-1}\) and b \(15\mathrm ~ab^{-1}\)

3.3 Properties determination of heavy colored particles

We now turn to the possibility of extracting properties of the heavy colored particle \(\mathcal {Q}\) from jet behavior within the monojet signature. As briefly discussed in the introduction, asserting color representation and spin of \(\mathcal {Q}\) is required to infer information on its mass from cross section measurements (or limits). Consequently, a single cross section measurement of a particular monojet signal region does not help much in determining the nature of \(\mathcal {Q}\). For example: in Fig. 3, one sees that the production cross section for the process \(p p \rightarrow \mathcal {Q}\overline{\mathcal {Q}}+ j\) with \(p_{T}^{\mathrm {gen}}=600\) GeV, for a top squark \({\tilde{t}}\) with \(m_{{\tilde{t}}} = \) 400 GeV, a gluino \({\tilde{g}}\) with \(m_{{\tilde{g}}} = \) 900 GeV, and a fermionic top partner \(T_p\) with \(m_{T_p} = \) 600 GeV are roughly within \(\varDelta \sigma \sim \) 5–10 fb of one another. However, despite this ambiguity, it is still possible to extract information from the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+ j\) cross section as a function of the leading jet \(p_T\), which can be measured directly, since it obeys a distinguishing pattern for each \(\mathcal {Q}\) hypothesis. That the nature of \(\mathcal {Q}\) is, in part, encoded in this observable reflects a nontrivial interplay between \(\mathcal {Q}\)’s mass, \(m_\mathcal {Q}\), its color representation and spin, and the dimensionless ratio \((m_\mathcal {Q}/p_{T}^{j_1})\). This interplay is what we now discuss.

The first discerning observation is that the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross sections do not depend on \(p_{T}^{\mathrm {gen}}\) in a universal manner. Keeping to Fig. 3, one sees that while \(\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; p_{T}^{\mathrm {gen}}=600\) GeV) are the same at \(\sqrt{s}=14\) TeV for the \((\mathcal {Q},m_{\mathcal {Q}})\) configurations under consideration, the relative size of \(\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j)\) changes with \(p_{T}^{\mathrm {gen}}\). In other words, while \(\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; p_{T}^{\mathrm {gen}})\) follows an anticipated power-law of \(\sigma (p_{T}^{\mathrm {gen}})\sim (p_{T}^{\mathrm {gen}})^{-\beta }\), with \(\beta >0\), the precise value of the exponent is dependent on the color and spin structure of \(\mathcal {Q}\). In a particular extreme, the gluino rate is the smallest (largest) of the configurations for \(p_{T}^{\mathrm {gen}}\) smaller (larger) than \(p_{T}^{\mathrm {gen}}=600\)  GeV, suggesting a smaller \(\beta \) than for other \(\mathcal {Q}\). Information on \(\mathcal {Q}\) can be extracted from \(\sigma (p_{T}^{\mathrm {gen}})\) by considering its ratio with respect to a benchmark \(\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j;p_{T}^{\mathrm {gen}})\). For example: for the benchmark process \({\tilde{t}}{\tilde{t}}^*+j\) with \(m_{{\tilde{t}}}=400\mathrm{~GeV}\), the \(\sqrt{s} = \) 14 TeV cross section ratios at \(p_{T}^{\mathrm {gen}}\) = 400 and 800 GeV are \(\sigma ({\tilde{g}}; m_{{\tilde{g}}}=900\mathrm{~GeV})/\sigma ({\tilde{t}})\) = 0.79 and 1.33, respectively, and similarly \(\sigma (T_p;m_{T_p}=600\mathrm{~GeV})/\sigma ({\tilde{t}})\) = 0.85 and 1.23. From this one can determine that the change in the cross section ratios over the range 400 GeV \(<p_{T}^{\mathrm {gen}}<\) 800 GeV, is \(\varDelta (\sigma /\sigma _\mathrm{Ref.})/\varDelta p_{T}^{j_1}\ \sim \mathcal {O}(14\%/100\mathrm{~GeV})\) for gluinos and \(\mathcal {O}(9.5\%/100\mathrm{~GeV})\) for fermionic top partners. Hence, cross section ratios for two different \(\mathcal {Q}\) hypotheses is crucially dependent on the choice of \(p_{T,\mathrm {cut}}^{j_1}\).

In addition, one can also think about the physics case where there are multiple heavy particles making up the total cross section consistent with a lighter particle. Take for example that the ratio of \(\sigma (T_p\bar{T}_p+j; m_{T_p}=800\mathrm{~GeV})/\sigma (T_p\bar{T}_p+j;m_{T_p}=600\mathrm{~GeV})\) is 0.21 at \(p_{T}^{\mathrm {gen}}= 400\mathrm{~GeV}\). Hence, five copies of \(T_p\) with mass \(m_{T_p}=800\mathrm{~GeV}\) can mimic the cross section of a single \(T_p\) with \(m_{T_p}=600\mathrm{~GeV}\) at \(p_{T}^{\mathrm {gen}}=400\mathrm{~GeV}.\) At the \(p_{T}^{\mathrm {gen}}=800\mathrm{~GeV}\), however, \(5\times \sigma (T_p\bar{T}_p+j; m_{T_p}=800\mathrm{~GeV})\) / \(\sigma (T_p\bar{T}_p+j, m_{T_p}=600\mathrm{~GeV})= 1.25.\) The one- and five-copy scenarios then predict  \(\varDelta (\sigma /\sigma _\mathrm{Ref.})/\varDelta p_{T}^{j_1}\) = \( \mathcal {O}(1\%/100\mathrm{~GeV})\) and \(\varDelta (\sigma /\sigma _\mathrm{Ref.})/\varDelta p_{T}^{j_1}\) = \( \mathcal {O}(6\%/100\mathrm{~GeV})\), respectively, thus providing a means to check this potential degeneracy. Likewise, four copies of \(T_p\) with \(m_{T_p}=800\mathrm{~GeV}\) can mimic the the cross section of \(m_{T_p}=600\mathrm{~GeV}\) at \(p_{T}^{\mathrm {gen}}= 800\mathrm{~GeV}.\) At \(p_{T}^{\mathrm {gen}}=400\mathrm{~GeV},\) however, one finds that \(4\times \sigma (m_{T_p}=800\mathrm{~GeV})/\sigma (m_{T_p}=600\mathrm{~GeV})= 0.84\). For the one- and four-copy scenarios, this leads to the predictions of \(\varDelta (\sigma /\sigma _\mathrm{Ref.})/\varDelta p_{T}^{j_1}= \mathcal {O}(1\%/100\mathrm{~GeV})\) and \(\varDelta (\sigma /\sigma _\mathrm{Ref.})/\varDelta p_{T}^{j_1}= \mathcal {O}(4\%/100\mathrm{~GeV})\), respectively.

The different dependence on \(p_{T}^{\mathrm {gen}}\) observed for gluinos, fermionic top partners, and stops arise from the fact that heavier particles give rise to harder, i.e. less steeply falling, \(p_T^{j_1}\) distributions. The benchmark cases having significantly different input masses, with \(m_{{\tilde{g}}}=900\mathrm{~GeV}\), \(m_{T_p} = 600\) and \(800\mathrm{~GeV}\), and \(m_{{\tilde{t}}}=400\mathrm{~GeV}\), and makes a significant impact on the \(p_{T}^{j_1}\) dependence. To isolate this behavior, in Fig. 5 we plot the cross section double ratio as a function of \(p_{T,\mathrm {cut}}^{j_1}\),

$$\begin{aligned}&P_{\mathcal {Q}}(p_{T,\mathrm {cut}}^{j_1},m_{\mathcal {Q}}) \quad = \frac{\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; m_{\mathcal {Q}}, p_{T,\mathrm {cut}}^{j_1}) / \sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; m_{\mathcal {Q}}, p_{T,\mathrm {cut}}^{j_1}=p_{T}^{*})}{\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; m_{\mathcal {Q}}=m_{*}, p_{T,\mathrm {cut}}^{j_1}) / \sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; m_{\mathcal {Q}}=m_{*}, p_{T,\mathrm {cut}}^{j_1}=p_{T}^{*})}. \end{aligned}$$
(13)

In the top (bottom) single ratio, both cross sections are with respect to the mass \(m_{\mathcal {Q}}~(m_*)\) but different \(p_{T,\mathrm {cut}}^{j_1}\). This has the effect of canceling overall color and kinematic factors while isolating logarithmic terms of the form \(\log (m_\mathcal {Q}/p_{T,\mathrm {cut}}^{j_1})^2\). The double ratio, then, is a measure of this logarithmic dependence with respect to a baseline mass \(m_*\) and minimum transverse momentum \(p_T^*\). For \(\mathcal {Q}=T_p\), we choose \(m_*=600\mathrm{~GeV}\) and \(p_T^*=200\mathrm{~GeV}\), and plot \(P_{\mathcal {Q}}(p_{T,\mathrm {cut}}^{j_1})\) for \(m_{T_p}=400\) and \(800\mathrm{~GeV}\). Quantitatively, one sees that the double ratio increases (decreases) by about 50% at \(p_{T,\mathrm {cut}}^{j_1}=800\) GeV for \(m_{T_p} =800~(400)\) GeV. This feature is universal for particles in the same color representation and follows from the nature of massless gauge boson emission in scattering processes.

To better understand this behavior, consider the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) process. The t-channel propagators gives rise, after phase space integration, to the aforementioned logarithms \(\log (m_\mathcal {Q}/p_{T,\mathrm {cut}}^{j_1})^2\). In the context of parton shower resummation, this dependence can be interpreted as the likelihood of emitting an additional QCD parton with transverse momentum \(p_{T,\mathrm {cut}}^{j_1}\), i.e., the differential probability is proportional to \(d\mathcal {P} = (1/\sigma ) d\sigma \propto \alpha _s(m_\mathcal {Q})\log (m_{\mathcal {Q}}^2/p_T^{j 2})\). Hence, a fixed probability \(\varDelta \mathcal {P}\) implies a fixed \((m_{\mathcal {Q}}/p_T^{j_1})^2\) ratio, and indicates that increasing \(m_{\mathcal {Q}}\) results in \(p_T^j\) increasing commensurately. Qualitatively, the emission of higher-\(p_T\) QCD partons becomes easier for heavier \(\mathcal {Q}\) because high-\(p_T\) emissions become relatively soft as \(m_{\mathcal {Q}}\) increases. This results in a rightward shift of the so-called Sudakov shoulder [92, 93]. For TeV-scale particles, the rightward shift of what constitutes “soft” is known to be large; see for example Refs. [94, 95].

The discrimination power of cross section ratios extends if one considers the additional dependency on a collider’s beam energy. In particular, we find that the cross section ratios shift from unity for \(p_{T}^{\mathrm {gen}}=600\) GeV at \(\sqrt{s}=14\) TeV, to \(\sigma ({\tilde{g}})/\sigma ({\tilde{t}}) = 2.1\) and \(\sigma (T_p) /\sigma ({\tilde{t}})=\) 1.35 at 27 TeV. Hence, a measurement of the ratio of the cross section at the \(\sqrt{s}=14\) and 27 TeV with 30% accuracy can resolve our benchmark gluino, fermionic top parter, and stop scenario. Moreover, increasing the beam energy can also significantly improve the signal-over-background ratio, thereby enabling measurements of the \(\mathcal {Q}\overline{\mathcal {Q}}+j\) cross section for different minimum \(p_T^j\) over a wide range of \(p_T\). Measuring the signal cross section for several \(p_{T,\mathrm {cut}}^{j_1}\) at \(\sqrt{s}=27\mathrm{~TeV}\) with 10% accuracy allows one to distinguish the benchmark scenario by narrowing down the mass of \(\mathcal {Q}\) without \(\sqrt{s}= 14\mathrm{~TeV}\) information. In light of this, it is necessary to emphasize that theoretical predictions on total and differential cross sections, as well as their ratios, must have the requisite accuracy to make these measurements.

Lastly, we note that the slope of \(\sigma (p_{T}^{j_1})\) with respect to \(p_{T}^{j_1}\) does not significantly depend on either the color structure or spin of \(\mathcal {Q}\). To show this, we display in Fig. 6 a second cross section double ratio,

$$\begin{aligned}&R_{\mathcal {Q}}(p_{T}^{\mathrm {gen}},m_{\mathcal {Q}}) = \frac{\sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; m_{\mathcal {Q}}, p_{T}^{\mathrm {gen}}) / \sigma (\mathcal {Q}\overline{\mathcal {Q}}+j; m_{\mathcal {Q}}, p_{T}^{\mathrm {gen}}=p_{T}^{*})}{\sigma ({\tilde{t}}{\tilde{t}}^*+j; m_{{\tilde{t}}}=m_{\mathcal {Q}}, p_{T}^{\mathrm {gen}}) / \sigma ({\tilde{t}}{\tilde{t}}^*+j; m_{{\tilde{t}}}=m_{\mathcal {Q}}, p_{T}^{\mathrm {gen}}=p_{T}^{*})}. \end{aligned}$$
(14)

The structure of \(R_{\mathcal {Q}}(p_{T}^{\mathrm {gen}},m_{\mathcal {Q}})\) here is analogous to \(P_{\mathcal {Q}}(p_{T,\mathrm {cut}}^{j_1},m_{\mathcal {Q}})\) in Eq. (13), but differs in that the normalizing process is fixed to \(pp\rightarrow {\tilde{t}}{\tilde{t}}^*+j\) with \(p_{T}^{\mathrm {gen}}=p^{*}_{T}=300\) GeV for \(\sqrt{s}=14\) TeV, and we vary \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) in the upper ratio. Subsequently, while overall color and kinematic multiplicative factors cancel, the relative dependence of \(q\overset{(-)}{q}\), \(\overset{(-)}{q}g\), and gg scattering within an individual process does not cancel and is inherently dependent on the color representation.

From Fig. 6, one can observe that fermionic top partners and stops with same mass have the almost same slope over a wide range of \(p_T\). On the other hand, because of color and matrix element effects from gluinos, the ratio \(R_{\tilde{g}}({p^{gen}_T})\) tend to decrease with \(p_T\) but again independent to the mass of the heavy particle. For validation, we consider the case for scalar gluons, known as sgluons [35,36,37], and overlay sgluon behavior on the same figure.Footnote 3 As the plot suggests, the distribution of the hard radiation associated with the colored particle production follows universally for a given color representation of \(\mathcal {Q}\).

4 Theoretical uncertainties of the monojet process

As investigated in Sect. 3.3, were one to discover new colored particles at the LHC, or a potential successor experiment, a measurement of the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section as a function of the leading jet’s minimum transverse momentum \((p_{T,\mathrm {cut}}^{j_1})\) can help establish the quantum numbers of \(\mathcal {Q}\). Ascertaining such information, however, requires accurate BSM signal predictions. For the benchmark scenarios listed in Table 1, we find one needs theoretical uncertainties no larger than \(\mathcal {O}(30\%)\) and \(\mathcal {O}(5\%)\), respectively, on the total inclusive cross section normalizations and on the change of the cross section for \(\sigma (pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j)\) as a function of \(p_{T,\mathrm {cut}}^{j_1}\) per 100 GeV. In this section, we discuss and quantify theoretical uncertainties associated with the monojet signal process. We particularly investigate (potential) sources of uncertainties when employing the state-of-the-art MC suites MG5_aMC@NLO +PY8 and Sherpa.

Fig. 5
figure 5

The double cross section ratio \(P_{T_p}(p_{T,\mathrm {cut}}^{j_1},m_{T_p})\), as defined in Eq. (13), as a function of \(p_{T,\mathrm {cut}}^{j_1}\) for fermionic top partners with \(m_{T_p} = 400\) and \(800\mathrm{~GeV}\), and normalization set at \((m_*,p_T^*)=(600\mathrm{~GeV},200\mathrm{~GeV})\)

Fig. 6
figure 6

The cross section ratio \(\sigma ( p p \rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j)/\sigma ( p p \rightarrow {\tilde{t}}{\tilde{t}}^*+j)\), both at NLO in QCD, as a function of a common \(p_{T}^{\mathrm {gen}}\), assuming \(m_\mathcal {Q}= m_{{\tilde{t}}}\) at \(\sqrt{s}=14\) TeV. The ratio is normalized to \(R(p_{T}^{\mathrm {gen}})=1\) at \(p_{T}^{\mathrm {gen}}=300\) GeV

To investigate theoretical uncertainties associated with the \(\mathcal {Q}\overline{\mathcal {Q}}+j\) process, we focus on fermionic top partners \(T_p\). We choose this benchmark because of wide implementability across different event generators as well as its comparability to \(t\bar{t}+j\) production in the SM. We compare several simulation setups and techniques at different levels of precision within QCD; see Sect. 2 for details. As the relevant signal topology in this work is characterized by the presence of hard QCD radiation recoiling against the \(T_p\bar{T}_p\) system, our primary benchmark observable is the \(pp\rightarrow T_p\bar{T}_p+j\) cross section as a function of the transverse momentum of the process’ leading, i.e., highest \(p_{T}\), jet \((p_{T}^{j_1})\).

Fig. 7
figure 7

Normalized scale uncertainty bands of the \(pp\rightarrow T_p\bar{T}_p+j\) cross section as a function of minimum \(p_{T}^{j_1}~(p_{T,\mathrm {cut}}^{j_1})\), for a the inclusive \(pp\rightarrow T_p\bar{T}_p\) process at NLO+PS (dash) and at LO+PS as calculated with up to two additional matrix element-level jets via multijet merging (dot); as well as b the inclusive \(pp\rightarrow T_p\bar{T}_p\) process at NLO+PS (dot) and the \(pp\rightarrow T_p\bar{T}_p+j\) process (no multijet merging) at NLO+PS (dash)

We begin with Fig. 7a, where we show the \(pp\rightarrow T_p\bar{T}_p+j\) cross section as a function of \(p_{T}^{j_1}\) as derived from the fully inclusive \(pp\rightarrow T_p\bar{T}_p\) calculation at NLO+PS (dash), as obtained from the MC@NLO formalism. We also show the calculation at LO+PS (dot), as obtained with up to two additional matrix element-level jets via multijet merging. The curves are shown with their factorization \((Q_F)\) and renormalization \((Q_R)\) variation envelopes and are normalized to the respective nominal prediction. As discussed in Sect. 2, the leading jet \(p_{T}\) distribution for both calculations is only LO accurate for \(p_{T}^j\) comparable to the scale of the hard process. Hence, the scale uncertainties for \(p_{T}^j\gtrsim Q_F,Q_R\) are large, asymmetric, equal for the two calculations, and span approximately \(+\,45\) to \(-\,30\%\). For smaller \(p_{T}^j\), namely \(p_{T}^j\lesssim m_{T_p}=600\mathrm{~GeV}\), one observes that the two uncertainty envelopes begin to differ. Whereas the uncertainty for the multijet calculation reduce only slightly for decreasing \(p_{T}^j\), the MC@NLO uncertainty reduces to roughly the \(\pm \,10\%\) level. The difference originates from virtual corrections present in the NLO+PS calculation, which soften dependencies on \(Q_F,Q_R\), but are obviously absent in the LO+PS calculation. The Sudakov-like factor in the multijet merging prescription only partially reduces the dependence on \(Q_F\) by matching low-\(p_{T}\) QCD emissions in the hard matrix element with those in the PDF. Hence, even for observables that are formally of the same precision, the presence of all \(\mathcal {O}(\alpha _s)\) terms in the NLO+PS calculation leads to a smaller scale dependence at low \(p_{T}^{j_1}\) than in the merged LO+PS calculation.

The small scale uncertainty observed for the lowest \(p_{T}^{j_1}\) suggests high theoretical precision is achievable with NLO+PS computations obtained via the MC@NLO prescription. However, unlike pure FO calculations, calculations matched to parton showers possess the additional dependence on the parton shower starting scale \(Q_S\). Within the MC@NLO formalism, \(Q_S\) controls whether the leading \(\mathcal {O}(\alpha _s)\) emission beyond the Born process is included in the FO matrix element or the all-orders parton shower. Loosely speaking, QCD radiations with \(p_{T}^j\) above (below) \(Q_S\) originate from the hard matrix element (parton shower). As pointed out in Refs. [96, 97], lowest order-accurate observables, e.g., the \(p_{T}^{j_1}\) distribution when the \(pp\rightarrow T_p\bar{T}_p\) process is evaluated at NLO+PS, and processes that possess large virtual corrections suffer from ambiguities when choosing \(Q_S\). This manifests as a strong dependence on \(Q_S\), and hence a large deviation.

Fig. 8
figure 8

The \(pp\rightarrow T_p\bar{T}_p+j\) cross section, derived from the inclusive \(pp\rightarrow T_p\bar{T}_p\) calculation at NLO+PS, as a function of \(p_{T,\mathrm {cut}}^{j_1}\), assuming multiplicative variations of the default parton shower scale. Rates are normalized with respect to the fNLO prediction

We assess this behavior in Fig. 8 by plotting the \(pp\rightarrow T_p\bar{T}_p+j\) cross section, derived from the inclusive \(pp\rightarrow T_p\bar{T}_p\) calculation at NLO+PS, as a function of \(p_{T,\mathrm {cut}}^{j_1}\). We assume multiplicative variations of the default parton shower scale, given in Eq. (7), but fix the factorization and renormalization scales to their central values. Rates are normalized to the fixed order NLO (fNLO) prediction. The \(Q_S\) variation of the NLO+PS result with respect to the default \(Q_S\) choice spans roughly \(\pm \,25\%\), with the dependence increasing (decreasing) for larger (smaller) values of \(Q_S\) over the range of \(p_{T}^{j_1}\) considered. The shower scale variation amounts to absolute deviations from the fNLO result up to about \(+\,60\) and \(-\,35\%\). For vanishing \(p_{T}^j\), one should take caution in interpreting the vanishing shower scale uncertainty. In this limit, the FO calculation is unphysical. The FO calculation possess an integrable singularity at \(p_{T}^j=0\) that leads to arbitrarily large cross sections. As a result, the ratio of the NLO+PS-level cross section, a finite and physical quantity, to the fNLO cross section, the unphysical quantity, vanishes as \(p_{T}^j/m_{T_p\bar{T}_p}\rightarrow 0\). For the curves with shower scale multiplier 0.5 and 1, we observe the curves are flat beyond \(\sim 300\) GeV and 600 GeV, respectively, as expected. However, the curves with shower scale multiplier larger than 1, we observe a crossing at \(p_T \sim 700\) GeV. The large shower starting scale introduces high \(p_T\) initial state radiations, which is neither soft nor collinear to the incoming parton. If the high \(p_T\) jets generated from the parton shower of Born-type S events immersed in this region, then the cross section is not reliable. Therefore, the curves (and the crossing) are not trustworthy beyond the top partner mass, or equivalently \(\sqrt{\hat{s}}/2\). A better way to estimate the \(pp\rightarrow T_p \bar{T}_p + j\) cross section (and associated shower scale uncertainty) would be to consider the process \(pp\rightarrow T_p \bar{T}_p+j\) at NLO+PS.

In light of the large theoretical uncertainties in the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section stemming from \(Q_F,~Q_R,\) and \(Q_S\), it is clear that the precision achieved with the aforementioned methods is insufficient for distinguishing \(\mathcal {Q}\) candidates. To explore if such precision is yet still possible with presently available general-purpose Monte Carlo technology, we consider the cross section obtained from the \(pp\rightarrow T_p\bar{T}_p+j\) process itself at NLO+PS. In Fig. 7b we plot the normalized \(pp\rightarrow T_p\bar{T}_p+j\) cross section with its \(Q_F,~Q_R\) envelope again as a function of \(p_{T,\mathrm {cut}}^{j_1}\), but as derived from the \(pp\rightarrow T_p\bar{T}_p\) (dot) and \(pp\rightarrow T_p\bar{T}_p+j\) (dash) calculations at NLO+PS. We observe that the uncertainty reduces to a largely uniform band that is just shy of \(\pm \,20\%\) for the \(T_p\bar{T}_p+j\) calculation. Both uncertainties are comparable for \( p_{T}^{j_1}\ll m_{T_p}\); for \(p_{T}^{j_1}\gg m_{T_p}\), however, the \(T_p\bar{T}_p+j\) uncertainty is about 40–60% smaller.

It is worth reiterating that the FO \(T_p\bar{T}_p+j\) calculation requires a minimum \(p_{T}^j\) selection to regulate matrix element poles and, since no Sudakov suppression is applied to low \(p_{T}\) QCD emissions in this case, to ensure perturbative stability, in the Collins-Soper-Sterman (CSS) sense [92]. We safeguard against the need for \(k_T\)-resummation on the leading jet by using the (somewhat conservative) CSS consistency requirement from Ref. [95]. For \(m_{T_p}=600\mathrm{~GeV}\), and hence \(m_{T_p\bar{T}_p}>1.2\mathrm{~TeV}\), one needs \(p_{T,\mathrm {cut}}^{j_1}\gtrsim 200-250\mathrm{~GeV}\) to force collinear logarithms from t-channel exchanges, which scale as \(\delta \sigma \sim \alpha _s\log (m_{T_p\bar{T}_p}/p_{T}^{j})^2\), to be much smaller than 1.

There is now an apparent conflict between theoretical needs and computational capabilities: While \(pp\rightarrow T_p\bar{T}_p+j\) at NLO+PS provides improved control and stability over \(Q_F\) and \(Q_R\), the calculation is only meaningful for sufficiently large \(p_{T}^{j_1}\). The \(pp\rightarrow T_p\bar{T}_p\) process at NLO+PS, on the other hand, extends the support for \(p_{T}^{j_1}\) to low \(p_{T}\) but at the cost of a larger uncertainties. Such demands, however, are precisely resolved with multijet merging at NLO+PS, i.e., MEPS@NLO. As the CKKW MEPS@NLO prescription is natively available in the event generator Sherpa, for the remainder of this section we use Sherpa to further quantify scale uncertainties in \(T_p\bar{T}_p\) production. However, before proceeding, we non-trivially demonstrate that at least up to the QCD order presently being investigated, both MG5_aMC@NLO +PY8 and Sherpa  predictions are in agreement with one another.

Fig. 9
figure 9

Comparison of the \(pp\rightarrow T_p\bar{T}_p+j\) cross section estimated using MG5_aMC@NLO +PY8 and Sherpa  with scale uncertainty band. a Comparison of the \(pp\rightarrow T_p\bar{T}_p+j\) cross section with \(Q_F,~Q_R\) uncertainty band as derived from the \(pp\rightarrow T_p\bar{T}_p+j\) process at NLO+PS (MG5_aMC@NLO +PY8) and from the \(pp\rightarrow T_p\bar{T}_p\) process with MEPS@NLO multijet merging to one additional jet (Sherpa). Curves are individually normalized to their central value. b The same quantities as a but overlaid with the individual \(pp\rightarrow T_p\bar{T}_p+j\) at NLO+PS samples (Sherpa) that constitute the MEPS@NLO result. Unmerged NLO+PS calculations are distinguished by their \(p_{T}^{\mathrm {gen}}\) selection on \(p_{T}^{j_1}\)

In Fig. 9a, we plot the \(pp\rightarrow T_p\bar{T}_p+j\) cross section and its scale uncertainty band as a function of \(p_{T,\mathrm {cut}}^{j_1}\) as derived from the \(pp\rightarrow T_p\bar{T}_p+j\) process at NLO+PS using MG5_aMC@NLO +PY8, as well as from the \(pp\rightarrow T_p\bar{T}_p\) process at NLO+PS merged with up to one additional NLO-accurate, matrix element-level jet using Sherpa. The curves are normalized to the Sherpa rate. For \(p_{T,\mathrm {cut}}^{j_1}>200\mathrm{~GeV}\), one sees that the Sherpa result, much like the MG5_aMC@NLO +PY8 result, possesses a stable uncertainty envelope spanning roughly \(\pm \,18\%\). Moreover, the central value of Sherpa remains only 3–11% below MG5_aMC@NLO +PY8, and is consistent with the difference choice of \(\alpha _s(Q_s)\) employed by Sherpa. For \(p_{T}^{j_1}<200\mathrm{~GeV}\), the MEPS@NLO calculation reveals that the \(Q_F\) and \(Q_R\) dependence remains smooth, and tapers down to \(\pm \,16\%\) as \(p_{T,\mathrm {cut}}^{j_1}\) is relaxed. In comparison to the inclusive \(pp\rightarrow T_p\bar{T}_p\) at NLO+PS in Fig. 7a, the uncertainty band here is \(\mathcal {O}(5\%)\) larger in both directions. The slightly larger uncertainty here is consistent with the presence of additional soft \(\mathcal {O}(\alpha _s^2)\) radiation terms whose scale dependence would otherwise be stabilized by \(\mathcal {O}(\alpha _s^2)\) virtual terms.

As described in Sect. 2, to construct the MEPS@NLO sample, various \(pp\rightarrow T_p\bar{T}_p+j\) NLO+PS computations at increasing \(p_{T}^{\mathrm {gen}}\) are necessary to populate the \(p_{T}^j\) tail for \(p_{T}^j\gg p_{T}^{\mathrm {gen}}\). To further demonstrate the comparability of the MG5_aMC@NLO +PY8 and Sherpa curves, we show in Fig. 9b the same quantities plotted in Fig. 9a for \(p_{T,\mathrm {cut}}^{j_1}>200\mathrm{~GeV}\) but overlaid with the individual, unmerged \(pp\rightarrow T_p\bar{T}_p+j\) samples. For the range of \(p_{T,\mathrm {cut}}^{j_1}\) considered, we find that the difference between the MG5_aMC@NLO +PY8 NLO+PS prediction and the Sherpa NLO+PS predictions for \(p_{T,\mathrm {cut}}^{j_1}>p_{T}^{\mathrm {gen}}\) differ only by about 0–15% again within the \(Q_F,~Q_R\) uncertainty.

A theoretical difficulty of the inclusive monojet signature is ensuring perturbative control over predictions, despite complicated, multi-scale requirements on jet momenta and particle multiplicity. Hence, it is important to ensure that sub-leading jets in the monojet processes can be modeled as well as the leading jet. To investigate whether such uncertainties for subleading jets are attainable, we compare MEPS@NLO uncertainties for the \(pp\rightarrow T_p\bar{T}_p\) process when merged with different jet multiplicities. As a control, we show in Fig. 10a the (normalized) \(pp\rightarrow T_p\bar{T}_p+j\) cross section and uncertainty band as a function of \(p_{T,\mathrm {cut}}^{j_1}\) for the \(pp\rightarrow T_p\bar{T}_p\) process at NLO+PS merged with up to one (dark curve) or two (light curve) additional NLO-accurate, matrix element-level jets. No significant difference between the two curves is observed. One does not expect to see such deviations for this observable as the two calculations are identical at this order of \(\alpha _s\). The difference appear at one order higher: At the matrix element level, the two calculations possess up to two and three hard QCD emissions, respectively. In the \(1j\) MEPS@NLO calculation, the second jet is LO+LL (alternatively, LO+PS) accurate. In the \(2j\) MEPS@NLO calculation, the second jet is NLO+LL accurate and the third jet is LO+LL accurate. Therefore, in Fig. 10b, we plot the \(pp\rightarrow T_p\bar{T}_p+2j\) cross section as a function of minimum \(p_{T}^{j_2}\) for the two calculations. In analogy to the comparison of \(T_p\bar{T}_p\) and \(T_p\bar{T}_p+j\) production at NLO+PS in Fig. 7b, we observe here that the 1j MEPS@NLO prediction suffers from an uncertainty spanning roughly \(+\,45\) to \(-\,30\%\) whereas the \(2j\) MEPS@NLO rate exhibits a largely uniform uncertainty of about \(\pm \,22\%\). For vanishing \(p_{T,\mathrm {cut}}^{j_2}\), the uncertainties become equal since the observable is no longer sensitive to such a high order of \(\alpha _s\). With successive MEPS@NLO matching to higher jet multiplicities, one expects a comparable reduction of dependence on \(Q_F,~Q_R\) in the corresponding cross section. However, one should not expect the uncertainty to drop below the observed \(\mathcal {O}(10-15\%)\) without first accounting for missing two-loop terms at \(\mathcal {O}(\alpha _s^2)\).

Fig. 10
figure 10

Factorization and renormalization scale uncertainty bands for a the \(pp\rightarrow T_p\bar{T}_p+j\) cross section as a function of \(p_{T,\mathrm {cut}}^{j_1}\) and b the \(pp\rightarrow T_p\bar{T}_p+2j\) cross section as a function of \(p_{T,\mathrm {cut}}^{j_2}\) assuming \(p_{T,\mathrm {cut}}^{j_1}=300\mathrm{~GeV}\), for the \(pp\rightarrow T_p\bar{T}_p\) process with MEPS@NLO matching up to one (dark curve) or two (light curve) additional jets. Curves are normalized to the \(2j\) MEPS@NLO rate

Fig. 11
figure 11

Parton shower scale dependence for a the \(pp\rightarrow T_p\bar{T}_p+j\) cross section as a function of \(p_{T,\mathrm {cut}}^{j_1}\) and b the \(pp\rightarrow T_p\bar{T}_p+2j\) cross section as a function of \(p_{T,\mathrm {cut}}^{j_2}\) assuming \(p_{T,\mathrm {cut}}^{j_1}=300\mathrm{~GeV}\), for the \(pp\rightarrow T_p\bar{T}_p\) process with MEPS@NLO matching up to one (upper) or two (lower) additional jets. Curves are normalized to the respective rate assuming the default shower scale

We now turn to the issue of parton shower scale uncertainties for processes evaluated with MEPS@NLO. Following the same variation procedure as before, we plot in Fig. 11a, b, respectively, the \(Q_S\) dependence for the \(T_p\bar{T}_p+j\) and \(T_p\bar{T}_p+2j\) cross sections as a function of minimum \(p_T\) for the leading \((p_T^{j_1})\) and subleading \((p_T^{j_2})\) jet. The upper (lower) panel corresponds to the \(1j\)  (\(2j\)) MEPS@NLO prediction. All curves are normalized to the appropriate cross section evaluated at the default scale choice. We find for both calculations that the uncertainty in the \(T_p\bar{T}_p+j\) cross section is very small, reaching maximally to \(\pm \,4\%\) over the range of \(p_{T}^{j_1}\) considered, and comparable to our MC, i.e., statistical, uncertainty. Relative to the \(T_p\bar{T}_p\) at NLO+PS uncertainty in Fig. 8, which exhibits deviations up to \(\pm \,25\%\), this represents significant reduction in \(Q_S\) dependence. For the \(T_p\bar{T}_p+2j\) cross section, we see a larger scale dependence in both the \(1j\) and \(2j\) calculations. The uncertainties for the two rates are comparable, reaching about \(+\,10\) to \(-\,5\%\), and again comparable to our MC uncertainty.

Fig. 12
figure 12

The normalized a PDF uncertainty and b MEPS@NLO jet merging scale \(Q_{\mathrm {cut}}\) dependence for the \(pp\rightarrow T_p\bar{T}_p+j\) cross section as a function of \(p_{T,\mathrm {cut}}^{j_1}\), as derived from the \(1j\) MEPS@NLO calculation

Having established that multijet-merged simulations with NLO matrix elements feature moderate-to-small factorization, renormalization, and shower scale dependencies, we lastly investigate the uncertainty associated with our PDF input and the MEPS@NLO merging scale. In Fig. 12a, we show again the normalized \(pp\rightarrow T_p\bar{T}_p+j\) cross section, as derived from the \(T_p\bar{T}_p\) process with MEPS@NLO merging to one additional jet, and the associated \(1\sigma \) PDF variation band. The uncertainty is derived from PDF replicas, following the procedure of Ref. [77]. For the range of \(p_{T,\mathrm {cut}}^{j_1}\) investigated, we find that the band is stable, symmetric, and spans roughly \(\pm \,6\%\) at low \(p_{T,\mathrm {cut}}^{j_1}\) to \(\pm \,10\%\) at high \(p_{T,\mathrm {cut}}^{j_1}\). At 27 TeV, we largely investigate a comparable range of Bjorken-x as we do at 13 TeV, and hence expect the PDF uncertainty to remain the same due to the stability of DGLAP evolution. In Fig. 12b, we show the uncertainty associated with the MEPS@NLO jet merging scale, \(Q_{\mathrm {cut}}\). Over a considerable range of \(p_{T,\mathrm {cut}}^{j_1}\), we observe variations below \(\pm \,5\%\) for \(p_{T,\mathrm {cut}}^{j_1}<400\mathrm{~GeV}\) and below \(10\%\) for \(p_{T,\mathrm {cut}}^{j_1}<1\mathrm{~TeV}\), and is comparable to our MC uncertainty.

Fig. 13
figure 13

a Upper: The \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section at \(\sqrt{s}=14\mathrm{~TeV}\) as a function of jet \(p_{T}\) selection criterion \((p_{T,\mathrm {cut}}^{j_1})\), for representative \((\mathcal {Q},m_{\mathcal {Q}})\) combinations. Lower: The same but normalized to the \((\mathcal {Q},m_{\mathcal {Q}})=(T_p,600\mathrm{~GeV})\) curve. b The \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section at \(\sqrt{s} =14\) and 27 TeV, with \(p_{T,\mathrm {cut}}^{j_1}\) = 500 GeV. The error bar reflects the renormalization and factorization scale variation

5 Outlook

In Sect. 3.3, we discussed the dependence of the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section on the mass, spin, and color representation of \(\mathcal {Q}\). A single signal cross section measurement does not, of course, constrain the nature of \(\mathcal {Q}\) uniquely because different spin and color hypotheses can lead to identical cross sections if the mass is chosen/tuned accordingly. For example: in the representative case above, a stop of mass 400 GeV, a fermionic top partner of mass 600 GeV, and a gluino of mass 800 GeV have practically the same \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section for \(p_{T,\mathrm {cut}}^{j_1}\approx 600\) GeV at \(\sqrt{s}=14\) TeV. The degeneracy can be resolved, however, through additional cross section measurements with larger \(p_{T,\mathrm {cut}}^{j_1}\) and/or at a higher center-of-mass energy. Due to the different masses of the representative cases, the cross sections scale differently with increasing jet transverse momentum and the energy. This leads to a break in the degeneracy, thereby enabling one to discriminate against the various \(\mathcal {Q}\) hypotheses.

In the left panel of Fig. 13, we show that all three cases could be distinguished through a secondary measurement of the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section at \( p_{T,\mathrm {cut}}^{j_1}=1 \) TeV were our predictions were precise at the 10% level. The nominal difference between the gluino and fermionic top partner case is roughly that magnitude, while for the case of scalar top, the difference could reach 50%. In light of the uncertainties estimated in Sect. 4, it is unrealistic to claim that one can distinguish a gluino from top partner in this manner, with renormalization and factorization scale uncertainties alone reaching about 20%. Additionally one can also make use of the higher energy collisions, viz. the potential 27 TeV upgrade of LHC, to distinguish these three particles. As shown in Fig. 13b, the degeneracy in cross sections observed for \(p_{T,\mathrm {cut}}^{j_1}\) = 500 GeV at \(\sqrt{s} \) = 14 TeV is lifted at \(\sqrt{s} \) = 27 TeV. Hence, additional information for discriminating against competing \(\mathcal {Q}\) hypotheses is made available by combining data at 14 TeV and 27 TeV.

Alternatively, instead of considering cross sections themselves, one could also consider ratios of the form

$$\begin{aligned} \frac{\sigma ( p_{T}^{j_1}> p_{T,\mathrm {cut}}^{j_1})}{\sigma ( p_{T}^{j_1}> p_{T,\mathrm {cut}}^{j_1}{}^\prime )}\,. \end{aligned}$$
(15)

As discussed in Sect. 3.3, ratios of cross sections measured at \(\sqrt{s}=14\,\)TeV and \(\sqrt{s}=27\,\)TeV can also effectively discriminate against various \(\mathcal {Q}\) hypotheses if theoretical uncertainties on these ratios are smaller than about 30%.

In ratios and double ratios, such as those discussed in previous sections, the normalization component of uncertainties cancel. As shown in Fig. 9a, for example, variations of the factorization and renormalization scales affect mainly the overall normalization, not the shape/\(p_{T}^j\) dependence. Therefore, uncertainties estimated in this way are expected to drop out in Eq. (15). To verify this, we calculate the scale uncertainty for the ratio of the cross section to the nominal scale choice:

$$\begin{aligned}&{{\mathcal {F}}}(p_{T,\mathrm {cut}}^{j_1}, p^{\star }, \mu _R, \mu _F) = \frac{\sigma (p_{T}^{j_1}> p_{T,\mathrm {cut}}^{j_1}, Q_R=\mu _R Q_0, Q_F=\mu _F Q_0)/\sigma (p_{T}^{j_1}> p^{\star }, Q_R=\mu _R Q_0, Q_F=\mu _F Q_0)}{ \sigma (p_{T}^{j_1}> p_{T,\mathrm {cut}}^{j_1}, Q_R=Q_0, Q_F= Q_0)/ \sigma ( p_{T}^{j_1}> p^{\star }, Q_R=Q_0, Q_F= Q_0)}. \end{aligned}$$
(16)

For the normalizing cross section, we consider the \(m_{T_p}=600\) GeV and vary the factorization and renormalization scales in the same manner and choose \(p_{\star }=300\) GeV. We observe that the effect of varying \(\mu _{R,F}\) between 1/2 and 2 with \(Q_0=H_T/2\) is marginal (only 2%) at 1 TeV. The PDF uncertainty, however, impacts the cross section more, though the effect remains no more than 10% at 1 TeV.

Fig. 14
figure 14

The effect of choosing an alternative scale choice \(Q_0 = m_{T_p}\) for \(T_p\bar{T}_p+j\) process using MG5_aMC@NLO +PY8 (left) and Sherpa (right). The envelopes denote scale uncertainty bands. The distributions are normalized to nominal scale \(Q_0 = H_T/2\)

In Eq. (16), we assume the scale choices are correlated for both the numerators and denominators. Fairly, one might argue that in order to estimate uncertainties on cross section ratios through the shape of the \(p_{T,\mathrm {cut}}^{j_1}\) dependence, one should consider different functional forms of \(Q_R\) and \(Q_F\) as well. In Fig. 14 we show that the effect of choosing a non-dynamical nominal scale can in fact induce significant slope variations. Unsurprisingly, deviations are observed for both MG5_aMC@NLO +PY8  and Sherpa events. The ratio of the cross section at 1 TeV and 300 GeV for \(Q_0=m_{T_p}\) is about 15% higher than \(Q_0=H_T/2\). The alternative scale choice \(Q_R = Q_F = m_{T_p}\) is, however, also quite extreme in the sense that it is completely static [98]. For example: the upper envelope for the scale variation (i.e., \(\mu _{R} = \mu _F = 0.5\)) is calculated with \(Q_R =Q_F = 300\) GeV, but the event with \(p_{T,\mathrm {cut}}^{j_1}=1\) TeV implies \(\sqrt{s}> 2\) TeV. Therefore, the ratio \(\sqrt{s}/Q_R\) is a factor of 7, which leads to large non-asymmetric effects as seen in the right panel of Fig. 14.

We conclude that cross section ratios, in general, will provide a more robust way for discrimination but also that the assessment of theoretical uncertainties on such ratios is a highly nontrivial issue. It is natural to expect a majority of the scale uncertainty will cancel in the cross section ratio. However, in general, an NNLO calculation can induce correction terms that scale with jet \(p_T\), and result in an effect may not be small. We have seen in Fig. 9a that different reference scale choices affect the cross section ratio at the level of 10%. The assessment of theoretical uncertainty on such ratio can only be estimated reliably by incorporating higher order calculations. A more quantitative study of this problem is, however, beyond the scope of this paper.

As a final comment, we note that more precise predictions for cross sections can in any case help to improve the discriminating power of the measurements proposed here. Presently, NNLO in QCD predictions for the relevant processes (which include an additional jet in the final state) are not available. For the inclusive process \(pp\rightarrow T_p\bar{T}_p\), however, NNLO predictions can be obtained from HATHOR [34] and top++ [22]. In Fig. 15 we show that the uncertainties on the total cross section are approximately halved when going from NLO to NNLO in \(\alpha _s\). Note that the uncertainty for the \(pp\rightarrow T_p\bar{T}_p+2j\) LO matched cross-section is around 25–30%, while the uncertainty reduces to 10% at NLO. Moreover, the uncertainty in NLO cross-section estimated by three independent sources also show excellent agreement. (Note that for Sherpa, we use the scale scheme \(Q_F = Q_R = m(T_p\bar{T}_p)/2\) to approximate the MG5_aMC@NLO scale scheme of \(Q_F = Q_R = H_T /2\).) The NNLO uncertainty on the total cross section is only around 5%.

Fig. 15
figure 15

Cross section for the \(pp\rightarrow T_p\bar{T}_p\) process with additional light QCD partons at the ME level at leading, next-to-leading and next-to-next-to-leading orders in \(\alpha _s\) provided by the available MC tools

Uncertainties on monojet signal cross section are expected to be larger than the uncertainty on the inclusive \(pp\rightarrow T_p\bar{T}_p\) process. In Fig. 7, we show that renormalization and factorization scale uncertainties are around 40% at LO and 20% at NLO, which is twice as high as for the inclusive \(pp\rightarrow T_p\bar{T}_p\) cross section. This is to be expected because the lowest order term of the total cross section is proportional to \(\alpha _s^2\) in the inclusive case, while the monojet signal cross section is proportional to \(\alpha _s^3\). Assuming that a comparable relative gain can be achieved for \(pp\rightarrow T_p\bar{T}_p+j\) when going to NNLO, it might be possible to get to a theoretical precision of about 10% when NNLO calculations become available. Considering the large momentum transfers involved, EW corrections might have to be added as well in order to capture the effect of (real and virtual) EW Sudakov logarithms [84, 99,100,101].

6 Summary

We have studied the monojet collider signature arising from the process \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\), where \(\mathcal {Q}\) is a heavy colored particle that decays into an invisible particle \(X\) with mass \(m_X\) close to the mass of \(\mathcal {Q}~(m_{\mathcal {Q}})\). In the limit where \((m_{\mathcal {Q}}-m_X)/m_{\mathcal {Q}}\) is small, the full process signature is monojet-like, and the discovery prospects for \(\mathcal {Q}\) and \(X\) become much more challenging, even at the HL-LHC or proposed HE-LHC, due to large SM backgrounds.

In this context, we have investigated the feasibility of extracting the properties of \(\mathcal {Q}\) were a monojet signature discovered, focusing on the observation that the monojet signal cross section \(\sigma (pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j)\) is sensitive to the mass, spin, color representation of the particle. Due to an interplay of these quantities, one cannot readily and uniquely determine the nature of the particle \(\mathcal {Q}\) from a single cross section measurement alone without assuming the mass of the particle.

We have studied several processes calculated at next-to-leading order in QCD with parton shower matching and multijet merging using the state-of-the-art Monte Carlo suites MG5_aMC@NLO +PY8 and Sherpa. It has been pointed out that the dependence on \(p_{T}^j\) in the \(pp\rightarrow \mathcal {Q}\overline{\mathcal {Q}}+j\) cross section, up to overall normalization, depends considerably on the mass of \(\mathcal {Q}\) but exhibits a much milder dependence on its color representation. Moreover, we find that there is little-to-no spin dependence in the shape of the cross section curves. The dependencies of the overall normalization and \(p_{T}^j\) on the underlying parameters can, in principle, allow efficient discrimination of the nature of \(\mathcal {Q}\), and hence the underlying physics behind the monojet signature.

Were the process discovered, there are many obstacles to extracting the cross section of the monojet signature and its \(p_{T}\) dependence. First, the signal is overwhelmed by the background. Therefore, large statistics and a precise understanding of the background is required to observe the signature at the HL-LHC for the parameter region that has not been excluded by current LHC data. Increasing the collider energy improves sensitivity significantly, and represents a promising possibility. In addition, theoretical predictions of the overall cross section normalization need to be sufficiently accurate in order to differentiate against various candidate scenarios. For the relevant model parameter regions investigated, we observe that the overall cross sections need to be predicted within an uncertainty of \(\varDelta \sigma /\sigma <30\)%, while \(p_{T}^j\) dependencies should be known within \(\delta \sigma / \delta p_{T}^j< 10\%/100\mathrm{~GeV}\). Finally, we find that the achievable precision for such quantities is comparable to the required precision.