Search for pair production of heavy vector-like quarks decaying into hadronic final states in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

A search is presented for the pair production of heavy vector-like quarks, $T\bar T$ or $B\bar B$, that decay into final states with jets and no reconstructed leptons. Jets in the final state are classified using a deep neural network as arising from hadronically decaying $W/Z$ bosons, Higgs bosons, top quarks, or background. The analysis uses data from the ATLAS experiment corresponding to 36.1 fb$^{-1}$ of proton-proton collisions with a center-of-mass energy of $\sqrt{s} = 13$ TeV delivered by the Large Hadron Collider in 2015 and 2016. No significant deviation from the Standard Model expectation is observed. Results are interpreted assuming the vector-like quarks decay into a Standard Model boson and a third-generation-quark, $T\rightarrow Wb,Ht,Zt$ or $B\rightarrow Wt,Hb,Zb$, for a variety of branching ratios. At 95% confidence level, the observed (expected) lower limit on the vector-like $B$-quark mass for a weak-isospin doublet ($B, Y$) is 950 (890) GeV, and the lower limits on the masses for the pure decays $B\rightarrow Hb$ and $T\rightarrow Ht$, where these results are strongest, are 1010 (970) GeV and 1010 (1010) GeV, respectively.


Introduction
Many theories beyond the Standard Model (SM) are motivated by the naturalness problem [1], and are intended to resolve the quadratic divergences in the radiative corrections to the Higgs-boson mass. Several extensions to the SM, such as Little Higgs [2,3] and Composite Higgs [4,5] models, have been proposed to address this issue. A common feature of these models is the existence of TeV-scale vector-like quarks (VLQs) that couple preferentially to third-generation SM quarks [6].
VLQs are spin-1/2 colored fermions with left-right symmetric transformation properties under the weakisospin SU(2) gauge group. Unlike chiral quarks, which obtain mass through electroweak symmetry breaking [7][8][9][10][11][12], VLQs can have a gauge invariant mass term mψψ. Therefore, VLQs are not subject to the constraints from Higgs production which highly disfavor additional chiral quarks [13][14][15][16]. VLQs also couple to flavor-changing neutral currents, so a charge1 +2/3 vector-like partner of the top quark, T, could decay2 into W b, Zt or Ht, while a charge −1/3 bottom quark partner, B, could decay into Wt, Z b, or Hb [17]. The branching ratios depend on the VLQ mass and weak-isospin multiplet. Vector-like T and B can occur alone in a singlet scenario. Doublet and triplet scenarios also allow for more exotic X and Y VLQs with charges +5/3 and −4/3, respectively. Charge conservation requires these to decay only via X → Wt and Y → W b. Because this search has not been optimized for X and Y vector-like quarks, they will not be discussed in this paper. 1 Electric charge is measured in units of e. 2 It is assumed that the VLQs decay only into SM particles, and couple to only third-generation quarks.
Many previous searches for pair-produced VLQs by ATLAS and CMS at √ s = 8 TeV [18][19][20][21][22][23] and √ s = 13 TeV [24][25][26][27][28][29][30][31][32] have focused on final states with one or more leptons. Additionally, previous results from CMS at √ s = 8 TeV have included fully hadronic as well as leptonic final states [33,34]. This analysis searches for heavy VLQs produced in pairs and decaying into fully hadronic final states. This channel is complementary to those used in previous ATLAS VLQ searches and is particularly powerful for the B → Hb decay mode, which is difficult to probe with leptonic final states.

ATLAS detector
The ATLAS detector [35] at the Large Hadron Collider (LHC) is centered on the collision point and covers nearly the entire solid angle.3 It consists of an inner tracking detector surrounded by a 2 T superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroid magnets. The inner detector, including the insertable B-layer installed in 2014 [36,37], provides charged-particle tracking information from a pixel and silicon microstrip detector in the pseudorapidity range |η| < 2.5 and a transition radiation tracker covering |η| < 2.0.
The calorimeter system covers the pseudorapidity range |η| < 4.9 and measures the positions and energies of electrons, photons, and charged and neutral hadrons. Within the region |η| < 3.2, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead and liquid-argon sampling calorimeters. The hadronic sampling calorimeter uses either scintillator tiles or liquid argon as active material and steel, copper or tungsten as absorber.
The muon spectrometer comprises separate trigger and high-precision tracking chambers measuring the tracks of muons in a magnetic field generated by superconducting air-core toroid magnets. The precision chamber system covers the region |η| < 2.7 , while the muon trigger system covers the range |η| < 2.4.
A two-level trigger system is used to select which events to save for offline analysis [38]. The first level is implemented in hardware/firmware and uses a subset of the detector information to reduce the event rate from 40 MHz to less than 100 kHz. This is followed by the software-based high-level trigger that reduces the event rate to approximately 1 kHz.

