1 Introduction

The standard model (SM) has been outstandingly successful in describing a wide range of fundamental phenomena. However, one of its notable shortcomings is that it does not provide a natural explanation for the Higgs boson (\(\mathrm {H}\)) [1,2,3] observed at 125\(\,\text {Ge}\text {V}\)  [4, 5] having a mass that is comparable to the electroweak scale. The suppression of divergent loop corrections to the Higgs boson mass requires either fine-tuning of the SM parameters or new particles at the \(\,\text {Te}\text {V}\) scale. Many theories of beyond-the-SM physics phenomena that attempt to solve this hierarchy problem predict new particles, which could be partners of the top and bottom quarks and thus cancel the leading loop corrections. Vector-like quarks (VLQs) represent one class of such particles among those that have fermionic properties. Their left- and right-handed components transform in the same way under the SM symmetry group \(\mathrm {SU(3)_C{\times }SU(2)_L{\times }U(1)_Y}\) [6]. This property allows them to have a gauge-invariant mass term in the Lagrangian of the form \(\overline{\psi }\psi \), where \(\psi \) represents the fermion field; hence, their masses are not determined by their Yukawa couplings to the Higgs boson. These quarks are not ruled out by the measured properties of the Higgs boson. They are predicted in many beyond-the-SM scenarios such as grand unified theories [7], beautiful mirrors [8], models with extra dimensions [9], little Higgs [10,11,12], and composite Higgs models [13], as well as theories proposed to explain the SM flavor structure [14] and solve the strong CP problem [15].

The VLQs can be produced singly or in pairs [6]. The cross section for single-quark production is model dependent and depends on the couplings of the VLQs to the SM quarks. On the other hand, pair production of VLQs occurs via the strong interaction, and its cross section is uniquely determined by the mass of the VLQ. Another characteristic of the VLQs is their flavor-changing neutral current decay, which distinguishes them from chiral fermions. The top and bottom quark VLQ partners \(\mathrm {T}\) and \(\mathrm {B}\) are expected to couple to the SM third-generation quarks [16], and decay via \(\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}, \mathrm {t}\mathrm {Z}, \mathrm {t}\mathrm {H} \) and \(\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}, \mathrm {b}\mathrm {Z}, \mathrm {b}\mathrm {H} \), respectively.

In this paper, a search for the production of \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) is presented, where at least one of the \(\mathrm {T}\) (\(\mathrm {B}\)) quarks decays as \(\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}\) (\(\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}\)), as shown in Fig. 1. The search is performed using events with two oppositely charged electrons or muons, consistent with coming from a decay of a \(\mathrm {Z}\) boson, and jets. The data were collected with the CMS detector at the CERN LHC in 2016, from proton–proton (\({\mathrm {p}}{\mathrm {p}}\)) collisions at \(\sqrt{s} = 13\) \(\,\text {Te}\text {V}\), corresponding to an integrated luminosity of 35.9\(\,\text {fb}^{-1}\).

Fig. 1
figure 1

Leading-order Feynman diagrams for the pair production and decay of \(\mathrm {T}\) (left) and \(\mathrm {B}\) (right) VLQs relevant to final states considered in this analysis

Searches for the pair production of \(\mathrm {T}\) and \(\mathrm {B}\) quarks have previously been reported by the ATLAS [17,18,19,20] and CMS [21,22,23] Collaborations. The strictest lower limits on the \(\mathrm {T}\) and \(\mathrm {B}\) quark masses range between 790 and 1350\(\,\text {Ge}\text {V}\), depending on the decay mode studied. The mass range for the \(\mathrm {T}\) and \(\mathrm {B}\) quarks studied in this analysis is 800–1500\(\,\text {Ge}\text {V}\).

2 The CMS detector and event simulation

The central feature of the CMS apparatus is a superconducting solenoid of 6\(\,\text {m}\) internal diameter, providing a magnetic field of 3.8\(\,\text {T}\). Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (\(\eta \)) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [24].

Events of interest are selected using a two-tiered trigger system [25]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100\(\,\text {kHz}\) within a time interval of less than 4\(\,\upmu \text {s}\). The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1\(\,\text {kHz}\) before data storage.

Monte Carlo (MC) simulated signal events of the processes \({\mathrm {p}}{\mathrm {p}}\rightarrow \mathrm {T} \overline{\mathrm {T}} \) and \({\mathrm {p}}{\mathrm {p}}\rightarrow \mathrm {B} \overline{\mathrm {B}} \) for \(\mathrm {T}\) and \(\mathrm {B}\) quark masses in the range 0.8–1.5\(\,\text {Te}\text {V}\) are produced in steps of 0.1\(\,\text {Te}\text {V}\). The events are generated with MadGraph 5_amc@nlo 2.3.3 [26], where the processes are produced at leading order (LO) with up to two partons in the matrix element calculations, using the NNPDF3.0 parton distribution function (PDF) set [27]. Showering and hadronization is simulated with pythia  8.212 [28] using the underlying event tune CUETP8M1 [29]. To normalize the simulated signal samples to the data, next-to-next-to-leading-order (NNLO) and next-to-next-to-leading-logarithmic (NNLL) soft-gluon resummation cross sections are obtained using the Top++ program (v.2.0) [30], with the MSTW2008NNLO68CL PDF set as implemented in the LHAPDF (v.5.9.0) framework [31].

The main background process is Drell–Yan (\(\mathrm {Z}\)/\(\gamma ^{*}\))+jets production, with smaller contributions from \(\mathrm {t}\overline{\mathrm {t}}\)+jets and \(\mathrm {t}\overline{\mathrm {t}}\mathrm {Z} \). Throughout the paper this background will be referred to as DY+jets. Other backgrounds, such as diboson, \(\mathrm {t}\mathrm {Z} \mathrm {q} \), \(\mathrm {t}\mathrm {W}\mathrm {Z} \), and \(\mathrm {t}\overline{\mathrm {t}}\mathrm {W}\) production, are considerably smaller. The DY+jets simulated background samples are generated in different bins of the \(\mathrm {Z}\) boson transverse momentum \(p_{\mathrm {T}}\), using the mc@nlo [32] event generator at NLO precision with the FxFx jet-matching scheme [33]. The \(\mathrm {t}\overline{\mathrm {t}}\)+jets events are generated using the powheg  2.0 [34,35,36] generator. The generated events are interfaced with pythia  8.212 [28] for shower modeling and hadronization, using the underlying event tune CUETP8M2T4 [37] for \(\mathrm {t}\overline{\mathrm {t}}\)+jets simulation and CUETP8M1 [29] for the DY+jets process. The SM diboson events are also produced using the same standalone pythia  8.212 generator. The production of rare single top processes \(\mathrm {t}\mathrm {Z} \mathrm {q} \) and \(\mathrm {t}\mathrm {W}\mathrm {Z} \), as well as a \({\mathrm {t}\overline{\mathrm {t}}}\) pair in association with a \(\mathrm {W}\)or \(\mathrm {Z}\) boson, are simulated with up to one additional parton in the matrix element calculations using the MadGraph 5_amc@nlo 2.3.3 [26] generator at LO precision and matched with the parton showering predictions using the MLM matching scheme [38].

