1 Introduction

The Standard Model of particle physics is a successful theory of nature, although it exhibits many conceptual issues, like the hierarchy problem or the strong CP problem, and practical limitations like the absence of a viable candidate for explaining the dark matter pervading our universe. As a result, it is commonly acknowledged as an effective theory stemming from a more fundamental setup that has still to be observed and confirmed experimentally. This effective description has been recently strengthened with the discovery of a Higgs boson with properties very similar to those expected from the Standard Model in 2012. The null results of all collider searches for new particles predicted by most of beyond the Standard Model theories are, however, at the same time pushing the limits on the mass of these potential new particles to higher and higher energy scales.

Among the viable options for new physics, many extensions of the Standard Model predict the existence of additional quark species that should be observable during the next runs of the Large Hadron Collider (LHC) at CERN. One of the common feature of such theories consists of the predicted vector-like nature of the additional quarks, i.e. their left-handed and right-handed components lie in the same representation of the electroweak symmetry group. These quarks appear for instance in models with extra space-time dimensions that exhibit an extended gauge symmetry or a new strong dynamics giving rise to massive composite states [1,2,3,4,5]. As a result, vector-like quark searches play an important role in the ATLAS and CMS experimental program.

Current searches, relying on signatures induced by both the vector-like quark pair-production and single-production modes, impose strong constraints on the masses of the heavy quarks that are now bounded to be above about 750–1500 GeV [6,7,8,9,10,11,12], this wide range reflecting the wealth of options for describing how the new states decay into a pair of Standard Model particles. Care must, however, be taken with the interpretation of these limits as they are extracted once various simplifying assumptions are accounted for the new physics signal. Most of the bounds indeed assume that the quark partners decay, with a branching ratio of 100%, into third-generation top or bottom quarks, while more general situations where decays into second or first generation quarks are possible are less explored yet. Sizeable couplings to light quarks are still allowed by indirect constraints [13,14,15], which could have a severe impact on electroweak vector-like quark processes at the LHC [16]. The latter, which induce vector-like quark decays, are driven by the couplings of the extra quarks to the weak and Higgs bosons. They admit a simple model-independent parameterisation [17] regardless of the representation of the new physics particles under the electroweak symmetry group, and this bottom-up approach is now adopted by the experimental collaborations. Most LHC searches for heavy quarks hence turn out to be agnostic of the ultraviolet completion of the model and can rather easily be reinterpreted in any framework, and options for combining different searches can also be considered.

Current vector-like quark searches, as well as most associated theoretical work, are, however, based on Monte Carlo simulations of the new physics signals where hard-scattering matrix-elements are evaluated at the leading-order (LO) accuracy in QCD. In addition, event samples featuring different final-state jet multiplicities are sometimes also merged in order to get a better control on the shapes of the key differential distributions. The formal precision of these calculations is nonetheless rather limited, which directly impacts the extraction of any limit on the properties of the new particle properties or the corresponding measurements in the case of a discovery. In this work, we build up a new procedure that allows us to make use of the MadGraph5_aMC@NLO framework [18] to compute total rates and differential distributions at the next-to-leading-order (NLO) accuracy in QCD for processes involving vector-like quarks. This more precisely concerns vector-like quark production on the one hand, and the production of Standard Model particles when vector-like quark diagrams contribute, regardless of the strong or electroweak nature of the Born process.

Our methodology relies on the joint use of the FeynRules [19] and NLOCT [20] packages, the latter making use of FeynArts [21], for automatically generating a UFO library [22] that contains all tree-level vertices and counterterms necessary for NLO QCD computations. This UFO library can then be further used by MadGraph5_aMC@NLO for event generation, at the LO and NLO accuracies in QCD as well as for loop-induced processes. Virtual loop contributions are numerically evaluated through the MadLoop module [23] and combined with the real-emission diagrams following the FKS subtraction method as implemented in MadFKS [24, 25], and the matching to parton showers is finally achieved according to the MC@NLO prescription [26].

In Sect. 2, we detail how we have modified the model-independent parameterisation of Ref. [17] to make it suitable for NLO calculations in QCD matched to parton showers for vector-like quark processes. We also define a series of benchmark scenarios for the phenomenological studies performed in Sect. 3. We investigate vector-like quark pair production and single production (in association with either a jet or a weak gauge or Higgs boson), as well as diboson production. We summarise our work in Sect. 4.

2 A model-independent parameterisation for vector-like quark models

2.1 Model description

Vector-like quarks appear in many extensions of the Standard Model. They are usually included as fields lying in the fundamental representation of \(SU(3)_c\) and carry colour charges similar to those of their Standard Model counterparts. However, they can lie in various representations of the weak interaction symmetry group \(SU(2)_L\) and be assigned different hypercharge \(U(1)_Y\) quantum numbers. Focusing on phenomenologically viable minimal models that comprise a single Standard Model Higgs field \(\varPhi \), only weak triplets, doublets and singlets of vector-like quarks are allowed [27]. Consequently, the particle content of the theory can solely include four species of extra quarks, which we denote by X, T, B and Y, their respective electric charges being \(Q=5/3\), 2 / 3, \(-1/3\) and \(-4/3\).

Although vector-like and Standard Model quarks with the same electric charge mix, the mixing pattern and the resulting phenomenology can be simplified when minimality requirements are imposed. Since we consider that the Higgs sector contains a single Standard-Model-like scalar Higgs field \(\varPhi \), quark mixings are solely generated by its Yukawa interactions. The mass splitting between the vector-like quarks is consequently also constrained to be small and connected to the Higgs vacuum expectation value v, so that the extra quarks will always directly decay into a gauge or Higgs boson and one of the Standard Model quarks [27]. A model-independent effective parameterisation apt to describe the phenomenology of this vector-like quark setup has been recently proposed [17], but it is not suitable for higher-order QCD calculations. The reason is that in the latter parameterisation, the strengths of the interactions of the vector-like and Standard Model quarks with a single Higgs boson depend on the masses of the model particles. As a result, the renormalisation of the quark masses and the one of the couplings are related, which prevents all ultraviolet divergences that arise at the NLO from cancelling. We therefore modify the modelling of Ref. [17] so that all the couplings of the vector-like quarks to a gauge or a Higgs boson are free parameters, and we have the following effective Lagrangian:

(1)

being supplemented to the Standard Model Lagrangian \(\mathcal{L}_\mathrm{SM}\). The terms in the first and second lines consist of gauge-invariant kinetic and mass terms for the vector-like quark fields (taken in the mass eigenbasis) after restricting the covariant derivatives to their QCD component,

$$\begin{aligned} D_\mu = \partial _\mu -i g_s T_a G_\mu ^a. \end{aligned}$$
(2)