Data and simulated events
The data analyzed correspond to pp collisions with a center-of-mass energy of √ s = 13 TeV recorded by the ATLAS detector in 2015 and 2016. Data quality requirements ensure that all components of the detector were functioning. The full data set corresponds to an integrated luminosity of 36.1 fb −1 .
The primary background for this search is multi-jet events, followed by tt events and minor contributions from single-top-quark and tt + X (X=W, Z, H) events. The multi-jet background is estimated using a data-driven method (Section 5.3), while signal events and other background contributions were simulated via Monte Carlo (MC) generation of LHC collisions that are then passed through a G 4 simulation [39] 3 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the center of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis.
of the ATLAS detector [40]. All simulated events are reconstructed using the same analysis chain as the data. In all MC samples, the top-quark and Higgs-boson masses were set to 172.5 GeV and 125.0 GeV, respectively, and the E G v1.2.0 program [41] was used to simulate the properties of bottom and charm hadron decays.
Simulated events of VLQ pair production, QQ, were produced with the leading-order (LO) generator P v2.2 [42,43] using the NNPDF2.3 LO parton distribution function (PDF) set [44] and passed to P 8.186 [45] for parton showering and fragmentation. The A14 [46] set of tuned parameters is used. VLQs were produced for the isospin singlet scenario with a narrow width and for masses between 700 and 1200 GeV in steps of 50 GeV, with additional events produced at 500, 600, 1300, and 1400 GeV. Additional samples were produced assuming a doublet scenario for VLQ masses of 700, 950 and 1200 GeV, in order to study differences from the different chirality of VLQs arising in singlet and doublet models.
The pair production cross section varies from 3.38 ± 0.25 pb (m Q = 500 GeV) to 3.50 ± 0.43 fb (m Q = 1400 GeV), computed using top++ v2.0 [47] at next-to-next-to-leading order (NNLO) in QCD, including resummation of next-to-next-to-leading logarithmic (NNLL) soft-gluon terms, and using the MSTW 2008 NNLO set of PDFs [48]. Theory uncertainties are estimated by variations of the factorization and renormalization scales and by taking uncertainties of the PDF and strong coupling constant, α S , into account. The latter two represent the largest contribution to the overall theoretical uncertainty in the predicted cross-section and are calculated using the PDF4LHC [49] prescription with the MSTW 2008 68% CL NNLO, CT10 NNLO [50,51] and NNPDF2.3 5f FFN PDF sets.
The tt events were generated using P -B v2 + P 8.210 [52,53] with the CT10 NLO PDF set and the Perugia2012 set of tuned parameters [54] for parton showering. The NLO radiation factor, h damp , was set to 1.5m top . The tt background is split into tt + light-flavor jets (tt + light) and tt + cor b-flavor jets (tt +HF). Single-top-quark production (Wt and t-channel) was generated using P -B v1 + P 6.428 [55][56][57] and the Perugia2012 set of tuned parameters for parton showering and the CT10 NLO PDF set. The tt +V (V=W, Z) and tt +H background was modeled using M G 5_aMC@NLO v2.3.2 [58] as the generator in LO precision with up to two additional partons and in NLO precision, respectively. The parton showering and fragmentation is performed using P 8.210 [45] (P 8.186) for tt +Z and tt + H (tt +W). The contribution from single-top-quark and tt + X events is less than 6% in all signal regions.
Finally, although a data-driven method is used to estimate the multi-jet background, a sample of simulated multi-jet events is used for the training of an algorithm employed to identify boosted objects (Section 4). The simulated multi-jet events were produced with P 8.186 using the A14 set of tuned parameters for the underlying event and the NNPDF2.3 LO PDFs. The renormalization and factorization scales were set to the average transverse momentum (p T ) of the two leading jets.