Backgrounds are normalized according to the theoretical predictions for the corresponding cross sections. The DY+jets production cross sections from the mc@nlo  [32] generator are valid up to NLO. Using a top quark mass of 172.5\(\,\text {Ge}\text {V}\), the \(\mathrm {t}\overline{\mathrm {t}}\)+jets production cross section at NNLO [30] is determined. Diboson production is calculated at NLO for \(\mathrm {W}\mathrm {Z} \) [39] and NNLO for \(\mathrm {Z} \mathrm {Z} \) [40] and \(\mathrm {W}\mathrm {W}\) [41]. The production cross sections for the rare processes \(\mathrm {t}\mathrm {Z} \mathrm {q} \), \(\mathrm {t}\mathrm {W}\mathrm {Z} \), and \(\mathrm {t}\overline{\mathrm {t}}\mathrm {W}\) are calculated at NLO [42].

A Geant4-based [43, 44] simulation of the CMS apparatus is used to model the detector response, followed by event reconstruction using the same software configuration as for the collision data. The effect of additional \({\mathrm {p}}{\mathrm {p}}\) interactions in the same or nearby bunch crossings (pileup) in concurrence with the hard scattering interaction is simulated using the pythia  8.1 generator and a total inelastic \({\mathrm {p}}{\mathrm {p}}\) cross section of 69.2\(\,\text {mb}\)  [42]. The frequency distribution of the additional events is adjusted to match that observed in data and has a mean of 23.

3 Event reconstruction

The event reconstruction in CMS uses a particle-flow (PF) algorithm [45] to reconstruct a set of physics objects (charged and neutral hadrons, electrons, muons, and photons) using an optimized combination of information from the subdetectors. The energy calibration is performed separately for each particle type.

The \({\mathrm {p}}{\mathrm {p}}\) interaction vertices are reconstructed from tracks in the silicon tracker using the deterministic annealing filter algorithm [46]. The \({\mathrm {p}}{\mathrm {p}}\) interaction vertex with the highest \(\sum p_{\mathrm {T}} ^{2}\) of the associated clusters of physics objects is considered to be the primary vertex associated with the hard scattering interaction. Here, the physics objects are the jets, which are clustered with the tracks assigned to the vertex using the anti-\(k_{\mathrm {T}}\) jet clustering algorithm [47, 48], and the missing transverse momentum \({\vec {p}}_{\mathrm {T}}^{\text {miss}}\), defined as the negative vector sum of the \({\vec {p}}_{\mathrm {T}}\) of those jets, with its magnitude referred to as \(p_{\mathrm {T}} ^\text {miss}\). The interaction vertices not associated with the hard scattering are designated as pileup vertices.

Electron candidates are reconstructed from clusters of energy deposited in the ECAL and from hits in the silicon tracker [49]. The clusters are first matched to track seeds in the pixel detector, then the trajectory of an electron candidate is reconstructed considering energy lost by the electron due to bremsstrahlung as it traverses the material of the tracker, using a Gaussian sum filter algorithm. The PF algorithm further distinguishes electrons from charged pions using a multivariate approach [50]. Observables related to the energy and geometrical matching between track and ECAL cluster(s) are used as main inputs. Additional requirements are applied on the ECAL shower shape, the variables related to the track-cluster matching, the impact parameter, and the ratio of the energies measured in the HCAL and ECAL in the region around the electron candidate. With these requirements, the reconstruction and identification efficiency of an electron from a \(\mathrm {Z} \rightarrow \mathrm {e}^+\mathrm {e}^- \) decay is on average 70%, whereas the misidentification rate is 1–2% [49]. Electrons with \(p_{\mathrm {T}} > 25\,\text {Ge}\text {V} \) and \(|\eta |<2.4\) are selected for this analysis. Further, electrons passing through the transition regions between the ECAL barrel and endcap sections, (\(1.444< |\eta | < 1.566\)), which are less well measured, are removed.

Muon candidates are identified by multiple reconstruction algorithms using hits in the silicon tracker and signals in the muon system. The standalone muon algorithm uses only information from the muon detectors. The tracker muon algorithm starts from tracks found in the silicon tracker and then associates them with matching tracks in the muon detectors. The global muon algorithm starts from standalone muons and then performs a global fit to consistent hits in the tracker and the muon system [51]. Global muons are used by the PF algorithm. Muons are required to pass additional identification criteria based on the track impact parameter, the quality of the track reconstruction, and the number of hits recorded in the tracker and the muon systems. Muons selected for this analysis are required to have \(p_{\mathrm {T}} > 25\,\text {Ge}\text {V} \) and \(|\eta |<2.4\).

Charged leptons (electrons or muons) from \(\mathrm {Z} \rightarrow \mathrm {e}^+\mathrm {e}^- \) or \(\mathrm {Z} \rightarrow \mathrm {\mu ^+}\mathrm {\mu ^-} \) decays, with the \(\mathrm {Z}\) boson originating from the decay of a heavy VLQ, are expected to be isolated, i.e., to have low levels of energy deposited in the calorimeter regions around their trajectories. An isolation variable is defined as the scalar \(p_{\mathrm {T}}\) sum of the charged and neutral hadrons and photons in a cone centered on the direction of the lepton, of radius \(\varDelta R \equiv \sqrt{\smash [b]{(\varDelta \eta )^2+(\varDelta \phi )^2}}\), with \(\varDelta R = 0.3\,(0.4)\) for electrons (muons). The \(p_{\mathrm {T}}\) contributions from pileup and from the lepton itself are subtracted from the isolation variable [49, 51]. The relative isolation parameter, defined as the isolation variable divided by the lepton \(p_{\mathrm {T}}\), is required to be less than 0.06 (0.15) for the electrons (muons), with corresponding efficiencies of 85 and 95%, respectively, based on simulation. The isolation requirement helps reject jets misidentified as leptons and reduce multijet backgrounds.

The anti-\(k_{\mathrm {T}}\) jet clustering algorithm [47, 48] reconstructs jets with PF candidates as inputs. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies. To suppress the contribution from pileup, charged particles not originating from the primary vertex are removed from the jet clustering. An event-by-event jet-area-based correction [52, 53] is applied to subtract the contribution of the neutral-particle component of the pileup. Residual corrections are applied to the data to account for the differences with the simulations [54].

Two types of jet are considered, distinguished by the choice of distance parameter used for clustering. Those clustered with a distance parameter of 0.4 (“AK4 jets”), are required to have \(p_{\mathrm {T}} > 30\,\text {Ge}\text {V} \), and those clustered with a value of 0.8 for this parameter (“AK8 jets”) must satisfy the condition \(p_{\mathrm {T}} > 200\,\text {Ge}\text {V} \), where the jet momentum is the vector sum of the momenta of all particles clustered in the jet. Both classes of jets must satisfy \(|\eta |<2.4\). A new value for \(p_{\mathrm {T}} ^\text {miss}\) is determined using the PF objects and including the jet energy corrections.

