Signals with six bottom quarks for charged and neutral Higgs bosons

In extensions of two Higgs doublet models with vectorlike quarks, the decays of vectorlike quarks may easily be dominated by cascade decays through charged or neutral Higgs bosons leading to signatures with 6 top or bottom quarks. Since top quark decays also contain bottom quarks, a common signature for many possible decay chains is 6 bottom quarks in the final state. We present a search strategy focusing on this final state and find the mass ranges of vectorlike quarks and Higgs bosons that can be explored at the Large Hadron Collider. Among other results, the sensitivity to the charged and neutral Higgs bosons, extending to about 2 TeV, stands out when compared to models without vectorlike matter.


Introduction
Among the simplest extensions of the standard model (SM) is the two Higgs doublet model (2HDM) featuring two neutral Higgs bosons, the CP-even, H, and the CP-odd, A, (if CP is conserved) and a pair of charged Higgs bosons, H ± , in addition to the neutral CP-even Higgs boson, h, which is experimentally constrained to have properties very close to the SM Higgs boson. However, the typically dominant decay modes of the additional Higgs bosons, tt, bb, tb andtb, in the type-II 2HDM are very challenging at the Large Hadron Collider (LHC) due to huge QCD backgrounds. For example, the limits on H + → tb currently exist only for small and large tan β (the ratio of vacuum expectation values of two Higgs doublets, that set the Yukawa couplings of Higgs bosons) and the LHC is not yet sensitive to this decay in most of the range, 2 tan β 35 [1,2]. Although the limits will certainly improve with larger data sets, even with 3 ab −1 of integrated luminosity it is not expected that the LHC would constrain the charged Higgs boson significantly in this range. 1 In extensions of two Higgs doublet models with vectorlike quarks, the decays of vectorlike quarks may easily be dominated by cascade decays through charged or neutral Higgs bosons leading to signatures with 6 top or bottom quarks [8], see figure 1. The heavy Higgs bosons are effectively pair produced with QCD size cross sections and lead to final states with very small irreducible SM background (the dominant background happens to originate from QCD multi-jet final states). Thus searching for them in top-and bottom-rich events presents a unique opportunity for the LHC. Depending on model parameters and especially JHEP07(2020)241 the hierarchies in the masses of new quarks and Higgs bosons, many decay chains are possible, and several of them can be simultaneously sizable. Optimal searches for individual possibilities could be designed but the large number of them makes this approach impractical. However, since top quark decays also contain bottom quarks, a common signature for many possible decay chains is 6 bottom quarks in the final state.
We present a search strategy focusing on the 6b final state and find the mass ranges of vectorlike quarks and Higgs bosons that can be explored at the LHC. Although the strategy is tailored for the b 4 → bH → bbb process we find it to be very effective for all the processes in figure 1. Among other results, the sensitivity to the charged and neutral Higgs bosons, extending to about 2 TeV, stands out when compared to models without vectorlike matter.
Complementary signatures of cascade decays of heavy neutral Higgs bosons through vectorlike and SM quarks, relevant when vectorlike quarks are lighter than the new Higgs bosons, were recently studied in ref. [9]. Besides the opposite required hierarchy in masses, the main difference from the decay modes discussed here is that the flavor violating decays H → t 4 t and b 4 b, resulting from Yukawa couplings that mix vectorlike and SM quarks, have to compete with the usual decay modes of heavy Higgs bosons, H →tt andbb. Therefore the branching ratios of heavy Higgses and thus the sensitivity to these channels highly depends on the size and the structure of Yukawa couplings of vectorlike quarks. On the other hand, the decays of the lightest vectorlike quark are necessarily flavor violating. Thus, if the decay channels through heavy Higgs bosons are kinematically open, these can dominate even for very small values of all Yukawa couplings. Large branching ratios for t 4 → Ht, H ± b and b 4 → Hb, H ± t, even close to 100%, are abundant in random scans of these couplings [8].
Similar signatures of heavy Higgses and vectorlike leptons were studied in refs. [6,[10][11][12][13]. 2 A simple embedding into grand unified theories suggests models with both vectorlike quarks and leptons. The supersymmetric extension with a complete vectorlike family also provides a possibility to understand the values of all large couplings in the SM from the IR fixed point structure of the renormalization group equations [17][18][19]. This possibility, that requires both vectorlike quarks and charged and neutral Higgs bosons near the TeV scale, possibly within the reach of the LHC, is the main motivation for the study presented in this paper.
The same or very similar signatures can also be found in composite Higgs models or models with various top partners. For recent analyses, see for example refs. [20][21][22][23]. A detailed study of the 6t final state, motivated by composite Higgs models, that focuses on the presence of up to three high transverse momentum leptons has been presented in ref. [24]. This signature is identical to one of the signatures of vectorlike quarks in 2HDM, see the bottom left process in figure 1. In addition, these signatures can arise in models with Z or W with the neural (charged) Higgs boson in figure 1 replaced by Z (W ), e.g. the Z with couplings to b 4 and b was suggested in connection with the Z-pole anomalies [25,26] or tensions in rare B decays [27,28]. However, the rates for the final states in other models JHEP07(2020)241 might not reach the rates that are possible in the 2HDM. Finally, we mention that the 6b final state has been studied in ref. [29] in the context of SM triple Higgs production at future colliders which has kinematics very different from the signatures considered in this work.
This paper is organized as follows. In section 2 we discuss the decay modes in figure 1 focusing on the branching ratios to various final states expected in the model. In section 3 we suggest search strategies for final states with six bottom quarks and estimate the reach of the LHC for heavy Higgs bosons and vectorlike quarks. We summarize the main results in section 4. Supplementary details of the analysis are provided in appendices.

Decays of vectorlike quarks through heavy Higgs bosons
Decays of vectorlike quarks through heavy Higgs bosons are especially well motivated, since in 2HDMs they can dominate for a generic set of even very small Yukawa couplings that mix vectorlike and SM quarks. This has been studied in detail in ref. [8], where a complete extension of the type-II 2HDM by vectorlike quarks with the same quantum numbers as SM quarks (and corresponding conjugate states) was considered. A description of the model, all the relevant formulas for couplings of Higgs and gauge bosons together with a detailed discussion of cascade decays of vectorlike quarks (assuming they dominantly mix with the third generation of SM quarks) can be found there. Here we summarize the findings at the level sufficient for the motivation of signatures and basic understanding of the results. In addition, we provide new results for scenarios where decay chains through all heavy Higgs bosons are simultaneously kinematically open.
In figure 1 we show representative cascade decays from pair production of a new bottom-like quark, b 4 , and top-like quark, t 4 , through charged and neutral Higgs bosons. These include decays of vectorlike quarks through two neutral Higgs bosons, two charged Higgs bosons, and cases when one branch of the decay proceeds through a neutral Higgs and the other through a charged Higgs (only one possibility for each quark flavor is shown). There are also processes where H in any branch is replaced with A. In addition, there are -3 -

JHEP07(2020)241
processes with any other allowed decay modes of the individual Higgs bosons. However, the decay chains shown are those that can have combined branching ratio close to 100% in the type-II 2HDM. For example, we do not show and will not consider cases such as t 4 → Ht → tbb or b 4 → Hb → ttb. In these decay chains, the branching ratios of vectorlike quarks and those of the neutral Higgs boson are large in different regions of tan β [8,9] and the combined branching ratio for the whole chain is not large for any tan β. Thus these topologies are suppressed compared to the ones shown. Note however that in other models such decay chains might be significant.
In order to understand why the decays of vectorlike quarks through heavy Higgs bosons can easily dominate, it is crucial to note that the decays of the lightest new quark are necessarily flavor violating. They can decay into a SM quark and W , Z, h, H, A, or H ± . The relevant couplings originate from Yukawa couplings that mix vectorlike (SU(2) doublet or singlet) quarks and SM quarks. The leading dependence of the t 4 and b 4 partial decay widths on the Yukawa couplings and tan β is summarized in table 2 of ref. [8] (the leading dependences of partial decay widths to A and H are identical). The important observation is that partial decay widths to heavy Higgs bosons are controlled by either a different Yukawa coupling than the partial decay widths to W , Z, h, or they have opposite tan β dependence. Thus, if the decay channels through heavy Higgs bosons are kinematically open, these can dominate even for very small values of all Yukawa couplings. Large branching ratios for t 4 → Ht, H ± b and b 4 → Hb, H ± t, even close to 100%, are abundant in random scans of these couplings [8].
The analysis in ref. [8] assumed that decays to only one of the heavy Higgs bosons were kinematically open. Since some of the topologies in figure 1 involve decays to two different Higgs bosons, we extend the analysis to the case when decays through all three heavy Higgs bosons are kinematically open. Such possibilities are expected to occur in supersymmetric extensions where heavy Higgs bosons are typically almost degenerate. In figure 2 we plot the allowed branching ratios of t 4 and b 4 into heavy charged and neutral Higgses vs. tan β resulting from a random scan over the parameter space described in ref. [8] with the additional assumption that m H = m A = m H ± = 1 TeV (the actual mass is not important as long as the decays remain kinematically open). In addition, in figure 3 we plot the generated scenarios in planes of various branching ratios. For a vast majority of the points the branching ratio to A is almost identical to the branching ratio to H, as is indicated in the top plots in figure 3 and thus we do not plot its tan β dependence. 3 The color coding of the points allows for easy understanding of numerical results and it is the same as in ref. [8].
In all cases, the behavior of the branching ratios and their dependence on tan β is straightforward to understand from their simple relationships to the couplings summarized in table 2 of ref. [8]. We find that, when decay channels through all three heavy Higgses are kinematically open and dominate, typical ratios of partial widths of t 4 to H, A and JHEP07(2020)241 In all panels, orange (green) points correspond to a parameter scan assuming couplings to H u (H d ) only; gray points correspond to a parameter scan when all couplings are allowed. In all figures, decay modes through any heavy Higgs boson are allowed and m H = m A = m H ± = 1 TeV is assumed. Other details of the parameter space scan are the same as in ref. [8].
In summary, we find that the branching ratios of t 4 → H ± b → tbb and b 4 → H(A)b → bbb can easily be close to 100% for any medium to large tan β even if all the decay modes are kinematically open. Thus these decay modes are especially well motivated. In the following section we focus on the 6b final state and estimate the sensitivity of the LHC to b 4 and H(A). We also use the same analysis to find sensitivities for other decays including

Search strategies and reach at the LHC
The final states shown in figure 1 are 6b (b 4 → H cascade), 6t (t 4 → H cascade), 4b2t and 4t2b (both requiring the presence of H ± ). In this section we develop an analysis strategy, based on the presence of four or five b-jets with high transverse momentum and large total transverse energy, which applies to all these final states. Obviously this strategy is expected -5 -JHEP07(2020)241 Figure 3. The allowed branching ratios of t 4 and b 4 into heavy charged and neutral Higgses in the general scenario with all couplings allowed. The new quark is 95% or more singlet-like (red), 50%-95% singlet-like (purple), 50%-95% doublet-like (cyan) or 95% or more doublet-like (blue). In all figures, decay modes through any heavy Higgs boson are allowed and m H = m A = m H ± = 1 TeV is assumed. Other details of the parameter space scan are the same as in ref. [8].
to be the most effective for the 6b final state, with slightly weaker sensitivity as the number of top-jets increases.
In figure 4 we show the LO production cross section of pair produced vectorlike quarks at the LHC as a function of their masses m t 4 ,b 4 . We show curves for 14 TeV and 13 TeV center of mass energies in green and blue, respectively.
To generate each signal, we implement the model discussed in ref. [8] into FeynRules [30] to produce a UFO file. Parton level events for signal and background are then generated with MadGraph5 [31] and subsequently showered and hadronized with Pythia8 [32,33]. Finally we used Delphes [34] (with standard settings) to simulate detector effects. The last step is imperative to our analyses as a crucial aspect of the search strategy relies on the estimation of the efficiency to observe multiple b-jets with high transverse momentum. σ (pp A useful kinematical quantity is the total transverse energy of reconstructed b-tagged jets: which for vectorlike quark masses above 1 TeV can easily exceed 2-3 TeV. Cutting on H T b allows us to strongly reduce multi-jet QCD backgrounds while, at the same time, preserving sizable signal efficiency. The dominant obstacle to an all-hadronic analysis is the estimation of the QCD background. In particular, besides irreducible backgrounds with at least four or five b quarks at the parton level (e.g. pp → 4b), there are multi-jet final states in which regular jets are mistagged as b (e.g. pp → 2b2j). The mistag rate receives a small contribution (1 − 2%) from jets which do not contain b-hadrons after hadronization but are nevertheless tagged by the algorithm, and a much larger one (we found about 7% from a study based on simulations with Pythia8) from jets in which a gluon splits into a bb pair. The parton level cross sections for processes in which a pair of b quarks is replaced with two jets receive an enormous combinatoric enhancement and, after detector simulation, dominate the background to a multi b-jets signal. The problem of separating genuine b-jets (which we refer to as "1b") from jets in which a gluon splits into bb (which we refer as "2b") has been already investigated in the literature [35][36][37][38]. For instance, the effect we are describing has been discussed in the theoretical investigation presented in ref. [35] in the context of background for pp → hh → 4b where an effective mistag rate of about 6% is found. On the experimental side, a recent search for a charged Higgs by CMS [38] found that the background from misidentified b-jets dominates over the genuine b-jets one. More interestingly, in ref. [36] the authors presented a tagging algorithm to discriminate between "1b" and "2b" jets and showed that, combining girth, charged track multiplicity and momentum fraction of the leading b-hadron, it is possible to achieve a sizable tagging efficiency for jets containing -7 -  one b-hadron ( 1b ) while having at the same time a large rejection of jets containing two bhadrons (¯ 2b ). In our analysis we adopt their tagging strategy but, instead of implementing directly their tagging algorithm, we simply adopt one average working point and re-weight signal and background events accordingly (details of this procedure are discussed in appendix A). Note that, according to ref. [36], the tagging algorithm becomes more efficient with increasing p T of the jets. Since our cuts enforce the presence of very high p T b-jets, we expect that our implementation is conservative.
In the next two subsections we discuss the reach of two possible search strategies based on tagging at least four and five b-jets. Both analysis strategies that we propose require b-tagged jets with ∆R > 0.5, |η| < 3, p T > p cut The actual values of the p T and H T b cuts are optimized for a given choice of vectorlike quark and heavy Higgs mass. Typical values are p cut T ∈ [100, 300] GeV and H cut T b ∈ [1, 3] TeV. We do not require the four (five) b-tagged jets to be the four (five) highest-p T jets in the event; in fact, the b-tagging efficiency decreases at large p T . We consider the four (five) highest p T b-jets in each event in order to take advantage of the combinatorics (all our signals have at least six b-jets) and thus achieve a larger signal efficiency. For the 4b and 5b analyses we adopt ( 1b ,¯ 2b ) = (0.8, 1/3) and ( 1b ,¯ 2b ) = (0.9, 1/2), respectively. In the 5b analysis we additionally require cuts on the invariant mass of two and three jets.
Finally, let us comment on the reason for adopting H T b rather than the total hadronic transverse energy H T . Large H T multi-jet backgrounds are very difficult to simulate if the minimum jet transverse momentum (p cut T ) is much smaller than H cut T . In fact, configurations with multiple relatively low-p T and well separated 4 jets are enhanced with respect to those 4 Splitting a jet leads to a much larger HT only if the two resulting jets have large ∆R. Shower generated -8 -

JHEP07(2020)241
with fewer higher-p T jets. As a result one needs to include final states where the number of well-separated jets is roughly up to H cut T /p cut T . Small scale computer simulations allow for up to 6 jets in the final states implying that QCD background with H cut T 1 TeV and p cut T 100 GeV is dangerously sensitive to higher order corrections. A possible solu tion to this problem is to refrain from a calculation and use signal-depleted regions to measure the multi-jet background. Another strategy, which we follow here in order to give an idea of the sensitivity that this kind of analyses can yield, is to replace H T with H T b . Finally, an alternative method is to perform a multivariate analysis which uses the transverse momentum of all identified b-jets as input.

4b analysis
The dominant backgrounds to an analysis based on four high-p T b-tagged jets are (2b2j, 4b, 4t, 2b2t) + n-jets (n = 0, 1). In figure 5 we show the H T b distribution of backgrounds and 6b final state signals for m b 4 = 1.5 and 2.5 TeV. The panel on the right shows the impact of adopting the 1b/2b tagging strategy, whose main effect is to suppress the 2bjj(j) background by an order of magnitude. In figure 6 we show the distributions of H T b versus the transverse momentum of the fourth-highest p T identified b-jet, P T (b 4 ), for the m b 4 = 1.5 TeV signal and 4b + n-jets background. The requirement of large H T b implies the presence of several high-p T b-tagged jets. The softest identified b-jet is usually much harder in our signal (due to the presence of heavy resonances) rather than in the QCD background. Therefore, since in the background large H T b is achieved mostly from the leading jets, cutting on P T (b 4 ) offers a good discrimination.
In tables 1 and 2 we present the result of the (p cut T , H cut T b ) optimization. For a range of m b 4 ,t 4 masses we give the optimal cuts, the fiducial signal cross section assuming 100% branching ratios and fiducial background cross sections. All cross sections are calculated for pp collisions at 14 TeV. The results in these tables can be easily converted into upper limits onto branching ratios for vectorlike quark decays into neutral and charged Higgses following the statistical procedure detailed in appendix B.
Before discussing the sensitivities we obtain, let us comment on the contributions of the 4b + n-jets background with n ≥ 2. We have studied the breakdown of the dominant background 4b(j) into 4b (without additional jets with p T (j) > 100 GeV and ∆R > 0.5) and 4b + j (with p T > 100 GeV and ∆R > 0.5). The latter receives contributions from qg initiated hard parton scattering which, after imposing the cuts mentioned above, accounts for about 60% of the total 4b+j background. The 4b+j cross section stemming from initial parton configurations that are shared with 4b is smaller than the 4b cross section. This suggests that the cuts we consider allow for a perturbative estimate of the relevant backgrounds; in particular, we expect the 4b + n-jets background with n ≥ 2 to be subdominant.
In figures 7 and 8 we present the results that we obtain for b 4 and t 4 cascade decays. In both figures blue (green) bounds correspond to 100 (3000) fb −1 of integrated luminosity at 14 TeV. For the LHC running at 13 TeV and current luminosity of 139 fb −1 we find that the curves would almost overlap with the bounds presented for 14 TeV 100 fb −1 .
jets do not capture these effects since QCD radiation is calculated in the collinear limit.    As anticipated the expected bounds weaken as the number of top-jets in the final state increases. If decays to multiple heavy Higgs bosons are simultaneously open, the upper limit on the total b 4 (t 4 ) branching ratio lies in between the solid and dashed lines but depends on the relative size of branching ratios into H, A and H ± . As an example, the dot-dashed lines show the cases in which the branching ratios into the three heavy Higgses, H, A, H ± , follow the typical patterns for singlet-like b 4 and t 4 , 1/4 : 1/4 : 1/2, assuming degenerate masses m H,A,H ± = 1 TeV. In this case, we also get contributions from processes in which the two vectorlike quarks decay to different Higgses. In particular this allows for the final states 4b2t and 4t2b from decays of b 4 and t 4 , respectively. This is especially relevant for H ± decays in b 4 cascades because the 4b2t final state is superior to the 4t2b one in the context of the all-hadronic analysis we are considering.     Figure 9. Left: H T b distribution of the 6b signal and backgrounds in the 5b analysis. We show the signal distribution for m b4 = 1.5 and 2.5 TeV in blue and cyan respectively. We show the dominant backgrounds, 4bj + 2bjj + 2bjjj in the red shaded region. We show the subdominant backgrounds, 6b and 2b2tj, in the green and yellow shaded regions respectively. Right: H T b distribution of the 6b signal and backgrounds in the 5b analysis after incorporating the 1b/2b-tagging strategy assuming ( 1b ,¯ 2b ) = (0.9, 1/2) efficiencies.    Figure 11. Expected 95% CL upper limits on the t 4 branching ratio into neutral or charged Higgses in the 5b analysis. Styling of all curves matches that in figure 8.

5b analysis
The dominant irreducible backgrounds to an analysis based on five high-p T b-tagged jets are 6b, 6t, 4b2t and 4t2b. On top of this we have contributions from (4b, 4t, 2b2t) + jet and (2b,2t) + 3 jets where the jets are mistagged. Due to the sizable mistag rate coupled with a large diagrammatic multiplicity, the latter sources of backgrounds are by far dominant. The situation can be ameliorated by adopting the additional 1b/2b tagging strategy. As mentioned above, in the 5b analysis we are guaranteed to be able to reconstruct one full vectorlike quark decay. For each choice of the heavy Higgs and vectorlike quark masses we additionally optimize the significance by requesting that the invariant masses of at least one pair and one triplet of b-jets (obtained by adding one of the other jets to the first pair) reconstruct the signal decay chain. For signal with larger masses we find optimal cuts with |m bbb − m b 4 | < 1 TeV and |m bb − m H | < 400 GeV, where the cut on the triplet invariant mass slightly decreases for smaller m b 4 .
In figure 9 we show the H T b distributions of signal and backgrounds. In tables 3 and 4 we present the result of the H cut T b optimization. The upper bounds on the total branching ratios of vectorlike quarks into heavy Higgses are presented in figures 10 and 11. The discussion of the tables and the bounds parallels the one presented in the previous section.

Discussion and summary of results
In the 5b analysis we find that the LHC with 3 ab −1 of integrated luminosity is sensitive to the b 4 → H(A)b → bbb decay for masses of b 4 up to about 2.2 TeV and masses of neutral Higgs bosons up to about 2 TeV (taking into account that the H → bb branching ratio is slightly less than 100% due to the presence of additional H decays). The same The 4b analysis has a slightly lower reach, but the lower b-jets multiplicity allowed us to perform an in depth study of multi-jet QCD backgrounds and to demonstrate that they can be brought under control. In particular, we have studied the breakdown of 4b + j and 2b + jj(j) backgrounds and their contributions to the 4b final state after shower effects. We found that a significant effect originates from jets in which a gluon splits into bb pairs. However, separation of genuine b-jets and 2b-jets can be used reduce the contribution of QCD mulit-jet backgrounds by about an order of magnitude (see figure 5) with only a mild effect on a genuine 4b signal. This makes the analysis robust against backgrounds with high jet multiplicity, while still remaining conceptually simple.
In both analyses, we incorporate a 1b/2b tagging efficiency that can be obtained in a simple cut-based analysis [36]. Though, improvements in the reach of vectorlike quark and heavy Higgs masses can be achieved with stronger rejection of 2b-jets. In figure 12 we show the impact of the 1b/2b tagging strategy on the reach of the 6b signal. Purple and green curves are the results obtained without and with the adoption of the 1b/2b discrimination.  Table 4. Signals pp → t 4 t 4 → 6t, pp → t 4 t 4 → 4b2t and background fiducial cross sections in units of fb in the 5b analysis. For each t 4 mass hypothesis, the optimized cuts are listed in the last row for integrated luminosity L = 3 ab −1 . The signal efficiency of the cuts ( ) is also presented. All results are presented using the 1b/2b-tagging working point ( 1b ,¯ 2b ) = (0.9, 0.5).  Figure 12. Left: impact of 1b/2b tagging efficiencies on the 4b analysis. In purple, we show the expected 95% CL upper limit on the b 4 → Hb → bbb branching ratio with m H = 1 TeV assuming no discrimination between genuine b-jets and bb-jets. In green we show the optimal tagging efficiencies used in the main results. We show the effect of perfect discrimination of 2b-jets in Cyan. Right: corresponding study of tagging efficiencies in the 5b analysis.

JHEP07(2020)241
The blue curves show the limiting case of perfect 2b rejection. Note that, for the 5b analysis, improvements in the 1b/2b separation have the potential to considerably improve the bounds. In any case, rejection of 2b-jets is critical for our confidence in the analyses and efficient rejection of multi-jet QCD background.

Conclusions
We discussed search strategies for cascade decays of pair-produced vectorlike quarks through charged or neutral Higgs bosons that lead to 6 top or bottom quarks in final states. In extensions of the type-II 2HDM with vectorlike quarks, such cascade decays can easily dominate. Depending on model parameters and especially on the hierarchies in the masses of new quarks and Higgs bosons, many decay chains are possible, and several of them can be simultaneously sizable. Among them, t 4 → H ± b → tbb and b 4 → H(A)b → bbb are especially well motivated, since the corresponding branching ratios can be close to 100% for any medium to large tan β, even if all the decay modes are kinematically open. The suggested search strategies, focusing on the 6b final state, are tailored for b 4 → bH → bbb process. However, we also find them very effective for t 4 → H ± b → tbb and other possible cascade decays of t 4 and b 4 . While the signals we consider are subject to large reducible backgrounds originating from QCD multi-jets, we find that a significant portion of these backgrounds can be rejected by separating genuine b-jets from those emerging from parton shower effects. This aspect of our analysis gives us confidence that our results are robust against QCD backgrounds with high jet multiplicity which otherwise would pose a formidable challenge to properly estimate.
Among the main results we find that the Large Hadron Collider with 3 ab −1 of integrated luminosity is sensitive to b 4 → H(A)b → bbb for masses of b 4 up to about 2.2 TeV and masses of neutral Higgs bosons up to about 2 TeV. The same analysis, applied for t 4 → H ± b → tbb, leads to a similar reach: the masses of t 4 up to about 2.2 TeV and the charged Higgs boson up to almost 2 TeV can be explored at the LHC. The sensitivity of the analysis based on the 6b final state gradually loosens as the number of top-jets increases. However, we find that even for the 6t final state from pair produced t 4 the sensitivity still extends to t 4 masses of about 2 TeV.
With already existing data accumulated at 13 TeV LHC corresponding to 139 fb −1 of integrated luminosity we estimate that the presented search strategy is sensitive to vectorlike quark masses up to about 1.8 TeV and charged and neutral Higgs bosons up to about 1.6 TeV.
The reach of the presented search strategy significantly exceeds the reach of separate searches for vectorlike quarks and charged and neutral Higgs bosons. Especially the sensitivity to the charged Higgs boson, extending to about 2 TeV, stands out when compared to models without vectorlike matter. Thus searching for bottom-rich events presents a unique opportunity to simultaneously discover vectorlike quarks and heavy charged and neutral Higgs bosons at the LHC.
For the 4b analysis, the additional efficiencies due to the 1b/2b discrimination can be estimated as: For the 5b analysis, the additional efficiencies due to the 1b/2b discrimination can be estimated as:

B Poisson statistics
For completeness we present the statistical framework we adopt to extract the bounds. Following ref. [39], the upper limit (at confidence level of 1 − α) of the number of signal events that can be extracted from observing n events over an expected number of background events b is: where F χ 2 is the χ 2 cumulative distribution. It is easy to show that, assuming no signal and Poisson distributed background, the median of the upper limits calculated over a sample of the background distribution is given by: where L is the integrated luminosity and ε sig is the signal efficiency of the adopted cuts. Finally, the signal cross sections that we consider are proportional to the square of the vectorlike quark branching ratio for which we calculate the expected upper limits (e.g. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.