Object definitions
The main objects used in this search are small-radius (small-R) jets reconstructed from clusters of energy deposited in the calorimeter. A variable-radius re-clustering algorithm [59,60] is then used to find groups of small-R jets that are consistent with the hadronic decays of high-momentum bosons and top quarks. To ensure orthogonality with ATLAS VLQ searches that include leptons [24][25][26][27][28], events containing electrons or muons with p T > 20 GeV are vetoed using the same tight object definitions as in those searches. For a given reconstructed event, the missing transverse momentum, with magnitude E miss T , is calculated from the negative vector sum of the p T of all reconstructed jets, and any reconstructed electrons and muons. A soft energy term is included to account for non-reconstructed particles originating from the hard scatter. It is calculated using only charged tracks matched to the primary vertex to reduce contamination from particles originating from other pp interactions in the same or nearby bunch crossings (pileup) [61].
Small-R jets are reconstructed from calibrated topological energy clusters in the calorimeter using the anti-k t algorithm [62, 63] with a radius parameter of 0.4. They are required to have p T > 25 GeV and |η| < 2.5. Low-p T jets produced in pileup interactions are suppressed using the Jet Vertex Tagger (JVT) algorithm [64]. A jet is removed from the event if it has p T < 60 GeV, |η| < 2.4, and a JVT value lower than 0.59. This requirement on the JVT value has an efficiency of 92% for jets of p T < 60 GeV and |η| < 2.4 originating from the primary vertex. In order to avoid misidentification and overlap of objects, a jet is removed from the event if an electron or muon selected with loosened identification criteria is found within ∆R = 0.2 or if a loosely selected muon is found in a jet that is not well matched to the primary vertex, as in Refs. [24,25].
A small-R jet is b-tagged if it satisfies the 77% working point criterion of the MV2c10 ATLAS b-tagging algorithm [65,66]. Working points are defined by a requirement on the output discriminant and are labeled by the b-jet efficiency they give on an inclusive tt sample. The 77% working point has rejection factors of 6.2 and 134 for jets containing charm hadrons (c-jets) and jets containing light-quark hadrons or gluons (light-jets), respectively. Correction factors are applied to the simulated event samples to correct for differences in the b-tagging efficiencies for b-jets, c-jets, and light-jets between data and simulation. In addition to using b-tagging to select events with the 77% working point, three other working points (60%, 70%, 85%) are used in the context of the boosted-object tagging as described later in this section.
Small-R jets are re-clustered [59] using the anti-k t algorithm with a variable cone size [60] to create variable-radius re-clustered jets (vRC jets). Constituent small-R jets are not allowed to be shared by multiple vRC jets. Because the small-R jets used in the re-clustering are already calibrated, the vRC jets are also calibrated and their uncertainties are obtained directly from the small-R jet uncertainties [67]. A requirement on a p T -dependent variable radius reduces the overlap of boosted objects in the highmultiplicity final state of this search and exploits the fact that the radius separation R between the decay products of a heavy, high-p T particle of mass m can be approximated with R ∼ 2m/p T . The radius parameter threshold is chosen to be R eff = ρ/p T , with ρ = 315 GeV, within the restriction of 0.4 ≤ R eff ≤ 1.2. This results in a good compromise between the accuracy and efficiency of the reconstruction for the objects considered in the final state. To reduce contributions from low-energy pileup, a trimming procedure [68] removes small-R jets from a vRC jet if their p T is less than 5% of the vRC jet p T . The vRC jets are required to have mass greater than 40 GeV, p T > 150 GeV, and |η| < 2.5.
A multi-class deep neural network (DNN) is trained to identify the most likely parent particle of the vRC jets, distinguishing between four categories: V-boson (W-or Z-boson), Higgs-boson, top-quark, and background jets. In simulation the label for a reconstructed signal jet (V-boson, Higgs-boson or top-quark jet) is obtained by matching the vRC jets to a hadronically decaying boson or top quark at generator level within a cone of ∆R = 0.75· ρ/p T . For the Higgs boson, only direct decays into quark pairs are considered. All vRC jets matched to multiple generator-level V bosons, Higgs bosons or top quarks are discarded. The background label is given to any vRC jets reconstructed from simulated multi-jet events. The DNN is trained using the mass, p T , and number of constituent jets of the vRC jet, as well as the four-momentum vectors and b-tagging information of the three highest-p T constituent small-R jets as input. and rectified-linear-unit activation functions, whereas the output layer uses a sigmoid function. As the performance of the DNN tagger is dependent on the architecture and training hyperparameters, DNNs with different number of layers, learning rate, L1 regularizer and batch size are tested to define the architecture and training hyperparameters. The p T distribution of the background is reweighted to match the signal distribution. In this way the DNN is prevented from learning the differences between the p T distributions of signal and background jets, while allowing for learning relations between the p T of the vRC jets and other input features. By using only properties of the calibrated small-R jets as input to the tagger, all jet-related systematic uncertainties can be propagated through the DNN by varying the corresponding properties of the small-R jets. To reduce the four-dimensional DNN output information (D DNN ), the outputs of the different classes are combined by building a discriminant function.
The discriminant function P for a V boson, Higgs boson and top quark is given by respectively. The relative weighting factors of 0.9 for background jets and 0.05 for V-boson, Higgs-boson or top-quark jets are chosen as a compromise between background rejection and the ability to discriminate amongst signal sources. For each signal discriminant P, an optimized working point is defined to obtain a boosted-object tagger with a specified signal efficiency. The discriminant functions and the corresponding thresholds for these working points are shown in Figure 1, where |η|, p T , and m refer to the pseudorapidity, transverse momentum and mass of the vRC jet. The multi-peak behavior in some of the discriminant functions is caused by the variations in important vRC jet features used as input to the DNN, such as the mass and number of constituent small-R jets variations depending on the p T range of the vRC jet or the number of b-tagged and non-b-tagged constituent small-R jets, that can cause one signal type to mimic another." The Vand Higgs-boson taggers use 70% working points, which correspond to the thresholds P(V) > −0.2 and P(H) > 0. 35. The top-quark tagger operates at a 60% working point using a threshold of P(t) > 0.1. The resulting signal efficiency and background rejection (estimated from simulated multi-jet events) for each boosted-object tagger is shown as a function of p T in Figure 2.
To handle the ambiguities due to multiple-tagged vRC jets, additional discriminant functions, shown in Figure 3, are defined. Optimized thresholds, shown in each sub-figure, are chosen to resolve doubletagged vRC jets. Higgs bosons are more frequently triple-tagged than V bosons or top quarks, so triple-tagged vRC jets are tagged as a Higgs boson. The vRC jets that are tagged as a V boson, top quark, or Higgs boson are referred to as V-tagged, top-tagged, and Higgs-tagged, respectively.
The shape of the vRC jet mass distribution before and after the final boosted-object tagging is shown in Figure 4 for each jet type. As expected, each tagger preferentially selects vRC jets with a mass near   the mass of the desired particle. For the top quarks, vRC jets with a mass near the W-boson mass are generally V-tagged (dominant at low p T ) and Higgs-tagged.