The combined secondary vertex \(\mathrm {b}\) tagging algorithm (CSVv2) [55] is used to identify jets originating from the hadronization of \(\mathrm {b}\) quarks. The algorithm combines information on tracks from the silicon tracker and vertices associated with the jets using a multivariate discriminant. An AK4 jet is defined as a \(\mathrm {b}\)-tagged jet if the corresponding CSVv2 discriminant is above a threshold that gives an average efficiency of about 70% for \(\mathrm {b}\) quark jets and a misidentification rate of 1% for light-flavored jets.

The signal events searched for in this analysis have two massive VLQs decaying to at least one \(\mathrm {Z}\) boson and either a \(\mathrm {Z}\), \(\mathrm {W}\), or Higgs boson and two heavy quarks. One \(\mathrm {Z}\) boson must decay leptonically, whereas the remaining \(\mathrm {Z}\), \(\mathrm {W}\), or Higgs boson is reconstructed using its hadronic decays into jets. Depending on the mass of the VLQ, the decay products can have a large Lorentz boost. In this case, the decay products of \(\mathrm {W}\rightarrow \mathrm {q} \overline{\mathrm {q}} ^\prime \) and \(\mathrm {Z} \rightarrow \mathrm {q} \overline{\mathrm {q}} \) (collectively labeled as \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \)), \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \), and \(\mathrm {t}\rightarrow \mathrm {q} \overline{\mathrm {q}} ^\prime \mathrm {b}\) may be contained within a single AK8 jet. These decays are reconstructed using a jet substructure tagger. The decay products of heavy bosons and top quarks that do not acquire a large Lorentz boost are identified by a resolved tagger using AK4 jets. Both types of taggers are described in the next section.

4 Event selection and categorization

For the dielectron (\(\mathrm {Z} \rightarrow \mathrm {e}^+\mathrm {e}^- \)) channel, event candidates are selected using triggers requiring the presence of at least one electron with \(p_{\mathrm {T}} >115\,\text {Ge}\text {V} \) or a photon with \(p_{\mathrm {T}} >175\,\text {Ge}\text {V} \). After passing one of the triggers, the triggering electron is also required to pass a set of criteria based on the electromagnetic shower shape and the quality of the electron track. A loose isolation criterion on the electrons is further required, as described in Sect. 3. One of the electrons is required to have \(p_{\mathrm {T}} >120\,\text {Ge}\text {V} \) in order to remain above the triggering electron \(p_{\mathrm {T}}\) threshold. Since the signal electrons originate from the decay of highly boosted \(\mathrm {Z}\) bosons, these selection criteria preserve the high signal efficiency, while reducing the number of misidentified electrons. The photon trigger helps to retain electrons with \(p_{\mathrm {T}} >300\,\text {Ge}\text {V} \) that would otherwise be lost because of the requirements on electromagnetic shower shape in the ECAL.

For the dimuon (\(\mathrm {Z} \rightarrow \mathrm {\mu ^+}\mathrm {\mu ^-} \)) channel, event candidates are selected using a trigger that requires presence of at least one muon with \(p_{\mathrm {T}} >24\,\text {Ge}\text {V} \). The trigger implements a loose isolation requirement by allowing only a small energy deposit in the calorimeters around the muon trajectory. After passing the trigger, one of the muons from the \(\mathrm {Z} \rightarrow \mathrm {\mu ^+}\mathrm {\mu ^-} \) decay must have \(p_{\mathrm {T}} >45\,\text {Ge}\text {V} \), which provides the largest background rejection that can be obtained without decreasing the signal efficiency for the VLQ mass range of interest. The trigger and lepton reconstruction and identification efficiencies are determined using a tag-and-probe method [56]. Scale factors are applied to the simulated events to account for any efficiency differences between the data and simulation.

The invariant mass of the lepton pair from the \(\mathrm {Z}\) boson leptonic decay must satisfy \(75< m(\ell \ell ) < 105\,\text {Ge}\text {V} \), to be consistent with the \(\mathrm {Z}\) boson mass, and have a total \(p_{\mathrm {T}} (\ell \ell ) > 100\,\text {Ge}\text {V} \), appropriate for the decay of a massive VLQ. Events must have exactly one \(\mathrm {e}^+\mathrm {e}^- \) or \(\mathrm {\mu ^+}\mathrm {\mu ^-} \) pair candidate consistent with a \(\mathrm {Z}\) boson decay.

Events are required to have at least three AK4 jets with \(H_{\mathrm {T}} > 200\,\text {Ge}\text {V} \), and \(H_{\mathrm {T}} \equiv \sum p_{\mathrm {T}} \), where the summation is over all jets in the event. The highest \(p_{\mathrm {T}}\) (leading) AK4 jet is required to have \(p_{\mathrm {T}} > 100\,\text {Ge}\text {V} \), the second-highest-\(p_{\mathrm {T}}\) (subleading) AK4 jet to have \(p_{\mathrm {T}} > 50\,\text {Ge}\text {V} \), and all other jets must satisfy the condition \(p_{\mathrm {T}} > 30\,\text {Ge}\text {V} \). The AK4 (AK8) jets j within \(\varDelta R(\ell , j) < 0.4\) (0.8) of either lepton from the \(\mathrm {Z}\) boson decay are not considered further in the analysis. At least one \(\mathrm {b}\)-tagged jet with \(p_{\mathrm {T}} > 50\,\text {Ge}\text {V} \) is required. The \(S_{\mathrm {T}}\) variable, defined as the sum of \(H_{\mathrm {T}}\), \(p_{\mathrm {T}} (\mathrm {Z})\), and \(p_{\mathrm {T}} ^\text {miss}\), must be greater than 1000\(\,\text {Ge}\text {V}\). The selection criteria are summarized in Table 1. The selections are optimized to obtain the largest suppression of SM backgrounds that can be achieved without reducing the simulated signal efficiency by more than 1%.

Table 1 Event selection criteria

The event topologies are different for \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) decays, and the product of the signal efficiency and the acceptance varies from 1.2 to 2.6% over the various signal channels. The \(\mathrm {T} \overline{\mathrm {T}} \) events are characterized by three heavy bosons and two heavy quarks in the decay sequence. The \(\mathrm {B} \overline{\mathrm {B}} \) events have two heavy bosons and two heavy quarks, hence more energetic final decay objects. Therefore, the analysis is optimized separately for the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) channels.

For both searches the decays of boosted \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \) and \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \) are reconstructed from AK8 jets, using the jet substructure tagger, and are referred to as \(\mathrm {V}\) and \(\mathrm {H}\) jets, respectively. As the Higgs boson mass is larger than \(\mathrm {W}\) and \(\mathrm {Z}\) boson masses, it requires a higher momentum for its decay products to merge into a single AK8 jet. Therefore, \(\mathrm {H}\) jets are required to have \(p_{\mathrm {T}} > 300\,\text {Ge}\text {V} \) and \(\mathrm {V}\) jets have \(p_{\mathrm {T}} > 200\,\text {Ge}\text {V} \). A jet pruning algorithm [57, 58] is used to measure the jet mass. The \(\mathrm {V}\) and \(\mathrm {H}\) jet candidates are required to have a pruned jet mass in the range 65–105 and 105–135\(\,\text {Ge}\text {V}\), respectively. The jet pruning algorithm reclusters the groomed jets [59] by eliminating low energy subjets subjets. In the subsequent recombination of two subjets, the ratio of the subleading subjet \(p_{\mathrm {T}}\) to the pruned jet \(p_{\mathrm {T}}\) must be greater than 0.1 and the distance between the two subjets must satisfy \(\varDelta R < m_\mathrm {jet}/2{p_{\mathrm {T}}}_\mathrm {jet}\), where \(m_\mathrm {jet}\) and \({p_{\mathrm {T}}}_\mathrm {jet}\) are the mass and \(p_{\mathrm {T}}\) of the pruned jet, respectively.

