Scale-invariant resonance tagging in multijet events and new physics in Higgs pair production

We study resonant pair production of heavy particles in fully hadronic final states by means of jet substructure techniques. We propose a new resonance tagging strategy that smoothly interpolates between the highly boosted and fully resolved regimes, leading to uniform signal efficiencies and background rejection rates across a broad range of masses. Our method makes it possible to efficiently replace independent experimental searches, based on different final state topologies, with a single common analysis. As a case study, we apply our technique to pair production of Higgs bosons decaying into $b\bar{b}$ pairs in generic New Physics scenarios. We adopt as benchmark models radion and massive KK graviton production in warped extra dimensions. We find that despite the overwhelming QCD background, the $4b$ final state has enough sensitivity to provide a complementary handle in searches for enhanced Higgs pair production at the LHC.


Introduction
Jets are a ubiquitous component of the LHC program, relevant for precision Standard Model measurements, Higgs boson characterization and Beyond the Standard Model searches [1,2]. In particular, searches for New Physics in multijet events are an important element of the LHC physics program. New resonances and contact interactions have been searched for by ATLAS and CMS in final states with two jets [3][4][5][6][7][8][9][10], four jets [11,12], six jets [13][14][15], eight jets [16] and up to ten jets for the semi-classical black holes searches [17,18].
A challenge in searches for new phenomena in multijet final states is the prohibitively large QCD multijet background. A range of techniques is then required in order to identify particular categories of jets, making it possible to reduce this background. Among those that have been validated and applied to searches, one can mention b-tagging [4,19,20], jet shapes for quark/gluon and other flavour identification studies [21,22] and jet substructure tools [23][24][25][26][27][28][29]. Stringent constraints on a variety of new physics models have been obtained this way, with many more expected with the full 2011-2012 dataset.
From the kinematic point of view, the most common scenario is that of a heavy resonance X produced in the s-channel which then decays back into a pair of quark or gluon jets. However, there is a large class of models where paired production of resonances dominates, processes of the form pp → X → 2Y → 4 partons, with Y being another massive particle. The mediator X of this production might be an exotic particle from a new strongly coupled sector, or a resonance from extra-dimensions, such as a massive graviton or a radion. The Y resonance could be either some BSM particle (sparticles in R-parity-violating supersymmetry, colorons [30], axigluons [31]) or some SM particle (W, Z or Higgs) that subsequently decays into quarks and gluons.
These generic four parton processes lead to very distinct final state signatures depending on the interplay between the masses of the two intermediate resonances, M X and M Y . If the mass ratio is large, M X M Y , the Y resonances will be produced very boosted, and typically the decay products of each of the two Y resonances will be collimated into a single fat jet. On the other hand, for M X ∼ 2M Y , the Y resonances will be produced nearly at rest, decaying into four well separated jets. Existing searches assume either the highly boosted or fully resolved regimes, and by doing so exclude a potentially large region of the New Physics parameter space.
It is the goal of this paper to design a jet reconstruction and analysis strategy that can be applied simultaneously to the boosted and resolved regimes. This will be achieved by merging the boosted-regime strategies, based on jet substructure techniques, with a suitable strategy for the resolved four-jet regime, based on dijet mass pairings, together with a smooth interpolation between the two limits. Such a strategy has the potential to make the experimental searches more efficient and allow a wider range of BSM models to be probed within the same common analysis.
The approach that we will present here is fully general and model-independent, assuming only that resonances are pair produced and then decay hadronically, with no constraint on the absolute masses: indeed, at the parton level the problem turns out to be scale independent, and the dynamics are completely determined by the mass ratio r M ≡ M X /2M Y . Of course additional QCD radiation and confinement on the one hand and experimental cuts on the other break this scale invariance, but we will see that the general qualitative results are robust.
To provide a realistic application of our technique, we will examine resonant Higgs pair production, recently studied as a promising probe of New Physics scenarios [32][33][34]. We will therefore derive model independent limits on BSM resonant Higgs pair production in the 4b final state. We will then apply these bounds in the context of warped extra dimensional models, where Higgs pair production is mediated by either a spin zero (radion) or spin two (massive Kaluza-Klein graviton) resonance. We will show that a wide range of the parameter space of the radion and massive KK graviton scenarios can be covered by present and future LHC data, and that despite the overwhelming QCD background, the 4b final state has enough sensitivity to provide a useful handle in searches for enhanced Higgs pair production at the LHC.
The outline of this paper is as follows. We begin in Sect. 2 by introducing the general search strategy for pair produced resonances that can be applied simultaneously to the boosted and resolved regimes. In Sect. 3 we review the theoretical models for resonant Higgs pair production in warped extra dimensions scenarios. Then in Sect. 4 we apply the jet reconstruction strategy both to signal events and to the QCD multijet background, and explore the potential for new physics searches in the 2H → 4b channel. In Sect. 5 we conclude and outline future developments.