Analysis strategy
The search presented in this paper focuses on all-hadronic final states with small E miss T , which allows it to be sensitive to all possible final states involving hadronic decays of W, Z, and Higgs bosons and top quarks. The key aspect of this search is to suppress multi-jet background and accurately model multi-jet events that satisfy the selection criteria. As a first step, the multi-jet background is reduced by requiring multiple high-p T and b-tagged small-R jets. As a second step, events are rejected if they do not contain vRC jets that originate from either a V boson, Higgs boson, or top quark as identified using the DNN boosted-object tagger. Events are then categorized according to the numbers of V-tagged, Higgs-tagged and top-tagged vRC jets and of b-tagged small-R jets and are divided into twelve non-overlapping signal regions, in order to be sensitive to all possible VLQ decays. Finally, multi-jet events are distinguished from signal events by calculating, for each signal region, a signal probability using the matrix element method [75]. This signal probability is then used in a binned profile-likelihood fit in order to extract the signal strength and improve the background modeling. The multi-jet background is estimated in each signal region using a bin-by-bin 'ABCD' method, which is described in Section 5.3. The analysis strategy is optimized while assuming pair production of VLQs and considering all possible fully hadronic decay modes.

Event selection and classification
Data were collected using a combined trigger that requires a single jet with p T > 100 GeV at the first trigger level and a total scalar sum of the transverse momenta of all track particles and energy deposits H T > 1000 GeV at the high-level trigger. An offline threshold of H T > 1250 GeV ensures that this trigger is fully efficient. Events are required to have exactly zero leptons and E miss T < 200 GeV to remove background and maximize the significance of the signal. Events enter the signal regions if they contain at least four selected small-R jets with descending p T thresholds of 300, 200, 125, and 75 GeV and at least two small-R jets that are b-tagged, where individual jets can satisfy one or both criteria. In addition, the events must have at least two vRC jets tagged as a V or Higgs boson and satisfy E miss T > 40 GeV. The E miss T requirement rejects significantly more background than signal. For example, the E miss T requirement is 71-82% efficient for the various decay modes of a signal with a mass of m V LQ = 1 TeV, but only 55% efficient for simulated multi-jet background events. Sources of E miss T in VLQ pair production can include true sources, such as Z → νν decays or leptonic decays of W bosons and top quarks with a soft or misreconstructed lepton, as well as E miss T from mismeasurement of high-energy jets.
The events are then classified into twelve different signal regions based on the number of Vand Higgs-tags (VV, V H, HH), top-tags (0, 1, ≥ 2), and b-tags (2, ≥ 3), as shown in Table 1. The regions are designed to cover all of the possible VLQ decays and enhance the ratio of signal events to SM background events. Figure 5 shows the fraction of events from each background source that contribute to each signal region after the full event selection and the background-only fit to data described in Section 7. In addition to the signal regions, nine validation regions are also defined in order to validate the multi-jet background estimation and evaluate a closure uncertainty for the method. The two regions that are used to validate the multi-jet background estimation are defined to have exactly two b-tagged jets, two Higgs-tags, and no top-tags. The seven regions that are used to evaluate the closure uncertainty require exactly one b-tagged jet and the same number of V-Higgs-and top-tags as in each of the signal regions. Table 1: Summary of the definition of the twelve signal regions in the analysis. The number of b-tags is based on all small-R jets, including those used to construct vRC jets with V-, Higgs-, or top-tags. The last two signal regions require two bosons of any type X (V or Higgs boson). The rightmost column lists the matrix element method (MEM) final states used to define the signal hypothesis in Section 5.2.
Region Name V-tags H-tags top-tags b-tags MEM final states  Figure 5: Relative size of all background contributions in each signal region after the background-only fit to data described in Section 7. 'Others' refers to backgrounds from single-top-quark and tt + X production.