The N-subjettiness algorithm [60] is used to calculate the jet shape variable \(\tau _{N}\), which quantifies the consistency of a jet with the hypothesis of the jet having N subjets, each arising from a hard parton coming from the decay of an original heavy boson. The \(\mathrm {V}\) and \(\mathrm {H}\) jets in the \(\mathrm {T} \overline{\mathrm {T}} \) (\(\mathrm {B} \overline{\mathrm {B}} \)) search are required to have an N-subjettiness ratio \(\tau _{21} \equiv \tau _{2}/\tau _{1} < 1.0\,(0.6)\). Both pruned subjets coming from the \(\mathrm {H}\) jet are required to be \(\mathrm {b}\)-tagged. This is done by using the above-mentioned CSVv2 b-tagging algorithm with a cut that gives a 70–90% efficiency for \(\mathrm {b}\) quark subjets, depending on the subjet \(p_{\mathrm {T}}\), and a misidentification rate of 10% for subjets from light-flavored quarks and gluons.

Boosted top quarks decaying to \(\mathrm {b}\mathrm {q} \overline{\mathrm {q}} ^\prime \) are identified (“\(\mathrm {t}\) tagged”) using AK8 jets and the soft-drop algorithm [61, 62] to groom the jet. This algorithm recursively declusters a jet into two subjets. It discards soft and wide-angle radiative jet components until a hard-splitting criterion is met, to obtain jets consistent with the decay of a massive particle. We use the algorithm with an angular exponent \(\beta = 0\), a soft cutoff threshold \(z_{cut} < 0.1\), and a characteristic radius \(R_0 = 0.8\). For top quark jets, the soft-drop mass must be in the range 105–220\(\,\text {Ge}\text {V}\) and the N-subjettiness ratio \(\tau _{32} \equiv \tau _{3}/\tau _{2} < 0.81\,(0.67)\) for the \(\mathrm {T} \overline{\mathrm {T}} \) (\(\mathrm {B} \overline{\mathrm {B}} \)) search, consistent with the expectation for three subjets from top quark decay. There are a total of five heavy bosons and quarks produced in \(\mathrm {T} \overline{\mathrm {T}} \) signal events, whereas there are only four in \(\mathrm {B} \overline{\mathrm {B}} \) events. Thus it is possible to apply a tighter N-subjetiness ratio criterion in the \(\mathrm {B} \overline{\mathrm {B}} \) analysis without a loss of signal efficiency.

Corrections to the jet mass scale, resolution and \(\tau _{21} \) selection efficiency for \(\mathrm {V}\) jets due to the difference in data and MC simulation are measured using a sample of semileptonic \({\mathrm {t}\overline{\mathrm {t}}}\) events [63]. For the correction to the jet mass scale and resolution, boosted \(\mathrm {W}\) bosons produced in the top quark decays are separated from the combinatorial \({\mathrm {t}\overline{\mathrm {t}}}\) background by performing a simultaneous fit to the observed pruned jet mass spectrum. In order to account for the difference in the jet shower profile of \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \) and \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \) decays, a correction factor to the \(\mathrm {H}\) jets mass scale and resolution [64] is measured by comparing the ratio of \(\mathrm {H} \) and \(\mathrm {V} \) jet efficiencies using the pythia  8.212 [28] and herwig ++ [65] shower generators. In addition, the corrections to \(\tau _{21} \) selection efficiency are obtained based on the difference between data and simulation [64] for \(\mathrm {H}\)-tagged jets. All these corrections are propagated to \(\mathrm {V}\), top quark and \(\mathrm {H}\) jets, respectively. For top quark jets, the corrections to the \(\tau _{32} \) selection efficiency are measured between data and simulation [63] using soft-drop groomed jets. To account for the misidentification of boosted \(\mathrm {V}\)-, \(\mathrm {H}\)-, and \(\mathrm {t}\)-tagged jets in the background samples, mistagging scale factors are derived from a region in the data enriched in \(\mathrm {Z} \)+jets events, which is constructed using the selection criteria listed in Table 1, with the exception that events must have zero \(\mathrm {b}\) jets. These mistagging scale factors are applied to the mistagged jets in simulated signal and background events.

In the \(\mathrm {T} \overline{\mathrm {T}} \) search, in addition to the jet substructure techniques, the \(\mathrm {W}\), \(\mathrm {Z}\), \(\mathrm {H}\), and top quark decays are reconstructed with a resolved tagger using AK4 jets, as described below. Only those AK4 jets that are a radial distance \(\varDelta R > 0.8\) from the tagged AK8 jets are considered in the resolved tagging algorithm. The resolved \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \) and \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \) candidates are composed of two AK4 jets \(j_{1}\) and \(j_{2}\) whose invariant mass must satisfy \(70< m(j_{1} j_{2}) < 120\,\text {Ge}\text {V} \) and \(80< m(j_{1} j_{2}) < 160\,\text {Ge}\text {V} \), respectively, and have \(p_{\mathrm {T}} (j_{1} j_{2}) > 100\,\text {Ge}\text {V} \). For \(\mathrm {H}\) candidates, at least one of the jets must be \(\mathrm {b}\) tagged. The resolved top quark candidate is composed of either three AK4 jets \(j_{1}\), \(j_{2}\), and \(j_{3}\) with an invariant mass \(120< m(j_{1} j_{2} j_{3}) < 240\,\text {Ge}\text {V} \) and \(p_{\mathrm {T}} (j_{1} j_{2} j_{3}) > 100\,\text {Ge}\text {V} \), or an AK4 jet \(j_{1}\) and an AK8 \(\mathrm {V}\) jet satisfying \(120< m(\mathrm {V} j_{1}) < 240\,\text {Ge}\text {V} \) and \(p_{\mathrm {T}} (\mathrm {V} j_{1}) > 150\,\text {Ge}\text {V} \). These selection criteria are derived from simulated \(\mathrm {T} \overline{\mathrm {T}} \) events, using MC truth information.