Scale-invariant resonance tagging
Multijet signatures have long been recognized as an important channel for Beyond the Standard Model searches at hadron colliders [35]. The main difficulty in these channels is how to tame the overwhelming QCD multijet background. Searches in multijet final states are commonly separated into boosted and resolved regimes. An example of the former arises when light partons are produced from the decay of a heavy resonance. Recently developed jet substructure techniques, reviewed for example in [36,37], make it possible to substantially improve the discrimination power in the boosted regime. At the LHC the advent of jet substructure methods has made it possible to study boosted production of the heavy Standard Model particles, like W and Z bosons and top quarks, for √ŝ , the centre-of-mass energy of the hard process, above 1 TeV. Searches for new physics such as resonant production of V V or tt [38][39][40][41] or searches for boosted supersymmetric particles and colored scalars [11,15] have also benefited from these developments.
In this section we introduce a general strategy for jet reconstruction designed for searches of pair-produced resonances in fully hadronic final states, which is simultaneously suitable for both the highly boosted and the fully resolved regimes and that smoothly interpolates between them. The generic process we are interested in is the s-channel production of a resonance X which then decays into a pair of resonances Y , which in turn each decay into a pair of light Standard Model particles, labeled z, The ratio between the masses of the X and Y resonances will determine the degree of boost of the Y resonances and consequently the angular distribution of their decay products z that will be observed in the detector. At parton level, neglecting the mass of the final state particles m z , the problem is scale invariant and is characterized by a single dimensionless variable, denoted by which is simply the boost factor from the Y rest frame to the X rest frame. In the highly boosted regime, r M → ∞, while in the fully resolved regime, where the intermediate resonances Y are produced at rest, r M = 1. Schematic diagrams for the boosted and resolved topologies are shown in Fig. 1. If we assume that the heavy resonance X is produced at rest, so that the laboratory and center-of-mass reference frames coincide, we can parametrize the four momenta of the X → Y Y decay with the convention that P = (p T,x , p T,y , p L , E). We then have where Y 1 and Y 2 are the two decay products of the X particle, θ * Y is the angle of Y 1 with respect to the beam, and φ * Y is the azimuthal angle. The boost parameters from the laboratory frame to the rest frame of the Y particles are given by As one can see, the boost of the Y particles, β Y , is independent of the absolute masses of X and Y , and depends only on their ratio. It is in this sense that we can consider that the problem at hand is scale invariant: γ Y does not depend on any absolute mass scale. Fig. 1 suggests that depending on the value of the mass ratio r M the search strategy should be different. For large r M , the resonances Y will be very boosted, and thus the angular distances of their decay products will be small, while for low r M the four final state particles will be well separated. Since we are mostly interested in the case in which the final state particles z are QCD partons, quarks or gluons, we will end up either with two fat jets (in the boosted regime), four well separated jets (in the resolved regime) or one fat and two separated jets (in the intermediate regime). Fat jets are jets for which the substructure pattern is unlikely to have arisen from QCD radiation.
Given that in general we do not have information on the masses of the intermediate resonances, we don't know a priori in which of the two regimes we will find ourselves, and it would be beneficial to have a search strategy that simultaneously explores all possibilities. It should exhibit reasonably homogeneous efficiencies and background mistag rates for any value the mass ratio r M within the physically allowed range. Below, we will present such a combined strategy that simultaneously explores the boosted and resolved regimes.
In order to validate the performance of the strategy that we will propose we have generated events for the generic process Eq. (1) using a toy Monte Carlo simulation. The heavier resonance X is assumed to be produced at rest in the laboratory frame, justified by the fall-off at large masses of parton luminosities [42], and to decay into the two intermediate resonances Y with a homogeneous angular distribution, as if it were a spinzero particle. The massless decay products of the Y resonance decay are also assumed to decay isotropically in the Y rest frame. In this toy simulation the possible widths of the intermediate resonances are neglected, as well as the masses of the final state particles m z .
In view of the later applications to Higgs pair production, we will set M Y = 125 GeV and vary M X in a wide range, although it should be clear that at parton level the event classification will depend only on the ratio r M . For each of the M X values in the range from 250 GeV (resolved regime, r M = 1) to 5 TeV (highly boosted regime, r M = 20), we have generated 50K toy MC events.
To study the performance of the jet reconstruction strategy for a realistic collider environment, the parton level events from the toy Monte Carlo have been showered and hadronized with Pythia8 [43], version 8.170. We have done this for LHC centre of mass energies of 8 and 14 TeV, and we include also underlying event and multiple interactions with the default tune 4C of Pythia8. Initial state radiation has been modeled assuming that the resonance is produced in the gluon-gluon channel 1 . Parton and hadron-level events were then clustered with the anti-k t jet algorithm [44] with radius of R = 0.5. Such small radii (R = 0.5 for CMS, R = 0.4 for ATLAS) are used in most experimental multijet analysis.
No additional cuts will be applied to the reconstructed jets at the parton level, so as to avoid introducing any explicit breaking of scale invariance. On the other hand, at hadron-level it becomes necessary to introduce additional kinematic cuts, which explicitly break scale invariance. In this section we will adopt the following set of basic kinematic cuts for jets in hadron-level events:  In the H T variable the sum goes over the four leading jets of the event above the p min T cut. These loose cuts have a very limited effect on the selection efficiencies except at the smallest values of r M . 2 In order to identify the three different regimes, boosted, resolved, and mixed, useful information is provided by evaluating the fraction of events with a given reconstructed jet topology. We show the relative fractions of the different jet topologies in Fig. 2 for both parton and hadron-level, as well as the sum of events with two, three and four jets. At parton level, by construction, events can only have between two and four jets, and we clearly see that four-jet events dominate at low r M , two-jet events in the boosted regime for large r M , with the three-jet case in between. 3 At hadron-level both the shower and the kinematic cuts break scale invariance, and now we can have events with fewer than two jets, at very low r M , due to the basic cuts Eq. (6), and with more than four jets due to parton radiation. Note though that still between 30% and 50% of the events, depending on the value of of r M , have between two and four jets.
Given that parton showering can significantly modify the number of jets, an event classification based on the number of jets is not optimal under realistic conditions. Instead, we use an alternative classification, based on the number of tagged jets per event, that is, jets that are found to have non-trivial substructure. We will proceed as follows: each of the two hardest anti-k t jets in the event is reclustered using the Cambridge/Aachen algorithm [45] with R sj = 1.3 (where the subscript sj means sub-jet), 4 and processed with the BDRS mass-drop tagger [46]. This tagger has two parameters µ and y cut . To determine if a jet arises from a massive object, the last step of the clustering for jet j is undone, giving two subjets j1 and j2, with m j1 > m j2 ; if both are significantly lighter than the parent jet, m j1 ≤ µ · m j and the splitting is not too asymmetric where ∆R j1,j2 is the angular separation between the two subjets, then j is returned as the tagged jet. Otherwise we replace j with j1 and apply the unclustering to the new j, repeating the procedure until we find a subjet for which the mass-drop and asymmetry conditions are both satisfied. If the procedure recurses to the point where it finds a single-particle jet, then the jet is considered untagged. We use the values µ = 0.67 and y cut = 0.09, as in the original BDRS paper [46]. Our strategy is separated into two parts: the analysis chain, which sets the flow of the event classification, and the quality requirements, which determine whether a given topology is classified as a signal event or rejected as a background event. We discuss these two parts in turn.
Analysis chain. We start by examining events with at least two jets after basic cuts. Summarized in the flow chart of Fig. 3, the analysis chain depends on the number of mass-drop jet tags, that is, jets that have been identified by the BDRS tagger to have an internal structure potentially not arising from QCD radiation.
• If the two hardest jets in the event are mass-drop tagged, we examine if these two jets can be identified as arising from the decay products of two boosted Y resonances. This is established by verifying if the two Y candidates satisfy the quality conditions on their mass difference and angular separation listed below, in which case the event is assigned to the 2-tag sample.
• If the event has a single mass-drop tag among the two hardest jets, or if the event had two mass-drop tags, but was not assigned to the 2-tag sample, then we examine whether the event can be classified as having an underlying three-jet topology, where the decays of one Y resonance are collected into a single jet but not those of the other. Events with fewer than three jets after cuts are discarded. If there is a single mass-drop tag, the second Y candidate is formed by adding the four-vectors of the other two hardest jets in the event. If there are two mass-drop tags but the event has been rejected in the 2-tag category, we examine combinations whereby one of the tagged jets is taken to correspond to a first Y candidate, while the other tagged jet is assumed to be a mistag and is combined with a third jet to make up the second Y candidate. If the jet mass and angular quality requirements listed below are satisfied, the event is classified into the 1-tag sample.
• If no mass-drop tags are found in the event, or tags have been found but the event has failed to be assigned in either of the above categories, we examine the possibility of an underlying four-jet parton kinematics. Discarding events with fewer than four jets passing the basic cuts, we select the jet pairing such that the combination ij and kl of the five hardest jets in the event leads to jet masses M ij and M kl that minimizes the difference |M ij − M kl |, and use this pairing to reconstruct the two Y candidates. 5 If these two candidates pass the mass and angular quality requirements given below, the event is classified as belonging to the 0-tag sample. 6 Quality requirements. To identify the event as arising from the decay of the X resonance, Eq. (2), additional mass and angular quality conditions are required, which are essential to further suppress the QCD background. Some of these requirements are designed so as to apply similar conditions to both the boosted and resolved topologies.
1. We require the masses of the two Y candidates to be the same up to a given mass tolerance f m , to account for experimental mass resolution, as well as mass smearing due to underlying event, hadronization and initial and final-state radiation: where m Y is the average mass of the two reconstructed Y resonances. We assume in this work a fixed value 7 for the detector mass resolution f m of 15% [4,6,47]. This requirement cannot be made too stringent otherwise a large fraction of signal events would be discarded.
2. In the case in which the mass of the Y resonance is known, the masses of the two Y candidates must lie in a mass window around M Y , where the width of the window is determined by the mass resolution of the detector.
Since we will be considering Higgs pair-production, we will set M Y = M H = 125 GeV in the following, though this requirement has a small impact in signal events, and is only relevant to suppress the QCD background.
3. The separation in rapidity of the two Y candidates must be smaller than some upper value, motivated by the fact that for a given mass of the Y 1, Y 2 system, background events, dominated by t-channel exchange, are enhanced in the forward region, while signal events, dominated by s-channel exchange, tend to be more central. We will take ∆y max = 1.3 in the following, a value optimized from the high mass dijet searches at the LHC [4,6].
4. Likewise, the separation in rapidity between the two jets of a Y candidate in the resolved case, y Y i,1 and y Y i,2 , with i = 1, 2, must also be smaller than some upper value, possibly different from before, since for these kind of topologies, signal events will be produced closer in rapidity that QCD multijet production. We will take ∆y res max = 1.5 in the following, and we discuss below the rationale for this choice.