Matrix element method
The matrix element method [75] has been utilized for measurements [76-78] and searches for SM physics processes [79][80][81][82][83][84][85]. This analysis applies the method to a search for physics beyond the SM. This method requires the calculation of an event-based probability density function P i (x|α) for a given physics process i described by the matrix element of the process and a set of theoretical and experimental parameters α: The numerical integration is performed over the phase space of the initial-and final-state particles and can be time consuming. In this equation, x and y represent the four-momentum vectors of all initial-and final-state particles at reconstruction and parton level, respectively. The Lorentz-invariant flux factor F 5 and phase-space element dΦ N describe the kinematics of the process. The transition matrix element M i is defined by the Feynman diagrams of the hard-scattering process. The functions f (p A ) and f (p B ) are the PDFs for the initial-state partons with momenta p A and p B . The transfer functions W (y| x) map the detector quantities x to the parton-level quantities y. Finally, the effective cross section σ eff i normalizes P i to unity taking acceptance and efficiency into account.
The reconstructed objects in an event can be combined to form multiple candidate VLQ final states. The process probability density is calculated for each allowed assignment permutation of the jets to the finalstate quarks and bosons. A process likelihood function is then built by summing the process probabilities of each allowed assignment permutation. The vRC jets are assigned to final states according to their DNN VLQ boosted-object-tag label (V-tagged, Higgs-tagged, or top-tagged) and are permuted if they have the same label. If more than two vRC jets are tagged as a boson (V-tagged or Higgs-tagged), only the two with the highest transverse momenta are used. If b-quarks occur in the hypothesized final state, up to five different b-tagged small-R jets are assigned to the final state and freely permuted. These b-tagged jets are allowed to overlap with the vRC jets and could have been used already in the reconstruction of the vRC jet.
The transition matrix element defines the hypothesis being tested and is calculated using M G 5 in LO precision. The VLQ pair-production matrix element calculation is performed using the Feynrules [86] model as defined in Ref. [87]. In this analysis, only probabilities of signal hypotheses are calculated, since the dominant background is from multi-jet processes, for which it is difficult to define a model in the matrix element method. The second most important tt +jets background was studied as a background hypothesis, but its inclusion does not improve the sensitivity of the search. Top quarks, V bosons, and Higgs bosons are assumed to be reconstructed as vRC jets and are hence not decayed in the matrix element calculation.
In each signal region, the signal hypothesis is computed from all Feynman diagrams of vector-like T or B pair production resulting in the same number of top quarks, V bosons, and Higgs bosons as defined in Table 1. Following the definition of the signal regions (XX,2t,2b) and (XX,2t,3b), all Feynman diagrams resulting in final states with two top quarks are used and no distinction is made based on the number of V and Higgs bosons. Combining these diagrams into a single hypothesis increases the performance significantly and allows mis-tags of the V and Higgs bosons. Because there is no direct decay of VLQs into a final state with two Higgs bosons and one top quark, the same diagrams as used for the (VH,1t,2b) and (VH,1t,3b) signal regions are used in the (HH,1t,3b) region taking mis-tags into account. Preliminary studies indicated that this analysis would be sensitive to VLQ masses around 900 GeV, therefore, in the calculation of the matrix elements, the masses of the vector-like B and T quarks are set to 900 GeV. The analysis sensitivity becomes slightly degraded when considering signal samples with a significantly higher VLQ mass.
The transfer functions are parameterized as single-Gaussian functions, which is a good compromise between separation power and reasonable integration time. For the modeling of the parton distribution functions, the CTEQ6L1 set from the LHAPDF package [88] is used. The integration is performed using VEGAS [89]. Due to the complexity and high dimensionality, adaptive MC techniques [90], simplifications and approximations are needed in order to perform the integration in a reasonable time. The matrix element calculation is accelerated by evaluating only the most significant helicity states, which are identified at the beginning of each integration. The dimensionality of the integration is reduced by assuming that the final-state object directions in η and φ are measured with negligible uncertainty. The total momentum conservation and the negligible transverse momentum of the initial-state partons allow further reduction. No change of integration variables is performed in order to allow a general treatment of all signal regions. The integration variables are the energies of the top quarks, b-quarks, V, and Higgs bosons according to their numbers as defined for each region. The total integration volume is restricted by requiring the difference between the parton-level quantities and the observed values to be within five standard deviations of the width of the transfer functions. Finally, the likelihood contributions of all allowed assignment permutations are coarsely integrated and sorted by their contribution, then the full integration is performed with a decreasing precision. The logarithm of the resulting signal likelihoods (signal LLH) is used in each signal region as the final discriminating variable. Normalized distributions of the signal LLH for the total background and signal simulations in the most sensitive signal regions assuming exclusive T → W b, B → Z b, and B → Hb decays are shown as examples in Figure 6. The binning in the signal LLH distribution is the same as that shown for the corresponding regions in Figures 10  and 12. The separation given in Figure 6 between signal and background is defined by the formula where S(x) and B(x) are the signal and background yields per bin and S and B are normalized to unity.