The \(\mathrm {T} \overline{\mathrm {T}} \) events are next classified based on the number of AK4 \(\mathrm {b}\)-tagged jets (\(N_{\mathrm {b}}\)), and number of \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \) (\(N_{\mathrm {V}}\)), \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \) (\(N_{\mathrm {H}}\)), and \(\mathrm {t}\rightarrow \mathrm {q} \overline{\mathrm {q}} ^\prime \mathrm {b}\) (\(N_{\mathrm {t}}\)) candidates identified using either the jet substructure or resolved tagging algorithms. In an event, \(N_{\mathrm {b}}\) can be 1 or \({\ge }2\), and \(N_{\mathrm {V}}\), \(N_{\mathrm {H}}\), and \(N_{\mathrm {t}}\) each can be 0 or \({\ge }1\). Thus, in total, \(2{\times }2{\times }2{\times }2=16\) categories of events are constructed. For simplicity, overlaps between candidates of different types are allowed, e.g., the same AK8 jet could be tagged as both a top quark and an \(\mathrm {H}\) candidate because of the overlapping mass windows. Such overlaps occur in a few percent of the signal events. However, by construction each event can belong to only one category. In the example above, the event would fall into a category with both \(N_{\mathrm {H}} \ge 1\) and \(N_{\mathrm {t}} \ge 1\) requirements satisfied. Further, the mistag rates and the relevant corrections to the jet mass scale and resolution are applied to the \(\mathrm {H}\) and \(\mathrm {t}\) candidates, based on MC truth information.

Table 2 The first four columns show different event groups used for the \(\mathrm {T} \) \(\overline{\mathrm {T}}\) search, classified according to the number of \(\mathrm {b}\)-tagged jets \(N_{\mathrm {b}}\) and the number of \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \), \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \), and \(\mathrm {t}\rightarrow \mathrm {q} \overline{\mathrm {q}} ^\prime \mathrm {b}\) candidates in the event, \(N_{\mathrm {V}}\), \(N_{\mathrm {H}}\) and \(N_{\mathrm {t}}\), respectively, identified using both the jet substructure and resolved tagger algorithms. The last three columns show the relative signal acceptance for a \(\mathrm {T}\) quark of mass 1200\(\,\text {Ge}\text {V}\) for decay channels tZtZ, tZtH and tZbW as described in text

Next, the event categories are sorted using the figure of merit \(S/\sqrt{B}\), where S and B are the expected \(\mathrm {T} \overline{\mathrm {T}} \rightarrow \mathrm {t}\mathrm {Z}\mathrm {t}\mathrm {Z}\) signal and background event yields, respectively, as determined from the simulation. The categories with similar figures of merit based on expected upper limits at 95% confidence level (\(\text {CL}\)) are grouped together, while the categories that are found not to add sensitivity to the \(\mathrm {T} \overline{\mathrm {T}} \) signal are discarded. A total of four event groups labeled A through D are selected, each with a different signal acceptance relative to the selection criteria described in Table 1 and depending on the \(\mathrm {T}\) decay channel. Table 2 shows the selections on these event groups, and the relative signal acceptances of the \(\mathrm {T}\) quark decay channels, namely, tZtZ, tZtH, or tZbW for a \(\mathrm {T}\) quark of mass 1200\(\,\text {Ge}\text {V}\). The decay channels are defined with a benchmark combination of branching fractions \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = 100\%\) (tZtZ), \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) = 50\%\) (tZtH), and \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 50\%\) (tZbW). Events from all the decay channels mainly contribute to groups A and B, whereas groups C and D have slightly lower acceptance depending on the decay channel. The fraction of the signal identified by the jet substructure and resolved taggers depends on the \(\mathrm {T}\) quark mass. For masses below 1200\(\,\text {Ge}\text {V}\), the two taggers are equally efficient in identifying signal events for all the channels. For \(\mathrm {T}\) quark masses above 1200\(\,\text {Ge}\text {V}\), the jet substructure tagger becomes more efficient. For example, for \(\mathrm {T}\) quark mass at 1800\(\,\text {Ge}\text {V}\), the jet substructure tagger selects twice as many \(\mathrm {T}\) quark candidates as the resolved tagger.

Table 3 The first four columns show different event categories used for the \(\mathrm {B} \overline{\mathrm {B}} \) search, classified according to the number of AK4 \(\mathrm {b}\)-tagged jets \(N_{\mathrm {b}}\) and the number of \(\mathrm {V} \rightarrow \mathrm {q} \overline{\mathrm {q}} \), \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \), and \(\mathrm {t}\rightarrow \mathrm {q} \overline{\mathrm {q}} ^\prime \mathrm {b}\) candidates in the event, \(N_{\mathrm {V}}\), \(N_{\mathrm {H}}\), and \(N_{\mathrm {t}}\), respectively, identified using the jet substructure algorithm. The last three columns show the relative signal acceptance for a \(\mathrm {B}\) quark of mass 1200\(\,\text {Ge}\text {V}\) for decay channels bZbZ, bZbH and bZtW as described in text

Because the event topology of \(\mathrm {B} \overline{\mathrm {B}} \) signal events is different from that of \(\mathrm {T} \overline{\mathrm {T}} \) signal events, as discussed previously, the \(\mathrm {V}\), \(\mathrm {H}\), and \(\mathrm {t}\) candidates in the \(\mathrm {B} \overline{\mathrm {B}} \) analysis are identified using only the jet substructure tagger. Events are then separated into five categories, labeled 1\(\mathrm {b}\), 2\(\mathrm {b}\), boosted \(\mathrm {t}\), boosted \(\mathrm {H}\), and boosted \(\mathrm {Z}\), based on the values of \(N_{\mathrm {b}}\), \(N_{\mathrm {V}}\), \(N_{\mathrm {H}}\), and \(N_{\mathrm {t}}\). Table 3 shows these categories, and the relative signal acceptances of \(\mathrm {B}\) quark decay channels, namely, bZbZ, bZbH, or bZtW for a \(\mathrm {B}\) quark of mass 1200\(\,\text {Ge}\text {V}\). The decay channels are defined with a benchmark combination of branching fractions \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = 100\%\) (bZbZ), \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) = 50\%\) (bZbH), and \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 50\%\) (bZtW).

5 Background modeling

The backgrounds from all sources are estimated using simulation, except for \(\mathrm {Z} \)+jets where corrections to the simulated events are applied using data, as described below. The modeling of simulated background events is validated using several control regions in the data, which are constructed by inverting one or more of the requirements listed in Table 1. The control region labeled CR0b+high-\(S_{\mathrm {T}}\) is constructed by requiring zero \(\mathrm {b}\) jets. The control region CR1\(\mathrm {b}\)+low-\(S_{\mathrm {T}}\) is constructed by inverting the \(S_{\mathrm {T}}\) requirement: \(S_{\mathrm {T}} \le 1000\,\text {Ge}\text {V} \). The control region CR0b is constructed by requiring zero \(\mathrm {b}\) jets and removing the \(S_{\mathrm {T}}\) requirement. Signal contamination from all channels in each of these control regions is less than 1%.

The AK4 jet multiplicity distribution is not modeled reliably in the \(\mathrm {Z} \)+jets simulation, and therefore it is corrected using scale factors obtained from data. Scale factors listed in Table 4 are determined using the CR0b control region, which is enriched with \(\mathrm {Z} \)+jets events. After applying these corrections, the distributions of kinematic variables in the control regions from the background simulations are in agreement with the data, as shown for example in Fig. 2 for the \(S_{\mathrm {T}}\) distributions.

Table 4 The scale factors determined from data for correcting the AK4 jet multiplicity distribution in the simulation. The quoted uncertainties in the scale factors are statistical only
Fig. 2
figure 2