5.
To prevent excessively asymmetric configurations, whenever we have two resolved jets that correspond to a given Y -candidate, one with p T and the other with p T (in either the 1-tag sample or the 0-tag sample), we require This cut plays a similar role as the asymmetry requirement in the BDRS mass-drop tagger, Eq. (7), but now in the case of resolved jets, and it helps reject events where a soft jet arises from final-state radiation (FSR). 8

6.
With a similar motivation, for each two resolved jets in a Y -candidate with mass m Y i , with i = 1, 2 we impose the mass-drop condition on the masses of these two resolved jets, m Y i,1 and m Y i,2 , as follows where µ is the same parameter as in the BDRS mass-drop tagger. Together with the asymmetry condition above, applying the mass drop requirement also in the resolved jets ensure that the same conditions hold for the three different possible topologies, from the highly boosted to the fully resolved regimes. Note, however, that for our default choice of µ = 0.67, the mass-drop cut has only very limited impact on the final reconstruction efficiency. 8 To see the equivalence with the cut of Eq. (7), note that for reasonably small ∆R, the mass of the Y candidate is m 2 T ∆R 2 12 . Making use of the fact that p (2) T , Eq. (7) reduces to p Jet Reconstruction R R sj R f n filt µ y cut 0. 5  The values of the parameters used in our implementation of the jet reconstruction strategy are summarized in Table 1.
In addition, in order to improve on resolution, jet masses are filtered [46] as follows: the constituents of each tagged jet are reclustered with a smaller radius R filt = min (∆R sj,sj /2, R f ), with R f = 0.3 and ∆R sj,sj the angular distance between the two subjets after mass-drop in the boosted case. Then only the three hardest subjets, n filt = 3, are retained to account for at least one QCD emission. The filtering procedure improves mass resolution of the reconstructed resonances [46,48] and makes the procedure more resilient to soft radiation from the underlying event and pile-up [49]. 9 This jet reconstruction strategy has been implemented in a code based on FastJet3 [51], and we have processed the parton and hadron-level toy Monte Carlo events through it. We show in Fig. 4 the efficiency of the resonance pair tagging algorithm as a function of resonances mass ratio r M for the parton and hadron-level toy Monte Carlo events. We show both the total efficiency and the breakup of the efficiencies corresponding to the 2-tag, 1-tag and 0-tag samples. The impact of the moderately loose selection cuts Eq. (6) on the parton-level efficiencies is negligible, and thus the differences between parton and hadron-level arise from initial and final-state radiation.
At parton level, at low r M , the 0-tag sample dominates as expected from the resolved regime, while for large r M , the boosted regime, it is indeed the 2-tag sample that dominates. The 1-tag sample is important at intermediate boosts. The combined efficiency is found to be rather flat in all the mass range, between 30% and 40% for all mass values, showing that we are able to obtain a reasonable tagging efficiency irrespectively of the degree of boost of the X resonance decay products. At hadron-level efficiencies are somewhat lower due to additional parton radiation and underlying event at low masses, but still we obtain a reasonable tagging efficiency of between 20% and 30% in all the relevant range, approximately constant for all topologies, except close to r M = 1. The 1-tag sample Let us mention that the production threshold region close to r M ∼ 1 is challenging from the jet reconstruction point of view. First of all, there will be a substantial degree of overlap between the decays products of the two Higgs bosons, since the two are at rest, which leads to wrong mass pairings. Second, it is quite frequent that large-angle initial-state radiation (ISR) appears as additional jets, again confusing the pairing of the original jets.
To help understand the efficiencies that we find, let us recall that the asymmetry cut y cut in the BDRS mass-drop tagger leads to an upper bound on the efficiencies of signal events in the boosted regime of approximately ∼ (1 − 2y cut /(1 + y cut )) (this result is exact for the two-prong decay of a highly boosted scalar). For the two tag sample at large r M , we therefore expect that at parton level the tagging efficiency is given by for our choice of parameters, where the last factor accounts for the contribution to the total efficiency from the cut in ∆y between the two Y -candidate jets. 10 This is exactly what is obtained in Fig. 4. At hadron-level the efficiency in the boosted regime is somewhat smaller due to the contamination from initial-state radiation and the underlying event.
When we have a resolved Y resonance candidate, the separation in rapidity between the two jets must be smaller than some upper value, ∆y res max . To determine the value of this cut, we note that in the small R limit, if such a cut is the only one applied to the final state, the efficiency of the fully resolved case for r M ∼ 1 is given by  so demanding that the efficiencies at low r M match the asymptotic large r M value, Eq. (14), we obtain ∆y res max = ln for the default value of y cut used in the mass-drop tagger algorithm. With this choice, we can achieve at low r M the same efficiency as at large r M , at least in the parton-level case without the basic kinematic cuts. Note that this cut ensures not just uniform signal efficiency, but it also is useful for background rejection, especially in scenarios where M Y /2 is substantially larger than the jet p min T cut. A final interesting comparison is that of the efficiencies between parton and hadronlevel, to gauge the robustness of our event classification based on giving priority to the mass-drop tags. This is useful in order to understand the impact of parton showering and underlying event, as well as of the basic kinematic cuts in the tagging of the heavy resonances. Let us recall that the only difference in the analysis chain between parton and hadron level events are the basic cuts in Eq. (6), and the fact that in the 0-tag case we study the mass pairings of the five leading jets. Results are shown in the two tag sample is very similar at parton and hadron-level, for all the values for r M . The efficiencies for 0 and 1-tags sample have a similar shape but a smaller magnitude, and the shape is somewhat shifted down to lower r M values. This small shift between the parton-level and hadron-level efficiencies is perhaps attributable to transverse boosts induced by initial-state radiation.
Since the jet reconstruction strategy that we advocate is approximately scale invariant, one also expects the results to be reasonably independent of the jet radius R used in the jet clustering: while the relative fraction of 2-tags, 1-tag and 0-tag events will of course vary with R, their sum should be stable. Indeed, we show in Fig. 6 the total efficiency in parton and hadron-level events for three different radii, R = 0.3, 0.5 and 0.8. At parton level, results in the boosted regime are strictly R-independent, as shown in Fig. 6. Except at very low masses, parton-level results are R-independent in all the mass range. For low r M , the degradation at parton level with increasing R arises in part because the likelihood that the decay products from different Y resonances end up in a single jet is higher for larger R. Also at hadron-level the total tagging efficiency is reasonably independent of R.
To summarize, in this section we have presented our general strategy for a resonance reconstruction analysis that can be applied simultaneously to the boosted and the resolved regimes, with a smooth transition between the two limits. It is clear however that some of the details of the strategy can be modified without affecting the general philosophy. One could study different ways of dealing with the four-jet events rather than selecting the pairings which minimize the relative dijet masses, like cuts in the angular distributions. It is also possible to extend the number of jets considered to build the resolved Y -candidate in the 1-tag case up to the fourth or the fifth jet, in analogy with the procedure used for the 0-tag case.
These modifications could lead to an overall improvement of the tagging efficiency, but the basic strategy would be left unaffected. Finally, other substructure taggers could be used to classify events, such as N-subjettiness [52] or pruning [53] among many others (see [54] for a recent systematic comparison). Note however that those taggers with an asymmetry cut, like mass-drop and pruning, are special, because that cut can be linked with 3-and 4-jet analysis parameters, as done in the present analysis. In this respect, N-subjettiness is quite different, because it is cutting on the radiation pattern in the jet.