Background estimation
The dominant multi-jet background is estimated using a data-driven double sideband method ('ABCD'). The four regions are orthogonal and there is no significant correlation between boson tagging and E miss T . The level of correlation is evaluated by checking the correlation factor between the two variables in simulated multi-jet events, which is found to be consistent with zero.
In the control regions A, B, and C, the non-multi-jet contributions are subtracted from the data using simulation. The relationship between the yields, N, in the signal region, D, and the control regions is given by N D = N C × (N B /N A ). This simple scaling is performed on a bin-by-bin basis in the signal LLH distribution to produce the expected multi-jet shape and normalization in the signal region. This procedure is followed separately for each of the twelve signal regions. Seven validation regions are also defined, with the same V-, Higgs-, and top-tagging requirements as the signal regions, but with exactly one b-tagged jet. These regions are used to evaluate a closure uncertainty, described in Section 6. Two examples of these validation regions can be seen in Figure 7, where the only uncertainties taken into account are those from statistical sources and related to the detector simulation.  The binning that is used for each region is determined by the number of events in the A, B, and C control regions. It is required that there are a sufficient number of events in each bin of the control regions (at least 50) to produce a sufficiently smooth distribution.
To evaluate the performance of the background estimation method with all uncertainties, two regions kinematically close to the signal regions, but with very low expected signal contribution, are also defined. These regions have two Higgs-tagged vRC jets, exactly two b-tagged small-R jets, and either zero or one top-tagged vRC jet. Good agreement is observed in these regions, as shown in Figure 8.
Standard Model backgrounds from tt, single-top-quark, and tt + X processes are estimated with simulated events, described in Section 3. The normalization and shape are taken directly from simulation for all of these processes.

Systematic uncertainties
The systematic uncertainties considered in this analysis arise from uncertainties in the treatment of the luminosity, object reconstruction and background modeling. Each source of uncertainty is treated as a nuisance parameter in the final likelihood fit, as described in Section 7. Different sources of uncertainty are assumed to be uncorrelated; however, a given uncertainty is assumed to be 100% correlated across all regions and samples. For each source of systematic uncertainty, the effect on the analysis is evaluated by propagating a ±1σ variation of the quantity in question.

Luminosity and pileup
The uncertainty in the integrated luminosity of the 2015 and 2016 data set is 2.1%. It is derived, following a methodology similar to that detailed in Ref.
[91], from a calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016. Because MC events are simulated with different pileup conditions than observed in data, the events are corrected to have the same pileup distributions as the data and an uncertainty is assigned to account for the uncertainty in the ratio of the predicted and measured inelastic proton-proton cross section [92].

Reconstructed objects
Several systematic uncertainties in the simulated background and the signal predictions arise from the reconstruction and identification of the selected reconstructed objects, as described in Section 4, due to the determination of correction factors applied to compensate for differences between data and predictions. The most important sources in this category are the uncertainties associated with jets, missing transverse momentum, and flavor tagging. Other sources, such as lepton reconstruction (affecting the lepton veto), are also considered, but have a negligible impact on the results. The impact on both shape and normalization is taken into account for the following uncertainties.
Jets In case of the small-R jet selection, uncertainties arise from the jet reconstruction, the jet energy and mass scale calibrations, the JVT requirement, and corrections to the jet energy and mass resolutions. The most significant uncertainties associated with small-R jets are from energy scale and energy resolution. The energy scale is determined using the transverse momentum balance between a jet and a reference object such as a photon, Z boson, or another jet [93]. The uncertainty in the energy scale ranges from less than 1% to around 5% for |η| < 0.8 and p T up to 500 GeV. Jets with higher |η| have an additional uncertainty of up to 2%. The jet energy resolution is measured by studying dijet events in data and simulation [94]. The jet energy resolution in data and simulation are found to agree within 10% and the differences are used to determine the relative systematic uncertainties, which range from 10% to 20%. Additional uncertainties are considered for the jet mass scale and mass resolution, but are found to have little impact on the search sensitivity. The uncertainties associated with vRC jets are inherited from the small-R jet uncertainties.