The \(S_{\mathrm {T}}\) distributions for the CR1\(\mathrm {b}\)+low-\(S_{\mathrm {T}}\) (left) and CR0\(\mathrm {b}\)+high-\(S_{\mathrm {T}}\) (right) control regions for the data (points) and the background simulations (shaded histograms) after applying the scale factors given in Table 4. The vertical bars on the points represent the statistical uncertainties in the data. The hatched bands indicate the total uncertainties in the simulated background contributions added in quadrature. The lower plots show the difference between the data and the simulated background, divided by the total uncertainty

6 Systematic uncertainties

Table 5 Summary of systematic uncertainties considered in the statistical analysis of \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) search on the background and signal events. All uncertainties affect the normalizations of the \(S_{\mathrm {T}}\) distributions. The tick mark indicates the uncertainties that also affect the shape, and the uncertainty range accounts for their effects on the expected yields across all the \(\mathrm {T} \overline{\mathrm {T}} \) groups or \(\mathrm {B} \overline{\mathrm {B}} \) categories. The \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) signal events correspond to the benchmark decay channels tZtZ and bZbZ, respectively, for \(\mathrm {T}\) and \(\mathrm {B}\) quark mass \(m_{\mathrm {T}}=m_{\mathrm {B}}=1200\,\text {Ge}\text {V} \)

The systematic uncertainties in the SM background rates are due to the uncertainties in the CMS measurements of \(\mathrm {d}\sigma /\mathrm {d}H_{\mathrm {T}} \) for \(\mathrm {Z} \)+jets  [66], \(\mathrm {d}\sigma /\mathrm {d}m_{{\mathrm {t}\overline{\mathrm {t}}}}\) for \(\mathrm {t}\overline{\mathrm {t}}\)+jets  [67], and \(\mathrm {d}\sigma /\mathrm {d}p_{\mathrm {T}} (\mathrm {Z})\) for diboson production [68]. They are estimated to be 15% in each case. The measured integrated luminosity uncertainty of 2.5% [69] affects both the signal and background rate predictions. The uncertainties associated with the measured data-to-simulation efficiency scale factors for the lepton identification and the trigger efficiencies are 3 and 1%, respectively.

The effect on the signal and background acceptance uncertainties due to the renormalization and factorization scale (\(\mu _f\) and \(\mu _r\)) uncertainties and the PDF choices in the simulations are taken into account in the statistical analysis. The influence of \(\mu _f\) and \(\mu _r\) scale uncertainties are estimated by varying the default scales by the following six combinations of factors, (\(\mu _f\), \(\mu _r\)) \(\times \) (1/2, 1/2), (1/2, 1), (1, 1/2), (2, 2), (2, 1), and (1, 2). The maximum and minimum of the six variations are computed for each bin of the \(S_{\mathrm {T}}\) distribution, producing an uncertainty “envelope”. The uncertainties due to the PDF choices in the simulations are estimated using the PDF4LHC procedure [27, 70,71,72], where the root-mean-square of 100 pseudo-experiments provided by the PDF sets represents the uncertainty envelope. The background and signal event counts are then varied relative to their nominal values up and down by a factor of two times the uncertainty envelopes. The impacts of these variations on the background and signal shape are also taken into account. The effect of the \(\mu _f\) and \(\mu _r\) scale uncertainties on the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) signal yield is \(<1\)%. However, this has the largest effect, amounting to as much as 36% on the background yield. The effect due to PDF choices amounts to a 3.2–9.5% change in the signal and background yields. The effect of the uncertainty in the pileup determination is estimated by varying the nominal \({\mathrm {p}}{\mathrm {p}}\) inelastic cross section by 4.6% [42], which has an impact of 1.5–3.6% on the signal yields. Differences between simulation and data in the jet multiplicity distributions in DY+jets background events, derived in the CR0b region as shown in Table 4, are taken as an estimate of the associated systematic uncertainty, which ranges from 4.0–11.5%

Several uncertainties are associated with the measurement of jet-related quantities. The jet energy scale and resolution uncertainties are about 1% [54, 73]. The AK8 pruned jet mass scale and resolution uncertainties are evaluated to be 2.3 and 18% [63], respectively. The effect of these uncertainties on the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) signal yields is 1.5–4.4% and 1.0–3.8%, respectively. These uncertainties, in addition to the uncertainties in the \(\tau _{21}\) (8%) and \(\tau _{32}\) (11%) selections [63], are applied for the \(\mathrm {V}\)-, \(\mathrm {H}\)-, and \(\mathrm {t}\)-tagged jets. The systematic uncertainties due to the jet shower profile differences between the jets in the \(\mathrm {W}\rightarrow \mathrm {q} \overline{\mathrm {q}} ^\prime \) and \(\mathrm {H} \rightarrow \mathrm {b} \overline{\mathrm {b}} \) processes are estimated from the difference observed between results obtained with the pythia  8 and herwig ++ generators and are applied to the \(\mathrm {V}\)- and \(\mathrm {H}\)-tagged jets. The overall effect of \(\mathrm {V}\), \(\mathrm {H}\), and \(\mathrm {t}\) tagging uncertainties on \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) signal yields is 0.2–8.4%. The uncertainties in the misidentification rates of boosted jets are 5, 14, and 7% for the \(\mathrm {W}\)-, \(\mathrm {H}\)-, and \(\mathrm {t}\)-tagged jets, respectively. They are used to derive the uncertainties in the estimates of the numbers of mistagged jets in the signal and background simulated events, which result in uncertainties in the \(\mathrm {B} \overline{\mathrm {B}} \) signal yields of up to 14%. The uncertainties in the \(\mathrm {b}\) tagging efficiency scale factors are propagated to the final result, with the uncertainties in the \(\mathrm {b}\)- and \(\mathrm {c}\)-flavored quark jets treated as fully correlated. These uncertainties are in the range 2–5% for \(\mathrm {b}\)-flavored jets, a factor of two larger for \(\mathrm {c}\)-flavored jets, and \(\approx \)10% for light-flavored jets. The uncertainties due to heavy- and light-flavored jets are considered uncorrelated. Table 5 summarizes the systematic uncertainties in the background and signal yields in the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) searches. The ranges correspond to the impact on event yields due to systematic uncertainties that affect both the rates and shapes across all the \(\mathrm {T} \overline{\mathrm {T}} \) groups or \(\mathrm {B} \overline{\mathrm {B}} \) categories. Here the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) signals correspond to the benchmark decay channels tZtZ and bZbZ, respectively, for a \(\mathrm {T}\) and \(\mathrm {B}\) quark mass \(m_{\mathrm {T}}=m_{\mathrm {B}}=1200\,\text {Ge}\text {V} \).

7 Results

7.1 \(\mathrm {T}\) quark search

The number of observed events for the \(\mathrm {T} \overline{\mathrm {T}} \) production search in the A, B, C, and D event groups are given for the electron and muon channels in Tables 6 and 7, respectively, along with the numbers of predicted background events. The expected numbers of signal events for \(\mathrm {T}\) quark masses of 800 and 1200\(\,\text {Ge}\text {V}\) are also shown in the same tables, for three different decay scenarios, with branching fractions \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = 100\%\) (tZtZ), \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) = 50\%\) (tZtH), and \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 50\%\) (tZbW). The predicted background and observed event yields agree within their uncertainties.