Resonant Higgs pair production in warped extra dimensional models
Now we discuss the benchmark models that we will consider for resonant Higgs boson pair production. These models are based on the warped extra dimensions scenario [55], where Higgs pair production is mediated by either a Kaluza-Klein (KK) graviton or by a radion. We will assume that the Higgs is the Standard Model boson [56], and consider its dominant decay into two pairs of bb quarks. Higgs pair production in the Standard Model has a small cross section [57] (approximately 18 fb at 14 TeV), but larger rates can be expected in New Physics models [34,58] like supersymmetry, composite models, and warped extra dimensions. With this motivation, in this section we review the theoretical expectations for resonant Higgs pair production in the context of warped extra dimensional models, keeping in mind that the strategy proposed in this paper is equally valid for any other Higgs pair production scenario. Due to Bose symmetry, only resonances of spin zero and spin two can decay on-shell into a pair of Higgs bosons. Both types are present in models with warped extra dimensions. They are referred to as radion and KK-graviton, denoted by φ and G respectively. These models can naturally explain the large hierarchy between the Planck and electroweak scales by introducing a nontrivial geometry in the extra dimension. The background metric for the case of a single extra dimension is given by where y refers to the coordinate in the 5th dimension and k is related to its curvature. The so-called ultraviolet (UV) and infrared (IR) branes are introduced at y = 0 and y = L, respectively. Depending on the scenario, SM fields can be localized in the IR brane or be allowed to explore the 5th dimension as well. At each position in the extra-dimension y * , the local cutoff is given by [59] Λ(y * ) = e −ky * Λ(y = 0) . g µν = e −2ky η µν → e −2(ky+F (x,y)) (η µν + G µν (x, y)) .
The fluctuation of the size of the extra dimension y is described by the 4D scalar radion field, denoted here by φ: where φ(x) is the 4D wave function and e 2ky the localization profile in the fifth dimension.
The fluctuations of 4D space-time are described by the graviton field G µν (x, y). The massless zero mode of this field corresponds to the usual graviton. The first massive excitation, the lightest KK-graviton which will focus on, is where J 2 is a Bessel function and m G the graviton mass. The mass of the KK-graviton is related to k/M P and to the ultraviolet mass scale of the theory Λ G by where x 1 = 3.83 is the first zero of the Bessel function J 1 . The three parameters, k, the UV mass scale Λ G and the reduced Planck mass M P are related by Λ G = e −kL M P . The mass scale Λ G is expected to take a value in the few TeV range. We neglect the effect of localized kinetic terms, which would change the value of x 1 [60]. We show the radion/graviton localization profiles in Fig. 7: both fields are localized towards the infrared brane, but the graviton localization is stronger.
The radion and KK-graviton couplings to the SM particles are fixed by the action where M, N are 5 D indexes and T M N is the 5D energy-momentum tensor involving all fields. After dimensional reduction, the effective coupling between the radion and KKgraviton lightest mode and the SM is given by: where the T i µν are the four-dimensional energy-momentum tensors of the Standard Model species i = b, f, V, H, ..., and V denotes a generic gauge boson. Here we are neglecting corrections depending on the fermion localization parameters c L and c R , that are small when the fermions are localized in IR brane such as the top quark and they are not large for the b quark [61]. It is important to notice that the radion couples to the trace of the energy-momentum tensor, which vanishes at the classical level for massless gauge fields. Note also that the radion scale is related to the KK-graviton scale by Λ φ = √ 6Λ G [62]. The coefficients c i and d i are proportional to the wavefunction overlap of the graviton/radion and the SM fields. For example, the Higgs is IR-localized, as its vacuum expectation value and mass are IR effects. The fact that the KK-graviton wave-function, Eq. (21), is more peaked towards the IR brane than the radion, Eq. (20), as shown in Fig. 7, translates into a stronger coupling of the graviton to the Higgs, beyond what is expected from the trivial rescaling of Λ φ and Λ G .
In the original Randall-Sundrum model (RS1) all SM fields are localized in the IR brane, so all couplings are c i O(1). More realistic models, consistent with experimental constraints, must have SM fields in the bulk, leading to different values of the couplings c i . A well motivated configuration, which we will refer to as bulk RS [63], predicts that the SM fields communicating to the EWSB sector are peaked towards the IR brane. This is the case of the Higgs and longitudinal W and Z bosons, and possibly the top quarks. Light fermions would be localized near the UV brane, whereas massless gauge bosons are de-localized. The graviton and radion would then couple preferentially to IR localized fields, namely h, W L , Z L and possibly t as well. The coupling to γ and g is suppressed by a volume factor 1/kL, and the coupling to light fermions (including the quarks in the proton) would be extremely suppressed. In summary, the two scenarios we are going to consider, RS1 and Bulk RS are defined by the following hierarchy of couplings of the KK graviton to SM particles: Bulk RS scenario: As an example, the stress-energy tensor for the Higgs field is given by which results in the following couplings, We do not consider a non-minimal Higgs coupling to gravity, which would require the use of an "improved" energy-momentum tensor [64] and would lead to a Higgs-radion mixing [65,66], since mixing is not relevant in our case where we require M φ ≥ 2M H . The Feynman diagrams relevant for Higgs pair production mediated by a radion φ and a KK-graviton G µν are schematically shown in Fig. 8. In the remainder of this section we will discuss the production rates of the radion and graviton in the RS1 and Bulk-RS scenarios. It is beyond the scope of this paper to review the experimental constraints on the parameter space of these models. A discussion of the implications of recent measurements for limits on extra dimension models can be found in [23,[67][68][69]. Note than when experimental limits arise from decays to photons, leptons and four-fermion operators involving light fermions, they can be interpreted only in the context of RS1 but not in the bulk RS scenario. Bounds on the radion mass as a function of the cutoff scale Λ φ have been compiled in Refs. [70,71].

Production rates at the LHC
The production rates of the radion and the graviton at hadron colliders will depend on the respective couplings to the light quarks and gluons in the incoming protons. In RS1, the graviton couples to light quarks and gluons with the same coefficient O(1), whereas in the Bulk-RS scenario the couplings to light quarks are very suppressed. The coupling of the graviton to gluons in the bulk RS model is given by whereas in RS1 c g = 1.
The coupling of the radion φ to gluons (and to photons) vanishes at tree level due to classical scale invariance. At 1-loop level it arises due to the trace anomaly, which is related to the beta function, and the top quark triangle diagram. We denote by κ φ g the coupling of the radion to gluons defined by where κ φ g is given by [61]: where we have neglected the top loop contributions. The coefficient of the QCD β function is b 3 = 8. The RS1 case corresponds to neglecting the volume suppressed term. Note that as compared to the graviton, for the radion production cross section there is less model flexibility in that the coupling to the gluons is fixed independently of the localization of the SM fields. In the following we will assume that gluon fusion is the dominant process for both radion and KK-graviton production at the LHC. This is certainly true for the radion, and also for the graviton in the bulk RS scenario, where the couplings to the light quarks in the proton are very suppressed. In the narrow width approximation, the production cross section via gluon fusion of a particle X with mass M is given by where here τ = M 2 /s and L gg is the gluon luminosity function.
We computed the production cross section of the processes pp → G, φ at leading order (LO) using Madgraph5 [72]. The results are shown in Fig. 9. We plot the production cross section, where we have factored out the trivial dependence on the coupling to gluons and the scale of dimension-five operators, Λ G and we also show the cross section for specific values of those parameters.
Notice that the KK-graviton production cross section is larger than the corresponding radion cross section due to the fact that the radion coupling to gluons is loop induced, whereas the KK-graviton has tree-level couplings to gluons. Also, the KK-graviton has five degrees of freedom, compared to the single degree of freedom of the radion.

Graviton and radion decays
In RS1, with all the SM fields localized on the IR brane, a heavy graviton would decay democratically to all degrees of freedom. In the bulk RS, the Higgs and fields associated with EWSB are still IR localized, and using the equivalence theorem, one can show that [73] In bulk RS, the width to gluons and photons is suppressed by the effective volume kL, The graviton would also couple to fermions localized near the IR brane. In many models, third generation quarks are pushed towards the IR brane via a localization parameter ν, which is a ratio of a 5D mass term, M f , and the curvature, ν = M f /k. The effect of ν is as follows: for ν = 1/2, the conformal value, the fermion zero mode is delocalized in the extra dimension, as the profile is flat and does not prefer a particular location inside the extra-dimension. For ν > 1/2, the fermion zero-mode will be localized towards the IR brane, whereas for ν < 1/2, the localization is near the UV brane.
The width to tops is given by The branching ratio of the graviton to the Higgs depends on the top localization as In Fig. 10 we show that the maximal branching ratio to a Higgs boson pair is 25% (when the Higgs is 1/4 of the IR degrees of freedom), and quickly decreases as one increases the branching ratio to top quarks (increasing the value of ν t , and therefore the localization towards the IR brane). The dominant decay modes of the radion are into pairs of massive gauge bosons, Higgs bosons and top quarks. Since the couplings are determined by the masses of the final state particles, and these masses arise from the TeV localized Higgs boson, the RS1 and bulk RS couplings are the same at leading order. The corresponding widths (for large m φ ) are: and Hence, for large radion masses the branching fraction to a pair of Higgs bosons is approximately 25%, independent of Λ φ , since the contribution from decays to top quarks can be neglected. We note that for the smaller m φ values that are relevant for phenomenology the decay into top quarks should in principle be taken into account, but in this work for simplicity we will assume that BR (φ → HH) = 25% independent of the radion mass.

Composite duals and model dependence
So far, we have described Higgs pair production via gluon fusion into a radion or KKgraviton in warped extra-dimensions. In this context, the graviton cross sections are larger than the radion by at least an order of magnitude, and there is little room for changing this hierarchy.
To test how robust this prediction is, we would like to approach this model building in extra-dimensions from the point of view of holography. In this approach, models in warped extra dimensions are an analogue computer for strong interactions. This duality between 4D strongly-coupled theories and 5D weakly-coupled theories with gravity was inspired by the AdS/CFT correspondence, but took hold on a more qualitative basis [74] and has been used to build models of QCD [75], technicolor [76], composite Higgs [77], and even condensed matter systems [78], with some success. In this context, the KK resonances, a consequence of compactification, are the holograms of massive resonances due to confinement.
The KK graviton is therefore the dual of a spin-two bound state in a strongly coupled theory, very much like the f 2 of QCD [79]. One could then wonder how different the coupling structure of the f 2 -like resonance would be with respect to the KK-graviton. As was shown in Ref. [80], the couplings we propose in Eq. (24) saturate the possibilities, once Lorentz, gauge and CP invariance are assumed. No other structures are allowed up to dimension-six operators. Hence, our KK graviton analysis can be directly generalized to strongly coupled sectors with spin-two resonances. Moreover, if the new strongly coupled sector participates in the electroweak symmetry breaking mechanism, a sizable coupling to Higgses would be expected.
The dual of the radion would be the dilaton, the Goldstone boson of scale invariance. 11 The dilaton couplings at tree level would be perfectly mimicked by the radion couplings, as the dilaton will couple to the trace of the stress tensor. This can be shown by writing down an effective theory where the dilaton is spurion of the scale symmetry [83]. Within this analysis, the loop contributions to the dilaton to massless gauge bosons will follow the same structure as the anomalies written in Eqs. (31). Therefore, our analysis of the radion couplings is also applicable to a dilaton in a composite sector.
In summary, the structure of couplings we describe for the KK-gravitons and radion would be the same for the bound state duals. The main difference between the analysis in extra-dimensions and composite theories is the strength of the coupling. For example, one could imagine a composite theory where the spin-two resonance is made up of colorless techni-quarks, hence there would be no tree-level coupling to gluons, whereas the dilaton couplings would be determined by the scale symmetry. In this case, one would expect a larger production of dilaton than spin-two resonances.

New Physics searches in the HH → 4b final state
In this section we apply the general resonance tagging strategy presented in Sect. 2 to a particular scenario, namely the resonant Higgs boson pair production pair which then decays into four b-quarks. The results presented here are model independent and can be applied to any generic BSM model with enhanced Higgs pair production [33,34,58], though we will provide an explicit interpretation of exclusion limits in terms of the radion and graviton couplings in the warped extra dimension models of Sect. 3.
First of all, we discuss the Monte Carlo event generation for the signal and background events with MadGraph and Pythia and evaluate the tagging efficiency as a function of r M , to compare with the approximate kinematics of the toy MC used in Sect. 2. We recall that the main differences between the toy MC and MadGraph are that the latter includes the rapidity distribution for the X resonance and the correct angular distributions of a spin-two particle in the case of the graviton. Then we present the different assumptions that underlie our implementation of b-tagging. A discussion of the background rejection capabilities of the tagging algorithm follows, where we show that the combination of the resonance tagging and b-tagging reduces the QCD multijet background by several orders of magnitude. In the last part of the section we present the implications in terms of model independent searches in the HH → 4b final state, and interpret these results in terms of exclusion ranges in the parameter space of warped extra dimension models. 12

Monte Carlo signal event generation
Our benchmark model is s-channel Higgs boson pair production mediated by a radion or a massive Kaluza-Klein graviton resonance in scenarios with warped extra dimensions. We have implemented these scenarios in the Madgraph5 Monte Carlo program [72]. 13 While the 12 The feasibility of the 4b final state to probe BSM resonant pair-production with jet substructure was also investigated in Ref. [94] in the context of composite octet searches. 13 Our results have been generated at leading order only; NLO corrections for resonant double-Higgs boson production have been calculated in the context of the minimal supersymmetric standard model, in the heavy top-mass limit, in Ref. [84] and were found to be substantial, giving a K-factor of order 2. Similar corrections are probably relevant to our radion case. However, given that other aspects of our study are probably not under control beyond a factor of two, e.g. the b-tagging assumptions for the background, we will conservatively not include the NLO signal enhancement. main motivation to study both radion and graviton simultaneously is to cover a wider range of the model parameter space, a useful by-product is to validate the jet finding strategy for two different angular distributions of the decay products. Indeed, from the kinematic point of view the radion and graviton cases are identical (for equal masses) except for the different angular decay distributions of spin-zero and spin-two particles. Note that in the radion case, since the radion is an scalar, the kinematics and angular distribution will be very close to those of the toy Monte Carlo of Sec. 2 used to validate the resonance tagging algorithm, with the only difference arising from the rapidity distributions of the radion.
We have followed Ref. [85] to model the radion couplings to the Higgs boson and to gluons using the FeynRules framework [86]. The implementation of the model has been based on the default MadGraph5 model with effective theory coupling of the Higgs to gauge bosons. In addition to the SM parameters, in the radion model we have four additional parameters: the radion mass, M φ , the ultraviolet mass scale of the theory, Λ φ , the radion-Higgs mixing parameter, ξ, and the compactification scale, kL. These parameters take the value Λ φ = 3 TeV and kL = 35, supplemented by the no mixing condition that reads ξ = 0. The absence of mixing is justified by the fact that the radion masses considered will be always much larger than the Higgs mass. Any modification of Λ φ translates into a trivial rescaling of the total rates.
To simulate the graviton production we have used the standard Randall-Sundrum model as implemented in MadGraph5. Here the relevant additional parameters are only the graviton mass M G and the ultraviolet mass scale of the theory, Λ G , chosen to be Λ G = 3 TeV. As pointed out on the Sect. 3, the mass scales Λ φ and Λ G of the radion and the graviton are theoretically related. However, from the practical point of view we select independently the parameters of the two models in order not to impose additional constraints on the search ranges. For both radion and graviton event generation the narrow width approximation has been assumed.
We have generated events for radion and graviton production for a range of masses between 250 GeV and 3 TeV. Higher masses lead to too small cross sections to be of any phenomenological interest. As in the case of the toy Monte Carlo events, Madgraph5 parton level events were showered and hadronized using Pythia8 with the same settings for underlying event and multiple interactions.
We have already discussed in Sect. 2 the tagging efficiency of the algorithm for the toy Monte Carlo kinematics, both at parton and at hadron level. However, the basic selection cuts did not match those of a realistic experimental analysis. We will use the following selection cuts instead in the rest of this paper: These cuts are inspired by typical trigger and angular acceptances of the LHC experiments [11,12]. We show in Fig. 11 the comparison between the hadron-level tagging efficiencies at LHC 8 TeV between the toy Monte Carlo events and the MadGraph radion and graviton events, as a function of the mass ratio r M . As we can see, the toy MC results agree well with the radion events, which is a non trivial cross-check that event generation is under control. Also, the efficiencies for the radion and graviton are very similar, showing that the spin-zero vs. spin-two angular distributions do not lead to any large differences at the level of the reconstruction. 14 Note that we have generated fewer mass points with MadGraph5 than with the toy MC, hence the somewhat less smooth distributions in the former case.
One significant difference between Fig. 11 and the results of Sec. 2, is the much lower efficiency in the low r M region. It is a consequence of the larger H T cut in Eq. (40) than in Eq. (6), which severely reduces the fraction of tagged events when m X 300 GeV. Insofar as the H T cut is present mainly to limit trigger bandwidth, one could also imagine lowering it and then controlling bandwidth by means of trigger-level b-tagging.

B tagging
The final state that we are interested in includes four b-quarks from the decays of the two Higgs bosons. Therefore, b-tagging will be an important ingredient to improve the signal over background ratio. We have adopted in this study a b-tagging scenario that we expect to be realistic (possibly conservative), inspired by the ATLAS and CMS capabilities [19,20,87,88]. The probability of tagging a b-quark is taken to be f b = 0.75, the mistag probabilities of c-quarks, f c = 0.10, and of light quarks and gluons, f l = 0.03. We apply the b-tagging conditions on the parton level events after showering but before hadronization, that is, we tag b quarks rather than B hadrons. We will require one b-tag per Higgs candidate. In detail it is implemented as follows: • Determine the number of b-quarks within each of the two Higgs candidates' jets. Such candidate jets can be a single anti-k T jet with radius R (in the boosted regime) or a jet composed by the sum of two different anti-k T jets (in the resolved limit).
• A Higgs candidate jet is considered to be b-tagged if it contains at least one b quark with p T,b ≥ p min T,b = 10 GeV. The b-tag efficiency is denoted by f b .
• A Higgs candidate jet which does not fulfill the previous condition, but contains at least one c quark with p T,c ≥ p min T,b , will be b-tagged with a mistag probability f c .
• A Higgs candidate jet which contains only light quarks and gluons will be b-tagged with a mistag probability f l . 15 • b-tagged events are those for which the two Higgs candidates' jets have been both b-tagged.
Therefore, events will be given different weights according to the number of b and c quarks present in each of the two Higgs candidate jets. For instance, if the two Higgs candidates' jets each contain at least one b quark, the event is assigned a weight f 2 b =0.56. This is the same probability for signal events and for QCD background events where two b quarks end up each in a Higgs candidate jet.
We could also have considered a more optimistic scenario for the b-tagging, in which each Higgs candidate is required to have two b-tags. In particular CMS has demonstrated the ability to tag pairs of B-hadrons even for angular separations ∆R bb < 0.4 [87], which suggests that this scenario could be viable also in the highly boosted limit where the B-hadrons are within a single anti-k T jet.
Results for the b-tagging efficiencies for graviton mediated Higgs pair production at the LHC 8 TeV as a function of r M are shown in Fig. 12. As we can see, for the relevant mass range we have approximately a 15% constant signal efficiency after taking into account the resonance tagging algorithm and the b-tagging.

QCD multijet background simulation
The dominant Standard Model background to multijet final states that leads to event topologies similar to the signal is QCD jet production. We have therefore produced a large sample of QCD multijets with Pythia8 [43], starting from dijet configurations and with the shower radiation taking care of generating the higher-order jet topologies. These events include a subset with two and also four B-hadrons in the final state. The resulting hadron level events are then processed through the same analysis chain as the signal events. There are several ways in which QCD radiation can mimic the conditions for resonance tagging: for example, fake mass drops can be generated from a sufficiently symmetric splitting of a quark or gluon. Note that while Pythia8 is known to underestimate the amount of QCD multi-jet topologies by a factor up to two with respect to experimental data [4], for the accuracy requested for this feasibility study we consider this precision to be sufficient. Similar results have been obtained with the Alpgen parton level event generator [89] matched to Pythia8 using the MLM matching [90]. 16 In Fig. 13 (left) we show the QCD dijet cross section obtained from the Pythia8 multijet sample at LHC 14 TeV. The dijet cross section is defined, for each mass point M , as the number of QCD events that survive the basic selection cuts Eq. In addition, we demand that the two leading jets are separated in rapidity by less than ∆y max . Note that the dijet cross sections flattens at small masses because there the selection cuts Eq. (40) have a sizable effect.
In order to achieve an efficient QCD multijet event generation, for any candidate resonance mass M we have generated dijet events with a generation cut of p T ≥ M/5, and no generation cut in rapidity. To motivate this choice, let us recall that the kinematics of massless jet pair production determine that the mass of the dijet will be given in terms of the p T of the jets and their rapidity separation in the laboratory frame ∆y by M = 2p T cosh (∆y/2). Therefore, to properly cover all phase space the generation cut for QCD dijets should be at least for any candidate resonance mass M . For the four jet configuration, it is reasonable to require a minimum p T value of half of that of above. Since we are using a rapidity cut of ∆y max = 1.3, we find that the Pythia8 minimum p T in generation should be p min T ∼ M/5. We have explicitly verified that the QCD dijet cross-section is not modified if looser generation cuts are adopted.
We also show in Fig. 13 (right) the background rejection factors, defined as the fraction of the QCD dijet events which are mistagged as arising from a heavy resonance, both with and without b-tagging. Note that the background rejection probability is approximately scale invariant: similar mistag probabilities are obtained for all values of the mass. It is clear that the QCD background cross sections is reduced by a combination of the resonance tagging and b-tagging by several orders of magnitude. For example, with b-tagging the mistag probability is about 10 −4 , constant to very good approximation in all the relevant mass range. This improvement is due to the requirement that each Higgs candidate should be associated with two identified b-hadrons, a topology that is less frequent in QCD multijets. The decomposition of the mistag probabilities according to the number of boosted object tags in shown in Fig. 14, and is qualitatively similar to what was seen for the signal in Fig. 4.
In summary, our study of the QCD background rejection confirm the consistency of the resonance tagging algorithm, since it makes possible to simultaneously explore the low mass and high mass region, achieving similar signal efficiencies and background rejection factors in all the mass range.

Model independent exclusion limits
We will now combine the results of the signal efficiencies and the multijet background rejection of the resonance tagging algorithm to derive model independent bounds on BSM scenarios with enhanced Higgs pair production in the 4b final state. This information is enough to derive the values of the cross section times branching fraction σ (pp → X) BR (X → HH) that can be excluded at the 95% confidence level from a measurement of the QCD b-tagged multijet cross sections, as a function of the mass of this hypothetical resonance. In the following, to compute the number of signal and background events, we will assume a total integrated luminosity of L = 25 fb −1 at 8 TeV and of 500 fb −1 at 14 TeV.
For each candidate resonance mass, M , we compute the number of background events in a mass window of width f m =15% around M . The local p-value for each mass point M based on the expected number of signal and background events, N s and N b respectively, in the mass window considered, is given by where Erf is the error function, and one assumes that the number of background events in each mass bin is a Poisson distribution with mean N b . 17 Then requiring the condition p = 0.05 determines the number of signal events N s that would allow an exclusion of the background-only hypothesis at the 95% confidence level, namely Using Eq. (44) to determine the value of N s in a given mass window, we can obtain the model independent bound on the combination σ (pp → X) BR (X → HH) by correcting the number of events for the signal tagging efficiency, with and without b-tagging, the Higgs to bb branching fraction and the assumed total integrated luminosity L, where the signal efficiency SignalEff(M ) is derived from the MadGraph radion and graviton samples, see Fig. 12. We have used BR H → bb = 0.577 from the Higgs Cross Section Working Group recommendations [95]. The 95% excluded model-independent cross sections times branching fractions are shown in Fig. 15. We see that that we are sensitive to cross sections as small as 200 fb 17 Eq. (42) is only valid where both Ns and N b are much larger than one, in the opposite case one has to use the corresponding discrete Poisson formula for the p-value, which involves the incomplete Gamma function. (50 fb) at M ∼ 500 GeV at LHC 8 TeV (14 TeV), while at higher masses, M ∼ 2 TeV, the 4b final state is sensitive to cross sections as small as 1 fb at both energies. Note that the increase in luminosity when going from 8 to 14 TeV is partially canceled by the corresponding increase of the high mass QCD dijet cross sections. On the other hand, signal cross sections in relevant models also increase when going from 8 to 14 TeV, so all in all we obtain a substantial improvement in exclusion reach when increasing the center of mass energy.
These results confirm that the 4b final state can be relevant for many new physics scenarios that lead to enhanced cross sections for resonant Higgs pair production, and the search strategy that we propose makes it possible to efficiently explore a wide range of resonance masses within a common analysis.

Graviton and radion searches in the 2H → 4b channel
Now we consider the specific benchmark scenarios for radion and graviton production introduced in Sect. 3. First of all, we summarize in Table 2 the model parameters that we adopt here. The mass scales and branching fractions to Higgs boson pairs are kept fixed, and only the couplings of the radion and graviton to gluons will be varied. For the graviton, we consider two different values of the coupling, c g = 1 (G-Brane) which corresponds to the RS1 model and c g = 1/35 (G-Bulk) as in the bulk models.
For the radion we study the nominal coupling κ φ g (R-Bulk), as well as the case in which this coupling is enhanced by a factor of ten by some unspecified mechanism, such as when the radion arises as a composite bound state (R-Comp). We use the same mass scale in the two cases, Λ φ = Λ G = 2 TeV: although the two scales are related, we prefer to explore independently the radion and graviton scenarios. Let us recall that the cross section scale as 1/Λ 2 φ,G , so any different choice of the mass scale Λ φ,G will lead to a trivial rescaling of the cross section.
First of all we evaluate the expected number of events for these four benchmark points at LHC 8 and 14 TeV, using the results of Sect. 3. We take into account the branching fraction of the Higgs bosons into bb pairs. As before, we assume total integrated lumi-radion Production Scenario  Table 2: Parameters of the benchmark scenarios for radion and graviton production. For the radion we consider both the nominal value of κ φ g (denoted by R-Bulk), and a coupling ten times larger that could arise for example in composite dual scenarios (denoted by R-Comp). For the graviton we consider two different values of the gluon-gluon-graviton coupling, c g = 1 as in RS1 (denoted by G-Brane settings) and c g = 0.02 as in bulk models (denoted by G-Bulk). nosities of L = 25 fb −1 at 8 TeV and of 500 fb −1 at 14 TeV. The model cross sections for the benchmark scenarios can be easily obtained from the results of Sect. 3, in particular from Fig. 9. The number of expected events is shown in Fig. 16, after accounting for the selection efficiencies from resonance and b-tagging analysis. At 8 TeV we expect just a handful of events at low masses for the R-Bulk and G-Bulk points, and about one thousand events (few tens of events) and low (high) masses for the R-Comp and G-Brane points. At LHC 14 TeV on the other hand we have a large enough number of events for all masses and all benchmark points, thanks to both the increased resonance production cross sections and the higher integrated luminosity. An illustration of the type of signal that one might observe is given in Fig. 17 for the case of an RS1 graviton with a mass of 1 TeV in 14 TeV pp collisions. Note that at this mass, the signal involves 0, 1 and 2-tag categories combined. Now in Fig. 18 we show the same excluded cross sections at the 95% confidence level of Fig. 15   compared to the specific cross sections of the four different model scenarios of Table 2. LHC 8 TeV we can explore a large part of the parameter range of the graviton models up to 2 TeV, but the default radion scenario seems to be out of reach, unless its cross section is enhanced by some mechanism, for instance as in the composite duals discussed in Sect. 3. At the LHC 14 TeV on the hand we are sensitive to R-Bulk scenario with κ φ g = 1, for most of the mass range up to M φ = 2 TeV. Likewise, we could exclude a Bulk graviton up to masses of 2.5 TeV. Therefore, after the energy increase to 14 TeV most of the parameter space of the radion and massive KK graviton models will become accessible in the 4b final state.
Using these results, it is also possible to determine the 95% confidence level exclusion ranges for some of the parameters of the benchmark scenarios. We can keep all the parameters as in Table 2 and determine the exclusion ranges for the couplings of the gluons to the massive KK graviton c g and to the radion κ φ g and scan the allowed values for Λ φ .
We show the results in Fig. 19. In the case of the graviton coupling, we see that at 14 intermediate mass regions, the exclusions are based on a large number of events, corresponding to small signal over background ratios. This can be understood from Eq. (44), which tells us that at the exclusion limit, Ns Graviton Mass ( TeV )  TeV the 2H → 4b final state can access essentially all the relevant range, from the RS1 value c g = 1 down to the bulk value of c g = 1/35. For the case of the radion coupling κ φ g , we see that at 8 TeV we are sensitive to values down to κ φ g 0.06 around 750 GeV, while at 14 TeV the LHC can exclude a bulk radion (with the default value for the coupling) for masses between 300 GeV and 2 TeV at least.
In summary, we have shown in this section that the 2H → 4b final state offers a promising channel to probe enhanced Higgs pair production at the LHC, despite the overwhelming QCD multijet background. The combination of jet substructure techniques and b-tagging makes it possible to probe a wide region in the parameter space of various benchmark models. Therefore, we advocate that the experiments explore this new channel in order to complement existing searches of new heavy resonances in other, more traditional, channels.

Conclusions and outlook
In this paper we have presented a new strategy for heavy-resonance searches in multijet final states, which attempts to unify in a single approach the techniques used in the boosted and resolved regimes. By classifying events as a function of the number of massdrop tags, we can smoothly interpolate between the boosted regime, where jet substructure techniques can be used, and the resolved regime, where the final state particles appear as well separated jets. In particular, we have considered the process X → Y Y → 4z in which the resonances Y are pair produced from the decay of a heavier resonance X and then decay into a pair of QCD partons, then observed as jets. We have shown that our strategy leads to approximately scale-invariant signal selection efficiencies and background rejection rates.
As a benchmark scenario, we have considered Higgs pair production in extra dimension models, where the Higgs pair is produced from the decay of a heavy graviton or radion, and then decays into four b quarks. Note however that the kinematical structure of the final state of the benchmark model holds for other scenarios, such as composite models, with more freedom on the couplings and therefore in the cross sections strengths. By comparing with the QCD multijet background, we have derived the model independent 95% confidence level exclusion ranges for the cross sections for σ (pp → X → HH) where the Higgs bosons decay into a bb pair, and showed that a substantial region of the parameter space of these models can be successfully explored in this final state with the tagging strategy that has been proposed.
In the particular case of graviton and radion production, it would be especially interesting to study the feasibility of radion/graviton searches in the bbγγ decay channel [96,97]. This is a cleaner channel than the all hadronically decaying case, since the two high p T photons substantially reduce the QCD background. In this final state, by varying r M one also moves from the boosted regime (with one single fat jet in the final state) to the resolved limit, with two well separated jets in the final state. Another interesting final state to search for enhanced Higgs pair production would be bbZZ.
The approach advocated in this paper could also be applied to other relevant problems, for example top quark pair production, again providing a smooth coverage across the transition between the resolved regime, relevant for SM measurements, and the boosted regime, where substructure techniques [98,99] are used to enhance the potential of new physics searches.