The coupling parameter \(g_s\) denotes the strong coupling constant, \(G_\mu \) the gluon field and T (and f for further references) the fundamental representation matrices (the structure constants) of SU(3). Although the electroweak pieces of the covariant derivative could have been included, they have been omitted in order to simplify our model description since they are model-dependent and are expected to yield a negligible effect with respect to their strong interaction part.

The third and fourth lines of Eq. (1) collects the effective interactions of the physical Higgs boson h with one Standard Model quark and its vector-like partner, generation indices being understood for clarity. As mentioned above, such interactions are yielded by the (flavour-changing) Yukawa couplings of the Higgs doublet \(\varPhi \) with the up-type (\(q_u\)), down-type (\(q_d\)) and vector-like quarks that induce a mixing of the Standard Model and the new physics quark sectors. The relevant elements of the mixing matrices have been included in the strengths of the effective interactions \(\hat{\kappa }\).

In the last six lines of Eq. (1), we include the weak interactions of the Z-boson and W-boson with one Standard Model quark and one vector-like quark. In our conventions, we have factorised out the weak coupling g which represents the overall interaction strength, and the \(\kappa \) and \(\tilde{\kappa }\) parameters include the relevant elements of the quark mixing matrices, as for the Higgs interactions. Moreover, the \(c_{\scriptscriptstyle W}\) parameter stands for the cosine of the electroweak mixing angle.

The \(\mathcal{L}_\mathrm{VLQ}\) Lagrangian above is equivalent, at the tree level, to the one of Ref. [17] once we impose

$$\begin{aligned} \begin{aligned} (\hat{\kappa }^{\scriptscriptstyle Q}_{\scriptscriptstyle L,R})_f&=\frac{\kappa _{\scriptscriptstyle Q}m_{\scriptscriptstyle Q}}{v} \sqrt{\frac{\zeta _{\scriptscriptstyle L,R}^f\ \xi ^{\scriptscriptstyle Q}_{\scriptscriptstyle H}}{\varGamma ^{\scriptscriptstyle Q}_{\scriptscriptstyle H}}}, \\ (\tilde{\kappa }^{\scriptscriptstyle Q}_{\scriptscriptstyle L,R})_f&=\kappa _{\scriptscriptstyle Q} \sqrt{\frac{\zeta _{\scriptscriptstyle L,R}^f\ \xi ^{\scriptscriptstyle Q}_{\scriptscriptstyle Z}}{\varGamma ^{\scriptscriptstyle Q}_{\scriptscriptstyle Z}}},\\ (\kappa ^{\scriptscriptstyle Q}_{\scriptscriptstyle L,R})_f&=\kappa _{\scriptscriptstyle Q} \sqrt{\frac{\zeta _{\scriptscriptstyle L,R}^f\ \xi ^{\scriptscriptstyle Q}_{\scriptscriptstyle W}}{\varGamma ^{\scriptscriptstyle Q}_{\scriptscriptstyle W}}}. \end{aligned} \end{aligned}$$
(3)

In this notation, f stands for a generation index and \(\varGamma ^{\scriptscriptstyle Q}_{\scriptscriptstyle X}\) denotes the kinematic factor of the partial decay width of the extra quark Q into a final state containing an X boson. The \(\kappa _{\scriptscriptstyle Q}\) parameter encodes the magnitude of the coupling of the extra quark Q to the different electroweak bosons, while the \(\zeta \) parameters refer to the mixing of the vector-like quarks with the Standard Model quarks,

$$\begin{aligned} \zeta _{\scriptscriptstyle L,R}^f = \frac{|(V_{\scriptscriptstyle L,R})_{4f}|^2}{\sum _{i=1}^3 |(V_{\scriptscriptstyle L,R})_{4j}|^2} \qquad \text {with}\quad \sum _{i=1}^3 \zeta _{\scriptscriptstyle L,R}^i = 1. \end{aligned}$$
(4)

In this expression, we have represented the left-handed and right-handed \(4 \times 4\) mixing matrices between the new quarks and the three Standard Model quarks by \(V_{\scriptscriptstyle L,R}\). Finally, the \(\xi \) parameters of Eq. (3) determine the relative importance of the various decay modes of the vector-like quarks, their sum being equal to one.

The calculation of differential and total cross sections for LHC processes at the NLO accuracy in QCD necessitates to evaluate, on the one hand, real-emission squared amplitudes and on the other hand, interferences of tree-level with virtual one-loop diagrams. The ultraviolet divergences that arise in the latter case are absorbed through the renormalisation of the fields and parameters appearing in \(\mathcal{L}_\mathrm{VLQ}\). This is achieved by replacing all fermionic and non-fermionic bare fields \(\varPsi \) and \(\varPhi \) and bare parameters y by the corresponding renormalised quantities,

$$\begin{aligned} \begin{aligned}&\varPhi \rightarrow \left[ 1 + \frac{1}{2} \delta Z_\varPhi \right] \varPhi ,\\&\varPsi \rightarrow \left[ 1+\frac{1}{2} \delta Z^L_\varPsi P_L + \frac{1}{2} \delta Z^R_\varPsi P_R \right] \varPsi , \\&y \rightarrow y + \delta y, \end{aligned} \end{aligned}$$
(5)

where we truncate the renormalisation constants \(\delta Z\) and \(\delta y\) at the first order in the strong coupling \(\alpha _s=g_s^2/(4 \pi )\). While the wave-function renormalisation constants of the Standard Model quarks are not modified by the presence of the vector-like quarks, the one of the gluon field is given, when the on-shell renormalisation scheme is adopted and when we include \(n_f=5\) massless flavours of quarks, by