Table 6 The number of observed events and the predicted number of SM background events in the \(\mathrm {T} \overline{\mathrm {T}} \) search using \(\mathrm {Z} \rightarrow \mathrm {e}^+\mathrm {e}^- \) channel in the four event groups. The expected numbers of signal events for \(\mathrm {T}\) quark masses of 800 and 1200\(\,\text {Ge}\text {V}\) for three different decay scenarios with assumed branching fractions \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = 100\%\) (tZtZ) , \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) = 50\%\) (tZtH), and \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 50\%\) (tZbW) are also shown. The uncertainties in the number of expected background events include the statistical and systematic uncertainties added in quadrature
Table 7 The number of observed events and the predicted number of SM background events in the \(\mathrm {T} \overline{\mathrm {T}} \) search using \(\mathrm {Z} \rightarrow \mathrm {\mu ^+}\mathrm {\mu ^-} \) channel in the four event groups. The expected numbers of signal events for \(\mathrm {T}\) quark masses of 800 and 1200\(\,\text {Ge}\text {V}\) for three different decay scenarios with assumed branching fractions \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = 100\%\) (tZtZ) , \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) = 50\%\) (tZtH), and \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 50\%\) (tZbW) are also shown. The uncertainties in the number of expected background events include the statistical and systematic uncertainties added in quadrature

To determine the upper limits on the \(\mathrm {T} \overline{\mathrm {T}} \) cross section, the electron and muon channels are combined, and a simultaneous binned maximum-likelihood fit is performed on the \(S_{\mathrm {T}}\) distributions in data for the four event groups. The measured \(S_{\mathrm {T}}\) distributions in data are shown in Fig. 3 for each of the event groups, along with the predicted background distributions and the expected signal distributions for \(\mathrm {T} \overline{\mathrm {T}} \rightarrow \mathrm {t}\mathrm {Z}\mathrm {t}\mathrm {Z}\) with \(m_{\mathrm {T}}=1200\,\text {Ge}\text {V} \). The impact of the statistical uncertainty in the simulated samples is reduced by rebinning each \(S_{\mathrm {T}}\) distribution to ensure that the statistical uncertainty associated with the expected background is less than 20% in each bin. There is no indication of a signal in the \(S_{\mathrm {T}}\) distribution of any of the event groups.

Fig. 3
figure 3

The \(S_{\mathrm {T}}\) distributions for groups A, B, C, D (left to right, upper to lower) from data (points with vertical and horizontal bars), the expected SM backgrounds (shaded histograms), and the expected signal, scaled up by a factor 2, for \(\mathrm {T} \overline{\mathrm {T}} \rightarrow \mathrm {t}\mathrm {Z}\mathrm {t}\mathrm {Z}\) with \(m_{\mathrm {T}}=1200\,\text {Ge}\text {V} \) (dotted lines). The vertical bars on the points show the central 68% \(\text {CL}\) intervals for Poisson-distributed data. The horizontal bars give the bin widths. The hatched bands represent the statistical and systematic uncertainties in the total background contribution added in quadrature. The lower plots give the difference between the data and the total expected background, divided by the total background uncertainty

The upper limits at 95% \(\text {CL}\) on the \(\mathrm {T} \overline{\mathrm {T}} \) cross section are computed using a Bayesian likelihood-based technique [74] with the Theta framework [75]. All the systematic uncertainties due to normalization variations described in the previous section enter the likelihood as nuisance parameters with log-normal prior distributions, whereas the uncertainties from the shape variations are assigned Gaussian-distributed priors. For the signal cross section parameter, we use a uniform prior distribution. The likelihood is marginalized with respect to the nuisance parameters, and the limits are extracted from a simultaneous maximum-likelihood fit of the \(S_{\mathrm {T}}\) distributions in all four groups shown in Fig. 3.

The upper limits on the \(\mathrm {T} \overline{\mathrm {T}} \) cross section are computed for different \(\mathrm {T}\) quark mass values and for the three branching fraction scenarios listed above. The upper limits at 95% \(\text {CL}\) on the \(\mathrm {T} \overline{\mathrm {T}} \) cross section are shown as a function of the \(\mathrm {T}\) quark mass by the solid line in Fig. 4. The median expected upper limit is given by the dotted line, while the inner and outer bands correspond to one and two standard deviation uncertainties, respectively, in the expected limit. The dotted-dashed curve displays the predicted theoretical signal cross section [30]. Comparing the observed cross section limits to the theoretical signal cross section, we exclude \(\mathrm {T}\) quarks with masses less than 1280, 1185, and 1120\(\,\text {Ge}\text {V}\), respectively, for the three branching ratio hypotheses listed above. The expected upper limits are 1290, 1175, and 1115\(\,\text {Ge}\text {V}\) for the respective scenarios.

Fig. 4
figure 4

The observed (solid line) and expected (dashed line) 95% \(\text {CL}\) upper limits on the \(\mathrm {T} \overline{\mathrm {T}} \) cross section as a function of the \(\mathrm {T}\) quark mass assuming (upper left) \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = 100\%\), (upper right) \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) = 50\%\), and (lower) \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) = \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 50\%\). The dotted-dashed curve displays the theoretical \(\mathrm {T} \overline{\mathrm {T}} \) production cross section. The inner and outer bands show the one and two standard deviation uncertainties in the expected limits, respectively

Figure 5 (upper) displays the observed (left) and expected (right) 95% \(\text {CL}\) lower limits on the \(\mathrm {T}\) quark mass as a function of the relevant branching fractions, assuming \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) + \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) + \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 1.0\). For a \(\mathrm {T}\) quark decaying exclusively via \(\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}\), the lower mass limit is 1280\(\,\text {Ge}\text {V}\).

Fig. 5
figure 5

The observed (left) and expected (right) 95% \(\text {CL}\) lower limits on the mass of the \(\mathrm {T}\) (upper) and \(\mathrm {B}\) (lower) quark, in \(\,\text {Ge}\text {V}\), for various branching fraction scenarios, assuming \(\mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}) + \mathcal {B} (\mathrm {T} \rightarrow \mathrm {t}\mathrm {H} ) + \mathcal {B} (\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}) = 1\) and \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) + \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) + \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 1\), respectively

7.2 \(\mathrm {B}\) quark search

The numbers of observed and predicted background events in the five event categories for the \(\mathrm {B} \overline{\mathrm {B}} \) search using \(\mathrm {Z} \rightarrow \mathrm {e}^+\mathrm {e}^- \) and \(\mathrm {Z} \rightarrow \mathrm {\mu ^+}\mathrm {\mu ^-} \) are given in Tables 8 and 9, respectively. The expected number of signal events in each category is also shown for \(\mathrm {B}\) masses of 800 and 1200\(\,\text {Ge}\text {V}\). The branching fraction hypotheses assumed for the three decay channels are \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = 100\%\) (bZbZ), \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) = 50\%\) (bZbH), and \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 50\%\) (bZtW). The numbers of observed and expected background events are consistent with each other for every event category. As with the \(\mathrm {T} \overline{\mathrm {T}} \) search, 95% \(\text {CL}\) upper limits on the \(\mathrm {B} \overline{\mathrm {B}} \) production cross section are determined using a simultaneous binned maximum-likelihood fit to the \(S_{\mathrm {T}}\) distributions for the different event categories, shown in Fig. 6.