Missing transverse momentum
The E miss T is sensitive to changes in the momenta of the reconstructed objects, namely the small-R jets, as well as the additional soft term that accounts for low-energy deposits not associated with a reconstructed object. Uncertainties from the reconstructed objects are already accounted for. A soft-term uncertainty is assigned to account for variations in the modeling of the underlying event that change the amount of unclustered energy. The uncertainties in the yields are in the range 0.0-18.7% for simulated samples and 0.0-8.2% for the multi-jet background.
Flavor tagging Uncertainties in the correction factors for the b-tagging identification response are obtained by comparing the simulated event samples with dedicated flavor-enriched samples in data [65]. An additional term is included to extrapolate the measured uncertainties to the high-p T region of interest. This term is calculated from simulated events by considering variations of the quantities affecting the b-tagging performance such as the impact parameter resolution, percentage of poorly measured tracks, description of the detector material, and track multiplicity per jet. The dominant effect on the uncertainty when extrapolating to high p T is related to the different tagging efficiency when smearing the track impact parameters based on the resolution measured in data and simulation.
Most of the vRC jet-tagger flavor-tagging uncertainties can be derived by propagating the small-R jet uncertainties through the DNN. An additional uncertainty associated with b-tagging is evaluated to take into account the use of b-tagging information in the vRC jet-tagger. This is a p T -dependent uncertainty in the vRC jet-tagging efficiency, considered separately for V-boson, Higgs-boson, and top-quark tagging. This uncertainty in the yields ranges from 4.0 to 11.9% for simulated samples and from 0.3 to 9.4% for the multi-jet background.

Background modeling
A theory cross-section uncertainty of 5.3% is taken for the combined small backgrounds, which are dominated by single-top-quark processes [95].

Multi-jet estimation
The dominant multi-jet background is estimated using a data-driven ABCD technique, as described in Section 5.3. To quantify a closure uncertainty for this method, the difference between the prediction and data in the one-b-tag validation regions is propagated as an overall normalization uncertainty to the corresponding two-and three-b-tag signal regions. The impact of including shape information in this uncertainty is negligible. To allow potential differences in performance as a function of jet multiplicity, the uncertainties are taken to be uncorrelated between the regions with exactly two b-tags and at least three b-tags.
Another uncertainty is taken from the impact on the multi-jet prediction of potential signal contamination in the validation regions. Detector-related uncertainties associated with all backgrounds estimated using simulation, as well as modeling uncertainties in tt processes, are also propagated through the multi-jet estimation via subtraction of non-multi-jet events in the validation regions. These uncertainties take into account differences in both shape and normalization.
In addition to the systematic uncertainties, each bin of the multi-jet prediction is assigned an uncertainty to account for statistical uncertainties in the CRs propagated through the ABCD method. Along with the statistical uncertainty of the data in the SR, these tend to have the largest impact on the sensitivity of the analysis.
tt modeling For the tt background, systematic uncertainties are considered for variations in initialand final-state radiation (ISR/FSR), choice of parton shower, and choice of matrix-element generator. Each of these sources of uncertainty are considered as separate nuisance parameters in the likelihood fit. These are evaluated using alternative simulated tt samples. The uncertainty in the treatment of radiative effects is estimated by varying the NLO radiation factor h damp and the factorization and renormalization scales in a correlated way to produce more or less radiation. Alternative samples produced with M G 5_aMC@NLO v2.3.3 + P 8.212 and P -B v2 + H 7.0.1 [96] are used to evaluate generator and shower model uncertainties, respectively. Due to the limited number of events in the alternative tt samples, the uncertainties are taken into account after merging signal regions with two and at least three b-tagged jets. These uncertainties are in the range 1.4-33% (13-51%) for the normalization of tt + light (tt + HF). Because the predicted cross sections of tt + light and tt + HF are not well known for the phase space of the signal regions, separate normalization factors are assigned to each of these two contributions and are allowed to float freely in the profile likelihood fit.

Statistical analysis
The statistical analysis quantifies the probability of compatibility between the measured data, expected SM background, and expected signal. The signal LLH distributions for the twelve signal regions are tested simultaneously for the presence of a VLQ signal. Hypothesis testing is performed using a modified frequentist method based on a profile likelihood, taking into account systematic uncertainties as nuisance parameters. The statistical analysis is based on a binned likelihood function L(µ, θ) constructed as the product of Poisson probability terms over all bins. In this function, µ is a multiplicative factor applied to the predicted production cross section times branching ratio of the signal, and θ is the set of nuisance parameters, implemented in the likelihood function as Gaussian or log-normal priors. In addition, there are two unconstrained parameters in the fit, corresponding to the total normalization of tt + light and tt + HF.
The test statistic q µ is defined as the profile likelihood ratio: q µ = −2ln(L(µ,θ µ )/L(μ,θ)),whereμ andθ are the values of the parameters that maximize the likelihood function (with the constraint 0≤μ ≤ µ), and θ µ are the values of the nuisance parameters that maximize the likelihood function for a given value of µ.
Upper limits on the signal production cross section for each of the signal scenarios considered are derived by using q µ in the CL s method [97,98], where CL s is computed using the asymptotic approximation [99]. For a given signal scenario, values of the production cross section that yield CL s < 0.05 are excluded at ≥95% confidence level (CL).