$$\begin{aligned} \delta Z_g= & {} -\frac{g_s^2}{24 \pi ^2} \sum _{q=t,T,B,X,Y}\left\{ - \frac{1}{3} + B_0(0,m_q^2,m_q^2) \right. \nonumber \\&\left. \phantom {\frac{g_s^2}{24 \pi ^2}}+\, 2 m_q^2 B_0'(0,m_q^2,m_q^2)\right\} , \end{aligned}$$
(6)

where the \(B_{0,1}\) and \(B_{0,1}'\) functions stand for standard two-point Passarino–Veltman integrals and their derivatives [28]. The left-handed and right-handed wave function renormalisation constants \(\delta Z_{\scriptscriptstyle Q}^{L,R}\) and the mass renormalisation constants \(\delta m_{\scriptscriptstyle Q}\) of a vector-like quark Q (with \(Q=T\), B, X, Y) are similar to the top-quark ones and read

$$\begin{aligned} \delta Z_{\scriptscriptstyle Q}^{L,R}= & {} \frac{g_s^2 C_F}{16 \pi ^2}[ 1 + 2 B_1(m_{\scriptscriptstyle Q}^2;m_{\scriptscriptstyle Q}^2,0)\nonumber \\&+\,8 m_{\scriptscriptstyle Q}^2 B_0'(m_{\scriptscriptstyle Q}^2; m_{\scriptscriptstyle Q}^2,0) + 4 m_{\scriptscriptstyle Q}^2 B_1'(m_{\scriptscriptstyle Q}^2; m_{\scriptscriptstyle Q}^2,0)], \nonumber \\ \delta m_{\scriptscriptstyle Q}= & {} \frac{g_s^2 C_F\ m_{\scriptscriptstyle Q}}{16 \pi ^2}[ 1 - 4 B_0(m_{\scriptscriptstyle Q}^2; m_{\scriptscriptstyle Q}^2,0) \nonumber \\&-2 B_1(m_{\scriptscriptstyle Q}^2; m_{\scriptscriptstyle Q}^2,0)]. \end{aligned}$$
(7)

In these two expressions, \(C_F=(n_c^2-1)/(2 n_c)\) is the quadratic Casimir invariant associated with the fundamental representation of SU(3), with \(n_c=3\). Finally, in order to fix the renormalisation group running of \(\alpha _s\) so that it originates from gluons and the \(n_f=5\) active light quark flavours, we renormalise \(\alpha _s\) by subtracting, from the gluon self-energy, the contributions of all massive particles evaluated at zero-momentum transfer,

$$\begin{aligned} \frac{\delta \alpha _s}{\alpha _s} = \frac{\alpha _s}{2\pi \bar{\epsilon }} \left[ \frac{n_f}{3} \!-\! \frac{11 n_c}{6}\right] \!+\! \frac{\alpha _s}{6\pi }\sum _{q=t,T,B,X,Y} \left[ \frac{1}{\bar{\epsilon }} \!-\! \log \frac{m_q^2}{\mu _R^2}\right] ,\nonumber \\ \end{aligned}$$
(8)

where the ultraviolet-divergent pieces are written in terms of \(\frac{1}{\bar{\epsilon }}=\frac{1}{\epsilon } - \gamma _E + \log {4\pi }\) with \(\gamma _E\) being the Euler–Mascheroni constant and \(\epsilon \) being connected to the number of space-time dimensions \(D=4-2\epsilon \). The first term in the right-hand side of Eq. (8) results from the Standard Model massless parton contributions, while the second term is connected to the massive states, namely the top quark and the four considered vector-like quark species.

In Sect. 3, we will compute predictions at the NLO accuracy in QCD for processes involving vector-like quarks. We will rely on a numerical evaluation of the loop integrals in four dimensions, which necessitates the calculation of rational terms associated with the \(\epsilon \)-dimensional pieces of the loop integrals. There exist two sets of such rational terms that are, respectively, connected to the loop-integral denominators (\(R_1\)) and numerators (\(R_2\)). While the former are universal, the latter are model-dependent and can be seen as a finite number of counterterm Feynman rules derived from the bare Lagrangian [29]. Starting from the \(\mathcal{L}_\mathrm{VLQ}\) Lagrangian of Eq. (1), several \(R_2\) counterterms with external gauge bosons are modified with respect to the Standard Model case,

$$\begin{aligned} R_2^{G G}= & {} \frac{i g_s^2}{96 \pi ^2}\left[ (9 p_2^2 \eta ^{\mu _1\mu _2} - 6 p_2^{\mu _1} p_2^{\mu _2})\phantom {\sum _q}\right. \nonumber \\&\left. +\, 2 \eta ^{\mu _1\mu _2}\sum _q\{p_2^2-6m_q^2\}\right] \delta _{c_1c_2}, \nonumber \\ R_2^{G G G}= & {} -\frac{g_s^3 f_{c_1c_2c_3}}{192 \pi ^2}\left( 33+\sum _q8 \right) V^{\mu _1\mu _2\mu _3},\nonumber \\ R_2^{G G G G}= & {} \frac{i g_s^4}{48 \pi ^2}\Big ( C^{(1)}_{c_1c_2c_3c_4} \eta ^{\mu _1\mu _2} \eta ^{\mu _3\mu _4} \nonumber \\&+\, C^{(2)}_{c_1c_2c_3c_4} \eta ^{\mu _1\mu _3} \eta ^{\mu _2\mu _4}+ C^{(3)}_{c_1c_2c_3c_4} \eta ^{\mu _1\mu _4} \eta ^{\mu _2\mu _3} \Big ),\nonumber \\ R_2^{WWGG}= & {} \frac{i g^2 g_s^2}{96 \pi ^2} V^{\mu _1\mu _2\mu _3\mu _4}\nonumber \\&\times \left( 3+\sum _{Q,f} \Big \{ (\kappa ^{\scriptscriptstyle Q}_{\scriptscriptstyle L})_f^2+ (\kappa ^{\scriptscriptstyle Q}_{\scriptscriptstyle R})_f^2\Big \}\right) \delta _{c_3c_4},\nonumber \\ R_2^{ZZGG}= & {} \frac{i g^2 g_s^2}{288 c_{\scriptscriptstyle W}^2 \pi ^2} V^{\mu _1\mu _2\mu _3\mu _4} \left( 8-18s_{\scriptscriptstyle W}^2 + 20 s_{\scriptscriptstyle W}^4\phantom {\sum _{Q,f}}\right. \nonumber \\&\left. +\,3 \sum _{Q,f} \Big \{ (\tilde{\kappa }^{\scriptscriptstyle Q}_{\scriptscriptstyle L})_f^2+ (\tilde{\kappa }^{\scriptscriptstyle Q}_{\scriptscriptstyle R})_f^2\Big \}\right) \delta _{c_3c_4},\nonumber \\ R_2^{GGhh}= & {} \frac{i g_s^2}{16 \pi ^2}\nonumber \\&\times \eta ^{\mu _1\mu _2} \left( -y_t^2 -2 \sum _{Q,f} \Big \{ (\hat{\kappa }^{\scriptscriptstyle Q}_{\scriptscriptstyle L})_f^2+(\hat{\kappa }^{\scriptscriptstyle Q}_{\scriptscriptstyle R})_f^2\Big \} \right) \delta _{c_3c_4},\nonumber \\ \end{aligned}$$
(9)

with

$$\begin{aligned}&V^{\mu _1\mu _2\mu _3}=(p_1-p_2)^{\mu _3}\eta ^{\mu _1\mu _2} + (p_3-p_1)^{\mu _2}\eta ^{\mu _3\mu _1}\nonumber \\&\quad +\, (p_2-p_3)^{\mu _1}\eta ^{\mu _2\mu _3}, \nonumber \\&V^{\mu _1\mu _2\mu _3\mu _4}=\eta ^{\mu _1\mu _2} \eta ^{\mu _3\mu _4} \!+\! \eta ^{\mu _1\mu _3} \eta ^{\mu _2\mu _4} \!+\! \eta ^{\mu _1\mu _4} \eta ^{\mu _2\mu _3}, \nonumber \\&C_{abcd}^{(1)}=\delta _{ad}\delta _{bc}+\delta _{ac}\delta _{bd}+\delta _{ab}\delta _{cd}\nonumber \\&\quad +\, 2 (\mathrm{Tr}[T_aT_cT_bT_d]+\mathrm{Tr}[T_aT_dT_bT_c])\left( 45+\sum _q11\right) \nonumber \\&\quad -\, 6 (\mathrm{Tr}[T_aT_bT_cT_d]+\mathrm{Tr}[T_aT_bT_dT_c])\left( 7+\sum _q2\right) \nonumber \\&\quad -\,6 (\mathrm{Tr}[T_aT_cT_dT_b]+ \mathrm{Tr}[T_aT_dT_cT_b])\left( 7+\sum _q2\right) , \nonumber \\&C_{abcd}^{(2)} = \delta _{ad}\delta _{bc}+\delta _{ac}\delta _{bd}+\delta _{ab}\delta _{cd} \nonumber \\&\quad -\, 6 (\mathrm{Tr}[T_aT_bT_dT_c]+ \mathrm{Tr}[T_aT_cT_dT_b])\left( 7+\sum _q2\right) \nonumber \\&\quad -\, 2 (\mathrm{Tr}[T_aT_cT_bT_d]+ \mathrm{Tr}[T_aT_dT_bT_c])\left( 24+\sum _q5\right) \nonumber \\&\quad +\,4(\mathrm{Tr}[T_aT_bT_cT_d]+ \mathrm{Tr}[T_aT_dT_cT_b])\left( 24+\sum _q5\right) , \nonumber \\&C_{abcd}^{(3)} = \delta _{ad}\delta _{bc}+\delta _{ac}\delta _{bd}+\delta _{ab}\delta _{cd} \nonumber \\&\quad -\, 6 (\mathrm{Tr}[T_aT_bT_cT_d] + \mathrm{Tr}[T_aT_dT_cT_b])\left( 7+\sum _q2\right) \nonumber \\&\quad -\, 2 (\mathrm{Tr}[T_aT_cT_bT_d]+ \mathrm{Tr}[T_aT_dT_bT_c])\left( 24+\sum _q5\right) \nonumber \\&\quad +\,4(\mathrm{Tr}[T_aT_bT_dT_c]+ \mathrm{Tr}[T_aT_cT_dT_b]) \left( 24+\sum _q5\right) .\quad \end{aligned}$$
(10)

In addition, the new \(R_2\) counterterms involving external vector-like quarks are given by

(11)

In the conventions of Eqs. (9) and (11), \(c_i\), \(\mu _i\), and \(p_i\) indicate the colour index (which can be associated either with the adjoint or the fundamental representation of SU(3)), the Lorentz index, and the four-momentum of the ith particle incoming to the \(R_2^{\ldots i\ldots }\) vertex, respectively. Moreover, an explicit summation upon q, Q and f implies a summation over all quark species, the extra quark species and the Standard Model quark species, respectively.

In the phenomenological study undertaken in Sect. 3, the virtual one-loop contributions to the NLO predictions are evaluated with the MadLoop module [23], and then combined with the real contributions by means of the FKS subtraction method [24] as implemented in the MadFKS package [25]. Both MadLoop and MadFKS being part of MadGraph5_aMC@NLO [18], the entire calculation is entirely automated from the knowledge of the bare Lagrangian of Eq. (1) and the specification of the process of interest [30]. Technically, the translation of the model Lagrangian into a UFO library [22] that contains ultraviolet and \(R_2\) counterterms and that could be used by MadGraph5_aMC@NLO is automatically performed with the FeynRules [19] and NLOCT [20] packages, the latter program taking care of the calculation of the one-loop ingredients of the model files. The corresponding FeynRules and UFO models have been made publicly available and can be downloaded from the webpage http://feynrules.irmp.ucl.ac.be/wiki/NLOModels.

2.2 Benchmark scenarios

Throughout our phenomenological analysis, we adopt several series of benchmark scenarios in which one single vector-like quark is light enough so that it could be reachable at the LHC. Moreover, for the sake of simplicity, we enforce its decay to proceed via a single channel. We denote each class of scenarios by the acronym QVi where the symbol Q can be either T, B, X or Y and refers to the nature of the relevant extra quark, the symbol V refers to the nature of the gauge boson which the vector-like quark Q decays into and i is a generation number related to the family which the quark Q mixes with. For instance, the scenario TW2 would correspond to a setup in which the Standard Model is supplemented by an extra up-type quark T that decays into a final state made of a W-boson and a strange quark with a branching ratio equal to 1.

These types of scenarios are motivated by several considerations. The mixings of the extra quark with the Standard Model sector are severely constrained by flavour-changing neutral current probes [17, 31, 32], LEP data [33, 34] and atomic parity violation measurements [35]. Taking the parameterisation of the \(\kappa \), \(\tilde{\kappa }\) and \(\hat{\kappa }\) parameters of Eq. (3), sizeable mixings with all three generations are only allowed when the \(\kappa _{\scriptscriptstyle Q}\) parameters are below \(10^{-2} \)\( 10^{-3}\) [17]. Those bounds can, however, be relaxed when the mixing pattern is restricted to involve one or two quark generations. In our study, we enforce the vector-like quark mixing to only deal with one specific generation of Standard Model quarks, and we fix the values of the \(\kappa _{\scriptscriptstyle Q}\) parameters to their current experimental limits of 0.07, 0.2 and 0.1 for mixings involving the first, second and third generation, respectively.

Fig. 1
figure 1

Strong (upper) and electroweak (lower) tree-level diagram contributions to the production of a pair of vector-like quarks, possibly carrying the same electric charge. The corresponding loop and real-emission diagrams are automatically generated by MadGraph5_aMC@NLO and contains one extra power of the strong coupling

3 LHC phenomenology

In this section, we compute total cross sections and differential distributions both at the LO and NLO accuracy for several processes involving vector-like quarks. We study the genuine effects of the NLO corrections, as well as the induced reduction of the theoretical uncertainties. We then investigate the effects of matching fixed-order calculations to parton showers, both at LO and NLO. In Sects. 3.1 and 3.2, we, respectively, focus on vector-like quark pair and single production. For these processes, we consider QCD and electroweak diagrams that arise from the Lagrangian of Eq. (1), and we calculate the NLO effects for both contributions and their possible interferences. In Sect. 3.3, we focus on a set of additional processes where vector-like quark effects could be relevant. We consider the production of a single vector-like quark in association with a Standard Model weak or Higgs boson, and of a pair of Standard Model Higgs or weak bosons in the presence of vector-like quarks.

For the considered processes, the central (total and differential) cross-section values are computed after setting the renormalisation and factorisation scales to the average transverse mass of the final-state particles and by using the NLO set of the NNPDF 3.0 parton density functions (PDF) [36] accessed via the LHAPDF 6 library [37]. Scale uncertainties are derived by varying both scales independently by a factor of two up and down, and the PDF uncertainties are extracted following the recommendations of Ref. [38], and both contributions to the theoretical uncertainties are added in quadrature.

3.1 Vector-like quark pair production at the LHC

Vector-like quark pair production is in general dominated by QCD contributions, which has the advantage to be independent of the model details. Model-dependent electroweak diagrams induced by the last four lines of the Lagrangian of Eq. (1) may, however, be non-negligible, in particular when the final state of interest can be produced from the scattering of one or two valence quarks. We focus on the production of a pair of vector-like quarks including the cases where they have the same electric charge,

$$\begin{aligned} p p \rightarrow Q \bar{Q},\quad p p \rightarrow Q Q\quad \text {and}\quad p p \rightarrow \bar{Q} \bar{Q}, \end{aligned}$$
(12)

with Q being either T, B, X or Y. While the first of these three subprocesses receives both strong (diagrams of the first line of Fig. 1) and electroweak (first diagram of the second line of Fig. 1) contributions, the latter two subprocesses can only be mediated by the t-channel exchange of a weak or Higgs boson (last two diagrams of the second line of Fig. 1).

NLO corrections to the strong contributions to the production of a pair of vector-like quarks (the diagrams of the first line of Fig. 1) can be automatically calculated within the MadGraph5_aMC@NLO framework. Thanks to the upgrading of the model to NLO as explained in the previous section, it is now sufficient to type in the program shell, taking the example of \(T\bar{T}\) production,

figure a

With this set of instructions, we first import the UFO model associated with the Lagrangian of Eq. (1) and then start the calculation of the cross section and the generation of Monte Carlo events, at the NLO accuracy in QCD, for the production of a pair of \(T\bar{T}\) quark and antiquark (whose UFO names are +tp+ and +tp +). Other vector-like quark processes can be obtained by replacing the +tp+ symbol by +bp+ (for a B quark), +x+ (for an X quark) and +y+ (for a Y quark). LO event generation can be achieved in the same way once the +[QCD]+ tag is omitted. We recall that the syntax is case insensitive.

For the electroweak channels (the diagrams of the second line of Fig. 1), the command to be typed in reads

figure b

still for the same example of \(T\bar{T}\) production. Other channels (including the production of a pair of heavy quarks carrying the same electric charge) can be addressed similarly. Care must, however, be taken as mixed electroweak and QCD loops appear at NLO. In our approach, we focus on NLO calculations in QCD, and not in QED or in the context of the electroweak theory. As a result, MadGraph5_aMC@NLO automatically discards Feynman diagrams tagged as an electroweak correction to a QCD graph, although in our case all diagrams should be kept for a proper cancellation of all divergences. One possible way to cure this issue would be to include in the UFO model library all ultraviolet and \(R_2\) counterterms that would be necessary for undertaking mixed NLO calculations in QCD and QED and to implement in MadFKS the necessary subtraction terms. This, however, goes beyond the scope of this work, so that in order to maintain automation from the user standpoint, we have instead released a public script that should be called for diagram generation in order to prevent MadGraph5_aMC@NLO from discarding any loop diagram that would be tagged as an electroweak correction to a QCD Born diagram. More precisely, the script allows for the inclusion of all box diagrams containing at least two strong interaction vertices, but it removes weak-boson or Higgs-boson loop contributions that consist of an electroweak correction to a QCD Born process. Additionally, diagrams exhibiting the t-channel exchange of a weak or Higgs boson but with an additional gluon are kept. Not using the script instead yields the removal of several necessary box diagram contributions.

Finally, QCD and electroweak diagrams can interfere. Because of the mixing of QCD and electroweak interaction orders at the one-loop level and the missing counterterms and subtraction terms that have been mentioned above, MadGraph5_aMC@NLO cannot currently be used for the calculation of these interferences beyond the LO accuracy. We therefore rely on LO simulations and reweight instead the results to include a K-factor approximating the effect of the QCD corrections, denoted by \(K_\mathrm{NLOQCD}^{(\mathrm int)}\), taken as the ratio of the NLO to the LO (differential) results. We choose this \(K_\mathrm{NLOQCD}^{(\mathrm int)}\) factor to be the geometrical average of the \(K_\mathrm{NLOQCD}^{(\mathrm QCD)}\) and \(K_\mathrm{NLOQCD}^{(\mathrm EW)}\) K-factors obtained in the context of the QCD and electroweak diagrams taken independently,

$$\begin{aligned} K_\mathrm{NLOQCD}^{(\mathrm int)} = \sqrt{K_\mathrm{NLOQCD}^{(\mathrm QCD)}\ K_\mathrm{NLOQCD}^{(\mathrm EW)}}, \end{aligned}$$
(13)

where our notations explicitly indicate the nature of the corresponding underlying Born process. We have additionally checked that, for the central value, our procedure yields numerical differences that are of at most 2%. The related MadGraph5_aMC@NLO command allowing for generating interference events is, again for the example of \(T\bar{T}\) production,

figure c

which is a standard command for LO event generation in MadGraph5_aMC@NLO.

All MadGraph5_aMC@NLO scripts necessary for differential and total cross-section calculations and event generation are available from the webpage http://feynrules.irmp.ucl.ac.be/wiki/NLOModels, together with the UFO and FeynRules models.

Fig. 2
figure 2

LO and NLO QCD inclusive cross sections for \(T\bar{T}\)/TT/\(\bar{T} \bar{T}\) pair production at the LHC with \(\sqrt{s}=13\) TeV. The QCD contribution results are presented together with the associated theoretical uncertainty bands, and we indicate the shifts in the bands that are induced by including weak or Higgs-boson exchange diagram contributions when they are non-negligible, which is only the case for the scenario TH1

In Fig. 2 and Table 1, we present LO and NLO total cross sections for the production of a pair of vector-like T quarks for the different scenarios introduced in Sect. 2.2, and we depict the dependence of the cross sections on the vector-like quark mass. We first focus on the pure QCD contribution that is independent of the vector-like quark nature (diagrams of the upper line of Fig. 1). The genuine NLO contributions are found to be important as they first induce a shift in the cross section of about 50% within the entire probed \(m_{\scriptscriptstyle Q}\) range and next reduce the dependence of the rate on the unphysical factorisation and renormalisation scales. The uncertainty band is indeed significantly reduced when NLO effects are accounted for, the scale dependence being reduced to the level of about \(\pm 10\%\) over the entire mass range. At the NLO, the scale dependence induced by the virtual contributions indeed partially compensates for the one stemming from the Born and real-emission diagrams. Our results for the pure QCD case (strong production Born diagrams only) agree with the literature, and we recall that this corresponds to the production of a top-quark pair with a different top-quark mass [39]. It is in contrast the first calculation including the impact of the electroweak diagrams at the NLO accuracy in QCD.

Table 1 LO and NLO QCD inclusive cross sections for \(T\bar{T}\)/TT/\(\bar{T}\bar{T}\) production at the LHC, running at a center-of-mass energy of \(\sqrt{s}=13\) TeV. The results are shown together with the associated scale and PDF relative uncertainties. For all scenarios but the TH1 one, the predictions match the pure QCD results

In the considered scenarios, electroweak contributions to the production of a pair of vector-like quarks possibly carrying the same electric charge are found to be important only for the TH1 scenario, the associated results being shown on Fig. 2 as dashed and dotted bands. These bands are, respectively, related to the production of a pair of vector-like quark–antiquark (dashed) and to the production of a pair of heavy quark regardless of their electric charge (dotted). In this last case, the contributions from the three processes of Eq. (12) are summed over. The TH1 scenario features an extra quark that mixes with the first generation of Standard Model quarks so that parton density effects could lead to an enhancement of the production rate due to quark–antiquark, quark–quark and antiquark–antiquark (electroweak) scattering diagrams involving one or two initial valence quarks. This is particularly pronounced for setups featuring heavy vector-like quarks that require one to probe large Bjorken-x phase-space regions. As a consequence, the central cross-section values and the scale and parton density uncertainties are different from the pure QCD context since non-QCD diagrams (featuring a different initial state) dominate. This is illustrated in Table 1 for a few mass choices. Whereas the production rates are always larger, the uncertainties can be either smaller or larger than in the QCD case.

We observe a huge gain in cross section for TH1 scenarios with a very heavy extra quark. This stems from Eq. (3) that shows that the coupling of the extra quark to the Higgs boson and a lighter quark has been taken proportional to the vector-like quark mass, and is thus enhanced for large values of \(m_{\scriptscriptstyle Q}\). In principle, such a coupling should also be proportional to the related mixing matrix element \(\zeta \) that compensates this enhancement, as shown in Eq. (3) and in Eq. (4). Setting \(\zeta =1\), this feature is translated into the adopted value for the \(\kappa _Q\) parameters.

Similar properties can be found in the context of the production of B, X and Y vector-like quarks, as illustrated in Appendix A.

Fig. 3
figure 3

Properties of the pair-produced vector-like B-quarks. We present the transverse-momentum (upper left) and pseudorapidity (upper right) distribution for the first B-quark, as well as the transverse-momentum spectrum of the B-quark pair (lower left for the [0, 1000] GeV transverse-momentum region and lower right for a zoom in the [0, 500] GeV region). We compare fixed-order QCD predictions at the LO accuracy (purple dashed curve), NLO accuracy (blue dashed curve), and after matching these two calculations to parton showers (green and red dashed curves for the LO and NLO cases, respectively), and the solid lines depict the results once the electroweak diagram contributions are included. We have fixed the B-quark mass either to 500 GeV (upper series of curves) or to 1500 GeV (lower series of curves)

Accurate differential distributions are often helpful for setting more precise exclusion limits and refine the experimental search strategies. Our implementation in the MadGraph5_aMC@NLO platform can be used to this aim, and we present in Fig. 3 several differential distributions in several observables, including NLO and parton-shower effects. We have chosen the BH2 class of benchmark scenarios with a vector-like quark mass set either to 500 GeV (upper series of curves on the figure) or 1500 GeV (lower series of curves on the figure). In our calculations, we have made use of the MadSpin [40] and MadWidth [41] programs for automatically taking care of the heavy quark decays in a way in which both off-shell and spin correlation effects are retained, matched the fixed-order calculation with the parton-shower and hadronisation infrastructure as modelled by the Pythia 8 package [42], and we have reconstructed all final-state jets by means of the anti-\(k_T\) algorithm [43] with a radius parameter set to 0.5 as implemented in FastJet [44]. As shown on Figs. 3 and 4, our predictions confirm the total cross-section findings of Table 4 (see Appendix A). The contributions of the electroweak diagrams are, respectively, negligible and significant for light and heavy vector-like quarks. Moreover, the parton density uncertainties dominate for setups exhibiting a large \(M_{\scriptscriptstyle B}\) value, rendering the theoretical predictions barely reliable.

Fig. 4
figure 4

Same as in Fig. 3 but for the distribution in the transverse momentum of the first three leading jets and the \(H_T\) variable defined as the sum of the transverse momentum of all the final-state jets and leptons

In the two upper panels of Fig. 3, we study the properties of the first produced B-quark and show its transverse momentum \(p_T(B_1)\) and pseudorapidity \(\eta (B_1)\) distributions. In the case of a light B-quark, the parton-shower effects (green and red solid lines) slightly affect the shapes of the spectra predicted by the fixed-order calculations (blue and brown solid lines), both at the LO and NLO accuracies. For heavier vector-like quarks, slight modifications can be observed in the small \(p_T\) and large \(|\eta |\) region although the accurate modelling of the first extra jet does not yield any impact at the level of the individual B-quarks. These differences are largely covered by the theoretical uncertainties stemming from the poor knowledge of the parton densities in the relevant phase-space regions. These PDF uncertainties are dominant so that the reduction at the level of the scale uncertainties has a small impact. However, PDF uncertainties are expected to improve in the coming years thanks to new LHC data, so that it is mandatory to have NLO calculations available to get a better control on the predictions. In contrast, parton-shower effects are directly visible when distributions related to the two B-quarks considered as a pair are considered (Fig. 3). Focusing on the related transverse-momentum distributions, the fixed-order predictions (for which only the \(p_T(B_1+B_2)=0\) bin is populated at LO) diverge at small \(p_T\) due to soft and collinear radiation giving rise to large logarithms that must be resummed to all orders to obtain reliable predictions. This resummation is effectively achieved by matching the fixed-order calculations to parton showers, and the resulting distribution exhibits a reliable behaviour with a peak for \(p_T(B_1+B_2)\) of about 10–20 GeV. Uncertainties originating from the choice of the shower algorithm and its inherent free parameters are, however, not estimated.

The magnitude of the electroweak diagrams is also studied (dashed lines). In the case of lighter B-quarks, the theoretical predictions are essentially driven by the QCD contributions so that no differences can be noticed. This contrasts with the heavy B-quark case where electroweak diagrams enhance the total rate by about 30% (see Appendix A) and also impact the differential distributions both in terms of normalisation and shape. This originates from the different topologies of the electroweak diagrams that feature a t-channel colourless boson exchange.

Fig. 5
figure 5

Representative Feynman diagram for single vector-like quark production in association with a quark. Other diagrams exist for final-state antiquarks, and all flavour assignments are understood. Virtual- and real-emission diagrams necessary for QCD correction calculations can be generated by MadGraph5_aMC@NLO

Fig. 6
figure 6

LO and NLO QCD inclusive cross sections for single B quark production at the LHC with \(\sqrt{s}=13\) TeV. The results are presented together with the associated theoretical uncertainty bands for the BZ1 (upper left), BZ2 (upper right), BW1 (lower left) and BW2 (lower right) scenarios

We now include the vector-like quark decays into a Higgs boson and a strange quark and reconstruct the final-state jets as detailed above. Considering only hard central jets for which \(|\eta |<2.4\) and \(p_T>30\) GeV, we present the transverse-momentum distributions of the three leading jets as well as the spectrum of the \(H_T\) variable defined as the scalar sum of the transverse momenta of all final-state jets and leptons in Fig. 4, the generated events being inclusive in the Higgs-boson decays. Leptons are included only if their transverse momentum is larger than 30 GeV, their pseudorapidity smaller than 2.4 in absolute value and if they are isolated from any jet by an angular distance in the transverse plane \(\varDelta R\) of at least 0.5. Focusing on fixed-order predictions matched to parton showers, we observe that the first two leading jets are in general hard, as they result from the decay of two heavy coloured particles. In contrast, the structure of the \(p_T\) dependence of the third jet is more representative of the one expected from a radiation jet, this jet being most of the time arising from initial-state or final-state radiation. Turning to the \(H_T\) distribution (lower right panel of the figure), one notices a peak for very small \(H_T\) values when the B-quark mass is fixed at 500 GeV. This feature arises from events where the two jets originating from the heavy quark decays are mis-reconstructed, the leading jet being thus the radiation jet so that the associated jet activity in the events is not significantly large.

Table 2 LO and NLO QCD inclusive cross sections for single B production at the LHC, running at a center-of-mass energy of \(\sqrt{s}=13\) TeV. The results are shown together with the associated scale and PDF relative uncertainties in the context of several benchmark scenarios

3.2 Single vector-like quark production in association with jets

Single vector-like quark production mechanisms are of a pure electroweak nature. The associated predictions are thus model-dependent as the sizes of the electroweak vector-like quark couplings are free parameters of the model. Comparing vector-like quark single and pair production, the latter gets an enhancement originating from the presence of strong diagram contributions (first line of Fig. 1) together with a phase-space suppression for large vector-like quark mass values. In contrast, electroweak single vector-like quark production is less suppressed for large vector-like quark masses, which could compensate the weakness of the involved interaction vertices and make this channel the main LHC discovery mode for a heavy vector-like quark. As a results, several ATLAS and CMS vector-like quark searches also target their single-production mode [10, 45,46,47,48].

In this section, we focus on vector-like quark single production in association with jets,

$$\begin{aligned} p p \rightarrow Q j\ \text {or}\ \bar{Q} j \quad \text {with}\ Q = T, B, X \ \text {or}\ Y. \end{aligned}$$
(14)

Other single-production mechanisms exist, with, for instance, a final-state gauge or a Higgs boson, but we refer to Sect. 3.3 for the latter. A representative set of Feynman diagrams related to single vector-like quark production with jets is shown on Fig. 5, and NLO cross-section calculation and event generation can be achieved with MadGraph5_aMC@NLO by typing in the program interpreter the command

figure d

for single T production. Other processes with a different final-state vector-like quark can be accounted for with a similar syntax, and LO event generation only necessitates to remove the +[QCD]+ tag. As for electroweak contributions to vector-like quark pair production, mixed QCD and electroweak loop diagrams appear at the NLO level and must be treated consistently for getting ultraviolet-finite results. This step being not automated, we provide scripts that steer the event generation process on the UFO model webpage.

In Fig. 6 and Table 2, we present LO and NLO total cross section for single B quark production for the BZ1, BZ2, BW1 and BW2 scenarios introduced in Sect. 2.2 for which single vector-like B quark production occurs via Z-boson or W-boson exchanges. This consists of the first NLO predictions for a single vector-like quark production process. Although all depicted total cross sections exhibit a similar order of magnitude, a small hierarchy is observed. It is driven by an interplay of the parton densities that enhance mechanisms involving first generation quarks and of the new physics coupling strengths that are much larger for vector-like quark mixings with second generation (\(\kappa _{\scriptscriptstyle Q}=0.2\)) than with first generation quarks (\(\kappa _{\scriptscriptstyle Q}=0.07\)). For vector-like B-quark masses smaller than 500 GeV, single-production cross sections are slightly suppressed by a factor of 2 or 3 with respect to the strong production of a pair of B quarks (see Appendix A), but feature a similar order of magnitude for \(m_{\scriptscriptstyle B} \in [500, 800]\) GeV. In contrast, single-B production always dominates for heavier vector-like quarks by up to two orders of magnitude for larger values of \(m_{\scriptscriptstyle B}\). For scenarios with smaller \(\kappa _{\scriptscriptstyle Q}\) values, the changes in the relative importance of the two production channels can, however, be shifted towards higher masses.

Fig. 7
figure 7

Differential distributions depicting the properties of a singly produced vector-like T (or \(\bar{T}\)) quark. We present its transverse momentum (upper), rapidity (lower left) and pseudorapidity (lower right) spectrum and compare fixed-order predictions at the LO accuracy (purple curve) and NLO accuracy (blue curve), as well as prediction including the matching of these two calculations to parton showers (green and red curves for the LO and NLO cases, respectively). We have fixed the heavy quark mass either to 500 GeV (upper series of curves) or to 1500 GeV (lower series of curves)

In Fig. 7, we turn to the study of differential distributions related to inclusive single T-quark production at the LHC, focusing on the TZ1 scenario as an illustrative benchmark point and for a vector-like quark mass \(m_{\scriptscriptstyle T}\) fixed to 500 GeV or 1500 GeV. Event generation and reconstruction is performed following the guidelines mentioned in Sect. 3.1, and we present the transverse momentum \(p_T(T)\), rapidity y(T) and pseudorapidity \(\eta (T)\) spectrum of the vector-like quark. We show predictions both at the fixed-order (purple and blue curves for the LO and NLO accuracy, respectively) and after matching the results to parton showers (green and red curves for the LO and NLO accuracy, respectively). Theoretical uncertainties originating from scale variations and parton densities are included for the matched predictions. In general, NLO effects only moderately affect the shape of the different spectra but in contrast drastically reduce the theoretical scale uncertainties and allow for a better prediction of the spectrum normalisation. Similarly to the pair-production case, matching to parton showers only mildly impacts fixed-order predictions for the properties of the produced heavy quarks. Regardless of \(m_{\scriptscriptstyle T}\), the transverse-momentum distribution of singly produced T quarks exhibits a typical steeply falling behaviour with increasing values of \(p_T(T)\) (with a peak at low \(p_T\), which is invisible due to the binning choice), and the vector-like quark is quite forward by virtue of the process topology, with \(|\eta |>2\) in average.

Fig. 8
figure 8

Differential distributions depicting the properties of the decay products of a singly produced vector-like T (or \(\bar{T}\)) quark. We present the transverse-momentum distributions of the three leading jets, as well as the one in the \(H_T\) variable defined as the scalar sum of all jet and isolated lepton transverse momenta. We compare LO (green) and NLO (red) predictions after matching the fixed-order calculations to parton showers. We have fixed the heavy quark mass either to 500 GeV or to 1500 GeV

In Fig. 8, we study the properties of the decay products of the heavy quarks that consist of a Z-boson and an up quark in the TZ1 scenario. We first present the transverse-momentum spectrum of the three leading jets with a pseudorapidity satisfying \(|\eta |<2.5\) and a transverse-momentum larger than 30 GeV. Being inclusive in the Z-boson decay, the jets can originate either from the heavy T-quark decay, from the Z-boson decay or from initial- or final-state radiation. Focusing on the leading jet \(p_T\) distribution, we observe that it peaks at half the T-quark mass, which shows that the leading jet is often issued from the heavy quark decay. In contrast, the second jet transverse-momentum distribution exhibit a plateau extending up to half the T-mass, which we conclude that it could alternatively originate either directly from the T-quark decay, or from the Z-boson decay. The third jet finally shows a different behaviour, the distribution peaking this time at a lower \(p_T\) value, so that it is likely to be connected to the hard process.

In addition to the leading jet transverse-momentum distributions, we also present the distribution in the \(H_T\) variable defined as the scalar sum of the transverse momentum of the final-state jets and isolated leptons, the latter being only considered if their pseudorapidity fulfils \(\eta |<2.4\) and their transverse momentum \(p_T>30\) GeV. Moreover, we require the leptons to be well separated from any jet, imposing the angular distance \(\varDelta R\) to be larger than 0.5. The \(H_T\) variable exhibits a non-trivial structure which peaks both at the T-quark mass \(m_{\scriptscriptstyle T}\) and at \(m_{\scriptscriptstyle T}/2\), the second peak arising from cases where the Z-boson decays invisibly.

In all cases, the shapes turn out to be only slightly affected by the NLO effects and the uncertainties are drastically reduced.

Fig. 9
figure 9

Representative Feynman diagram for Higgs-boson pair production that depend on vector-like quark exchange (two leftmost diagrams) and for single vector-like quark production in association with a Higgs boson (two rightmost diagrams). Similar diagrams exist when Higgs bosons are replaced by weak gauge bosons. With the second diagram, we include LO loop-induced contributions to Higgs-boson pair production that can only be calculated at the LO accuracy with MadGraph5_aMC@NLO

Table 3 LO and NLO QCD inclusive cross sections for Higgs-boson (upper) and Z-boson (lower) pair production at the LHC, running at a center-of-mass energy of \(\sqrt{s}=13\) TeV. The results are shown together with the associated scale and PDF relative uncertainties in the context of the TH1 (upper) and TZ1 (lower) class of benchmark scenarios

3.3 Other processes impacted by the presence of vector-like quarks

Vector-like quarks can also affect Standard Model processes due to additional Feynman diagram contributions featuring a virtual heavy quark. For instance, the production of a pair of Higgs bosons could proceed via the new physics diagrams shown in Fig. 9, similar diagrams existing for the diboson case. MadGraph5_aMC@NLO can be used for event generation at the NLO accuracy, once topologies featuring intermediate heavy quark resonances are treated accordingly. This is achieved with the command

figure e

for di-Higgs and diboson production, respectively, using the $$ symbol to remove any possible intermediate resonance. In the last case, the pair of symbols v v stands either for w+ w- or z z according to the weak boson under consideration. LO event generation can finally be achieved by removing the +[QCD]+ tag.

Additional single vector-like quark production processes where the heavy quark is produced in association with a Standard Model boson can be considered as extra mechanisms useful for seeking for vector-like quarks (diagrams shown in the rightmost part of Fig. 9). Such processes can be simulated with MadGraph5_aMC@NLO, by typing in the commands

figure f

for \(T/\bar{T}\) production, as an example. Once again, the $$ symbol is used in order to avoid intermediate resonances. VQ production can be undertaken similarly.

Fig. 10
figure 10

LO and NLO QCD inclusive cross sections for Higgs-boson (left) and Z-boson (right) pair production at the LHC with \(\sqrt{s}=13\) TeV. The results are presented together with the associated theoretical uncertainty bands for the TH1 (left) and TZ1 (right) scenarios

As an illustrative example, we show in Table 3 and Fig. 10 total rates at LO and NLO for Higgs-boson (upper) and Z-boson (lower) pair production, respectively, for scenarios featuring a T quark interacting with the first generation of Standard Model quarks and either the Higgs boson or the Z-boson. We observe gigantic K-factors in the case of di-Higgs production. This enhancement is connected to an interplay of two effects. Turning back to the observations made in Sect. 3.1, the vector-like quark coupling to light quarks and a Higgs boson has strength that is proportional to the heavy quark mass, and thus becomes stronger and stronger with increasing vector-like quark masses. In addition, a new channel where the initial state is comprised of a gluon and a quark opens up at NLO. This component of the NLO cross section turns out to dominate due to the large gluon density in the proton. As a result, the total cross section for vector-like-quark-mediated di-Higgs production is more or less constant with the heavy quark mass at the NLO QCD accuracy, which contrasts with the LO case.

4 Conclusions

We have modified a previously introduced model-independent parameterisation suitable for the study of the vector-like quark phenomenology at the LHC so that it is now suitable for NLO calculations in QCD matched to parton showers within the MadGraph5_aMC@NLO framework. We have illustrated its usage in the context of vector-like quark pair production and vector-like quark single production in association either with a jet or with a weak or Higgs boson. For all showcased processes, we have considered QCD and electroweak diagram contributions and investigated NLO and parton-shower effects on the normalisation and shapes of the associated kinematical distributions.

Fig. 11
figure 11

LO and NLO QCD inclusive cross sections for B pair production at the LHC, for \(\sqrt{s}=13\) TeV. The QCD contribution results are presented together with the associated theoretical uncertainty bands, and we indicate the shifts in the bands that are induced by weak diagram contributions in a scenario in which they are non-negligible

Table 4 LO and NLO QCD inclusive cross sections for B pair production at the LHC, running at a center-of-mass energy of \(\sqrt{s}=13\) TeV. The results are shown together with the associated scale and PDF relative uncertainties in the context of several benchmark scenarios. For all non-indicated scenarios, the results match the pure QCD one given in Table 1

We have found that NLO K-factors are important, both globally (at the total-rate level) and at the differential distribution level and could hence potentially impact limits currently extracted from vector-like quark search results of the ATLAS and CMS collaborations. We have in particular noticed the existence of potentially huge K-factors for new physics scenarios involving the coupling of a heavy vector-like quark to first generation quarks and a Higgs boson due to new production channels that open at the NLO accuracy. This motivates further investigations, in particular to assess how an experimental analysis including detector effects could benefit from the gain in cross section stemming from the new topologies that dominate at NLO and to determine the impact on the current vector-like quark limits and LHC discovery potential.