Table 8 The numbers of observed events and the predicted number of SM background events in the \(\mathrm {B} \overline{\mathrm {B}} \) search for the five event categories using \(\mathrm {Z} \rightarrow \mathrm {e}^+\mathrm {e}^- \) channel. The expected numbers of signal events for \(\mathrm {B}\) masses of 800 and 1200\(\,\text {Ge}\text {V}\) with branching fraction hypotheses for the three decay channels, \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = 100\%\) (bZbZ), \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) = 50\%\) (bZbH), and \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 50\%\) (bZtW) are also shown. The uncertainties in the number of expected background events include the statistical and systematic uncertainties added in quadrature
Table 9 The number of observed events and the predicted number of SM background events in the \(\mathrm {B} \overline{\mathrm {B}} \) search for the five event categories using \(\mathrm {Z} \rightarrow \mathrm {\mu ^+}\mathrm {\mu ^-} \) channel. The expected numbers of signal events for \(\mathrm {B}\) masses of 800 and 1200\(\,\text {Ge}\text {V}\) with branching fraction hypotheses for the three decay channels, \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = 100\%\) (bZbZ), \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) = 50\%\) (bZbH), and \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 50\%\) (bZtW) are also shown. The uncertainties in the number of expected background events include the statistical and systematic uncertainties added in quadrature
Fig. 6
figure 6

The \(S_{\mathrm {T}}\) distributions for the 1\(\mathrm {b}\), 2\(\mathrm {b}\), boosted \(\mathrm {t}\), boosted \(\mathrm {H}\) and boosted \(\mathrm {Z}\) (left to right, upper to lower) event categories for the data (points with vertical and horizontal bars), and the expected background (shaded histograms). The vertical bars give the statistical uncertainty in the data, and the horizontal bars show the bin widths. The expected signal for \(\mathrm {B} \overline{\mathrm {B}} \rightarrow \mathrm {b}\mathrm {Z}\mathrm {b}\mathrm {Z} \) with \(m_{\mathrm {B}} = 1200\,\text {Ge}\text {V} \) multiplied by a factor of 5 is shown by the dashed line. The statistical and systematic uncertainties in the SM background prediction, added in quadrature, are represented by the hatched bands. The lower panel in each plot show the difference between the data and the expected background, divided by the total uncertainty

The upper limits at 95% \(\text {CL}\) on the \(\mathrm {B} \overline{\mathrm {B}} \) cross section are shown by the solid line in Fig. 7. As before, the inner and outer bands give the one and two standard deviation uncertainties, respectively, in the expected upper limits. The dotted curve displays the theoretical signal cross section. Comparing the observed cross section limits to the signal cross section, we exclude \(\mathrm {B}\) quarks with masses less than 1130, 1015, and 975\(\,\text {Ge}\text {V}\) in the bZbZ, bZbH, and bZtW channels, respectively. The corresponding expected values are 1200, 1085, and 1055\(\,\text {Ge}\text {V}\).

Fig. 7
figure 7

The observed (solid line) and expected (dashed line) 95% \(\text {CL}\) upper limits on the \(\mathrm {B} \overline{\mathrm {B}} \) production cross section versus the \(\mathrm {B}\) quark mass for (upper left) \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = 100\%\), (upper right) \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) = 50\%\), and (lower) \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) = \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 50\%\). The dotted-dashed line displays the theoretical cross section. The inner and outer bands show the one and two standard deviation uncertainties in the expected limits, respectively

Figure 5 (lower) displays the observed (left) and expected (right) 95% \(\text {CL}\) lower limits on the \(\mathrm {B}\) quark mass as a function of the relevant branching fractions, assuming \(\mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}) + \mathcal {B} (\mathrm {B} \rightarrow \mathrm {b}\mathrm {H} ) + \mathcal {B} (\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}) = 1.0\). For a \(\mathrm {B}\) quark decaying exclusively via \(\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}\), the lower mass limit is 1130\(\,\text {Ge}\text {V}\).

8 Summary

The results of a search have been presented for the pair production of vector-like top (\(\mathrm {T}\)) and bottom (\(\mathrm {B}\)) quark partners in proton–proton collisions at \(\sqrt{s} = 13\,\text {Te}\text {V} \), using data collected by the CMS experiment at the CERN LHC, corresponding to an integrated luminosity of 35.9\(\,\text {fb}^{-1}\)  . The \(\mathrm {T} \overline{\mathrm {T}} \) search is performed by looking for events in which one \(\mathrm {T}\) quark decays via \(\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}\) and the other decays via \(\mathrm {T} \rightarrow \mathrm {b}\mathrm {W}\), \(\mathrm {t}\mathrm {Z}\), \(\mathrm {t}\mathrm {H} \), where \(\mathrm {H}\) refers to the Higgs boson. The \(\mathrm {B} \overline{\mathrm {B}} \) search looks for events in which one \(\mathrm {B}\) quark decays via \(\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}\) and the other via \(\mathrm {B} \rightarrow \mathrm {t}\mathrm {W}\), \(\mathrm {b}\mathrm {Z}\), or \(\mathrm {b}\mathrm {H} \). Events with two oppositely charged electrons or muons, consistent with coming from the decay of a \(\mathrm {Z}\) boson, and jets are investigated, and are categorized according to the numbers of top quark and \(\mathrm {W}\), \(\mathrm {Z}\), and Higgs boson candidates. These categories are individually optimized for \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) event topologies.

The data are in agreement with the standard model background predictions for all the event categories. Upper limits at 95% confidence level on the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) production cross sections are obtained from a simultaneous binned maximum-likelihood fit to the observed distributions for the different event categories, under the assumption of various \(\mathrm {T}\) and \(\mathrm {B}\) quark branching fractions. Comparing these upper limits to the theoretical predictions for the \(\mathrm {T} \overline{\mathrm {T}} \) and \(\mathrm {B} \overline{\mathrm {B}} \) cross sections as a function of the \(\mathrm {T}\) and \(\mathrm {B}\) quark masses, lower limits on the masses at 95% confidence level are determined for different branching fraction scenarios. In the case of a \(\mathrm {T}\) quark decaying exclusively via \(\mathrm {T} \rightarrow \mathrm {t}\mathrm {Z}\), the lower mass limit is 1280\(\,\text {Ge}\text {V}\), while for a \(\mathrm {B}\) quark decaying only via \(\mathrm {B} \rightarrow \mathrm {b}\mathrm {Z}\), it is 1130\(\,\text {Ge}\text {V}\). These lower limits are comparable with those measured by the ATLAS Collaboration [20], also using the \(\mathrm {Z}\) boson dilepton decay channel. The results of the analysis presented in this paper are complementary to previous CMS measurements [21,22,23], and have extended sensitivity in reaching higher mass limits for \(\mathrm {T}\) and \(\mathrm {B}\) quarks.