Results
Following the prescription described in Section 7, the profile likelihood fit for the background-only hypothesis is performed simultaneously in all signal regions. The post-fit event yields are given in Table 2 and Figure 9 shows a comparison between the predicted and observed numbers of events in all signal regions both before and after the fit.
The most notable shift in the post-fit yields is that of the tt + HF normalization. The overall change in normalization is by a factor slightly greater than two, which is achieved in the fit through a shift of the tt + HF normalization factor, as well as through pulls of systematic uncertainties, such as tt modeling and jet energy resolution uncertainties. The post-fit distributions of the signal LLH from each signal region are shown in  No significant excess of signal-like events is observed, and the analysis proceeds to set upper limits on the production cross section of TT and BB events in various scenarios. The sensitivity is mainly limited by the statistical uncertainty in the signal regions and in the control regions for the ABCD method. For example, if only statistical uncertainties and normalization factors are taken into account, the expected (observed) cross-section limit for BB → HbHb with m B = 1 TeV only changes by 5% (11%).
In a given scenario, a lower limit on the VLQ mass can be obtained by comparing the cross-section limits with the predicted cross section as a function of mass [47]. Figure 13 shows the expected and observed upper limits on the TT and BB cross section at 95% CL as a function of the VLQ mass in the scenario where the VLQ decays purely via the Higgs decay mode (TT → HtHt or BB → HbHb), as well as in the benchmark scenario of the (B, Y ) doublet. In this scenario, a B VLQ will decay almost equally into Z b and Hb, although the exact branching ratios depend on mass. For example, for m B = 1 TeV, B(B → Z b) = 0.51 and B(B → Hb) = 0.49 [42]. Only contributions from the B VLQ are considered, so the limit is conservative. In the case of a (B, Y ) doublet, B masses below 950 GeV are excluded at 95% CL.
To evaluate the level of sensitivity of the results to the weak-isospin of the VLQ, samples of VLQ events with masses of 700, 950, and 1200 GeV were generated for an SU(2) doublet T and B quark and compared with the SU(2) singlet samples. Small differences between the limits are observed in decay modes with bottom quarks in the final state, where the SU(2) singlet produces a slightly weaker limit due to the slightly Others  lower average momenta of decay products in the singlet final state. Thus, limits on the (B, Y ) SU(2) doublet, which are taken from scaling the SU(2) singlet samples to doublet branching ratios, represent a slightly conservative limit. In final states with top quarks, the SU(2) singlet produces a slightly stronger limit, due to a slightly higher efficiency in VLQ DNN top-tagging. Therefore, a limit on SU(2) doublets with T VLQs are not included here.
The largest difference between the observed and expected limits is for B(B → Hb) = 1 with a VLQ mass around 950 GeV. This results from a deficit in data in the final two bins of the (HH,0t,3b) signal region and the fact that the matrix element calculation for final states with two bottom quarks has its maximum sensitivity for masses near 900 GeV.
By reweighting the relative fractions of the three T (B) decay modes, it is possible to test all combinations of branching ratios. Figure 14 shows  Table 3.       Contour lines, alternating solid and dashed lines, are provided to show sensitivity to different VLQ masses across the planes. Signal hypotheses are considered in a mass range of 500-1400 GeV, so the white space on the observed limit figures corresponds to branching ratios where there is no observed exclusion above a mass of 500 GeV.

Conclusion
A search for pair production of vector-like quarks in the all-hadronic final state is presented using 36.1 fb −1 of collision data collected by the ATLAS detector at the LHC in 2015 and 2016. The analysis selects events with high-p T small-R jets and multiple b-tags. Small-R jets are combined using a variable-R clustering algorithm and then classified with a neural network as a V-boson, Higgs-boson, top-quark, or background jet. A signal log-likelihood calculated via the matrix element method is used as the final discriminant across multiple categories based on the number of V/H-tags, top-tags, and b-tags. The analysis targets all third-generation decays of VLQs, but it is particularly powerful for the B → Hb decay mode, which is difficult to probe with leptonic final states. The observed data are consistent with expected background and a 95% CL limit is placed on VLQ pair production as a function of the hypothetical VLQ mass.