Pulling out all the stops: searching for RPV SUSY with stop-jets

If the lighter stop eigenstate decays directly to two jets via baryonic R-parity violation, it could have escaped existing LHC and Tevatron searches in four-jet events, even for masses as small as 100 GeV. In order to recover sensitivity in the face of increasingly harsh trigger requirements at the LHC, we propose a search for stop pairs in the highly-boosted regime, using the approaches of jet substructure. We demonstrate that the four-jet triggers can be completely bypassed by using inclusive jet-HT triggers, and that the resulting QCD continuum background can be processed by substructure methods into a featureless spectrum suitable for a data-driven bump-hunt down to 100 GeV. We estimate that the LHC 8 TeV run is sensitive to 100 GeV stops with decays of any flavor at better than 5σ-level, and could place exclusions up to 300 GeV or higher. Assuming Minimal Flavor Violation and running a b-tagged analysis, exclusion reach may extend up to nearly 400 GeV. Longer-term, the 14 TeV LHC at 300 fb−1 could extend these mass limits by a factor of two, while continuing to improve sensitivity in the 100 GeV region.


I. INTRODUCTION
Well over one hundred null searches for supersymmetric signals at the LHC have now been completed, leaving us to ponder the fate of naturalness after the Higgs boson discovery. As we await the energy upgrade, much attention is being focused on the possibility that supersymmetric particles have been produced at the LHC, but for one reason or another are buried amidst the copious Standard Model backgrounds. In this paper, we turn to one of the simplest of these possibilities: the lightest superpartner is a stop eigenstate, and decays promptly to two quarks through a baryonic R-parity violating (RPV) coupling.
Such a situation is in fact quite well-motivated within the context of "effective" or "natural" supersymmetry [1][2][3][4][5], where only the superparticles required to regulate the Higgs mass need to be within the immediate reach of the LHC, and the detailed spectrum and dynamics above the multi-TeV scale can be left unspecified. In this framework, third generation squarks can be amongst the lightest superparticles, forcing us to take seriously the possibility that a stop eigenstate sits at or near the bottom of the SUSY spectrum. R-parity conservation is still often assumed in natural SUSY, but its violation must also be considered. The main motivations for R-parity are stabilization of a supersymmetric dark matter candidate and suppression of proton decay. The dark matter motivation is attractive, but the issues of dark matter and electroweak naturalness might easily be decoupled. The proton decay motivation does not require R-parity per se, but only that baryon-number violation (BNV) and lepton-number violation (LNV) are not simultaneously active. Moreover, the motivation for R-parity is further weakened by the presence of R-parity-conserving operators at dimension-5 with simultaneous BNV and LNV. In a relatively model-agnostic approach with a multi-TeV cutoff, these operators in any case force us to consider proton-stabilizing symmetries other than R-parity. RPV is also sometimes considered dangerous given precision tests, in particular K −K and n −n mixing in the BNV case. However, couplings that are small enough to avoid these constraints are still usually large enough to allow for prompt RPV decays from the perspective of the LHC [6] (see [7] for a complete list of constraints).
The idea that SUSY might be both "natural" and violate R-parity has received increasing attention in the past few years [4,[12][13][14][15][16][17][18][19][20][21][22][23], exactly because the relatively small production cross sections and non-canonical decay topologies lend themselves to evading LHC searches. Many spectra and choices of active RPV operators exist, and the nontrivial flavor structures of the couplings must be specified. In some cases, especially those involving leptonic RPV, exclusion or discovery prospects are quite good (see [12] for proposals and recast limits, and e.g. [24,25] for some recent experimental searches). Successful baryonic RPV searches have also been performed, such as LSP gluinos decaying fully hadronically into multiple jets [26][27][28][29]. However, an LSP stop decaying promptly to two jets via baryonic RPV has proven particularly evasive. Stop pairs are produced at much smaller rates than gluino pairs of comparable mass, and the lower-multiplicity 4-jet final state is even more difficult  [8][9][10][11] to account for stop acceptances relative to coloron or hyperpion acceptances.
to disentangle from the pure QCD backgrounds. Another major complicating aspect at the LHC is the multijet triggers, which can heavily prescale-away the signatures of stops lighter than several hundred GeV. Some of the best current direct limits actually come from LEP, which rules out mt < ∼ 90 GeV [30]. A recent search at the Tevatron extends this limit up to only about 100 GeV [31]. However, so far, direct searches for pair-production of dijet resonances at the LHC have failed to reach the sensitivity necessary to place constraints for any stop mass [8][9][10][11]. A snapshot of the current situation can be seen in Fig. 1. In fact, the inevitable rise of trigger thresholds with instantaneous luminosity and beam energy leaves us to wonder whether the LHC will ever be sensitive to this signal. At the very least, this trend suggests that masses near the current limit of 100 GeV might be left unexplored. 1 One way around these difficulties is to search for the stop as a dijet resonance produced in the decays of heavier colored superparticles, such as gluinos [33] or sbottoms [6] (or possibly the heavier stop eigenstate), or to simply set bounds using the associated leptonic activity and high H T of these decays [34][35][36][37]. Naturalness suggests that these colored superparticles should also not be far above 1 TeV, and might be produced with observable rates. It is also possible to invoke Minimal Flavor Violation (MFV), which suggests that stops dominantly decay (with a branching ratio≃ 95%) intobs orbd [13]. It was pointed out in [38] that incorporating b-tagging into the triggering might allow the direct stop pair signal to write to tape with higher efficiency, and subsequent kinematic analysis can discriminate it from flavored backgrounds. The use of b-tagging is also usually considered to help identify stop production in gluino decay. However, at present, there seems to be no strategy at the LHC that guarantees full coverage of the generaltt * → 4j signal, exploiting neither spectrumdependent nor flavor-dependent features of the theory.
Rather than consider the generic signal lost, we here exploit what is now a well-worn trick [39][40][41]: we focus on "boosted" production at high transverse momentum, and apply the methods of jet substructure to "stop-jets," each containing an entire stop decay. Indeed, an analysis of this type has already been applied in the search for light RPV gluinos [27], and has been studied for color-octet scalars [42]. The advantages of a boosted/substructure analysis are manyfold, including better S/B, automatic resolution of combinatoric ambiguities, and improved mass resolution due to more complete decay radiation containment and rejection of uncorrelated soft radiation. But probably the primary advantage for a stop pair search is in how a substructure treatment can process the continuum QCD background. The standard resolved 4-jet analyses partition the leading jets into two pairs so as to form two stop candidates, either minimizing the candidate mass asymmetry (CMS) or minimizing a ∆R measure for the two decays (ATLAS). The stop is searched for as a bump in the average pair-mass spectrum. In either case, the requirement of four well-separated and high-p T jets passing the triggers, as well as other downstream analysis cuts, leads to a highly-shaped background QCD spectrum with a pronounced "trigger turn-on" peaking at 100-200 GeV. As a consequence, the most recent CMS [9] (ATLAS [11]) analysis of 2011 data does not even search for masses below 250 GeV (150 GeV).
By contrast, as we will see, a carefully constructed jet substructure approach potentially leads to a largely featureless average-mass spectrum for the QCD background from a few 10's of GeV up. This is because jet substructure, unlike traditional jet analyses, gives us far more flexibility in assigning sprays of hadrons to individual hard "quarks" or "gluons." In particular, we can make these assignments in a way much closer to QCD itself, which is approximately a scale-invariant theory. To best exploit this feature, we will also capitalize on the summed jet-H T trigger rather than triggers that count specific numbers of jets above some threshold. The jet-H T trigger, which exists in some form in both CMS and ATLAS, probably has the least sensitivity to precisely how the event activity groups into standard LHC jets.
With an unbiased background spectrum, a bump-hunt becomes much more tractable, regardless of what stop mass we consider. Either a parametrized spectral fit or kinematic sideband approaches like the ABCD method can be employed to determine the QCD background, and we explore the viability of several such approaches, including some novel ones. We find that 100 GeV stops should be visible in the 2012 data set with better than 5σ statistical significance. Exclusion reach should extend up to more than 300 GeV. We also consider what might be possible in the MFV case, using a b-tagged version of the analysis. S/B improves dramatically, for example yielding 10σ-level statistical sensitivity to 100 GeV stops, and exclusion up to almost 400 GeV. For the longer-term Run II of the LHC, with a projected 300 fb −1 at 14 TeV, we estimate untagged exclusions extending up to 650 GeV, with discovery-level sensitivity or better between 100 GeV and 500 GeV.
Our paper is organized as follows. In the next section, we discuss our substructure procedures, and explore the response of the signal and QCD continuum in simulation data. We present our data analysis techniques and final sensitivity estimates in Section III. We conclude in Section IV. Two appendices contain more details of our simulations and more in-depth substructure studies.

II. JET SUBSTRUCTURE TECHNIQUES AND ANALYSIS CUTS
The fundamental obstacle to any multijet search at the LHC is the overwhelming production rate for continuum QCD backgrounds. The dynamics of QCD is approximately scale-invariant at high energy. Convolved with smoothly-falling parton distribution functions, a good stage is set for extracting signals with sharply-localized mass features, such as RPV stop pairs. However, the standard approaches to reconstructing QCD events are far from scale-invariant. "Jets" as usually defined are tied to a dimensionful p T threshold, as well as a dimensionless jet radius R. For example, in a multijet search with p T (j) > 100 GeV and R = 0.5, there is an absolute minimum invariant mass between jet pairs of approximately 50 GeV. Any QCD background invariant mass spectrum constructed from such jets will start at 50 GeV, rise at higher invariant masses as the efficiency turns on, and only then turn over into a smoothly-falling shape. Searches for low-S/B features in the broad turn-on region or at the peak are often not considered, as the precise signal and background shapes in these regions have high sensitivity to reconstruction uncertainties.
Indeed, these kinds of considerations have limited the range of applicability of current LHC stop searches based on the multijet approach, as seen in Fig. 1. To control triggering rates as the LHC continues to increase instantaneous luminosity and energy, the jet p T thresholds are gradually increasing. For example, to keep 4-jet triggering rates at their 2012 levels during the projected 300 fb −1 run at 14 TeV, offline p T (j) thresholds would need to be increased to 150-200 GeV. Stops sitting near the current exclusion limit of m = 100 GeV might therefore never be visible. Even with 2011 data, only the very low-luminosity ATLAS search [10] was capable of probing masses near that limit, but was not sensitive to stop pair cross sections. Despite the fact that the cross section is becoming quite large near 100 GeV, basic triggering and analysis cuts at nominal LHC luminosities are carving away the signal.
As has been pointed out in [43,44], reconstruction-induced biases can be highly ameliorated by applying modern approaches to jet substructure on "fat-jets." Jet substructure can be less tied to fixed p T or ∆R thresholds for the individual reconstructed hard partons, and can put more focus on dimensionless quantities such as ratios of masses or p T 's. The "price" for applying a substructure-based search is that we must work in the boosted region of production phase space, which might still represent only a small fraction of the total signal. Nonetheless, as we will see, the tradeoff is more than worthwhile. In any case, existing 4-jet searches are already capitalizing on the boosted or semi-boosted region of phase space, which offers the additional advantages of better S/B and much-reduced combinatoric ambiguities.
We are in the midst of an ongoing boom of jet substructure ideas, many of which can have applicability to the boosted stop pair signal. We here focus on the relative-p T declustering method used in the JHU top-tagger [45] and the diboson-jet tagger of [46], which is itself based on the "BDRS" method of [41]. These methods identify localized clusters of energy within the fat-jets as "subjets," and subsequently treat these objects similar to ordinary QCD jets. However, before proceeding, we emphasize that ideas such as (but not limited to) pruning [47], N-subjettiness [48], template overlaps [49], shower deconstruction [50], and Q-jets [51] could all worth more dedicated study in this context, and might lead to further improvements.
For our fat-jet clustering radius, we choose a very large value of R = 1.5, roughly giving us hemisphere-sized jets. The motivation for such a large radius is twofold. First, the large catchment area gives us a broader reach in stop masses. For given momentum thresholds, heavier stops will be able to pass with less boost, and therefore will have more widelyseparated decay products. Heavier stops also have much smaller cross section, and the very highly-boosted portion of phase space may be too poorly-populated to be exploited. In fact, for much of the mass range probed in this paper, the quarks in the stop decays are well-separated enough to be reconstructed as individual jets of normal radius. However, we find it interesting that these masses can alternatively be covered using substructure procedures, raising the question of whether a traditional analysis is even necessary. The second motivation for the large radius is more subtle, and has to do with how the continuum QCD background is processed. Interplaying with momentum and substructure cuts, the jet radius sets a maximum mass (analogous to the minimum mass for normal jet-pair masses described above). Larger radii push the turn-off of the mass spectrum out to higher values, creating a less steeply-falling background. Larger radii also somewhat improve S/B. This is because enlarging the jet radius increases sensitivity to wider-angle radiation, and provides a primitive form of color discrimination. At high momenta, stop pair production is dominated by qq →tt * , whereas the QCD background contains more highly-colored processes such as qg → qggg. (For a comparison of different R's, see Appendix B.) Fixing R = 1.5, we cluster jets in each event using the Cambridge/Aachen (C/A) algorithm [52,53] as implemented in FastJet3 [54], and study the leading two jets within |η| < 2.5. Then we iteratively undo the clustering stages individually within each jet. Now viewed in reverse as a splitting, at each stage we can look at the two branches "a" and "b" and decide whether the splitting should be considered "hard" or "soft" according to some prescription. For example, a soft splitting might resolve a low-p T collection of radi-ation at the jet's periphery. In our specific implementation, derived directly from the first declustering stage of the diboson-jet-tagger of Ref. [46], a splitting is considered hard if both a and b carry appreciable p T relative to the original fat-jet, and if their individual m/p T ratios are not too large. Explicitly, min Otherwise the splitting is considered soft, in which case the lower-p T branch is thrown away, and the declustering is continued along the surviving branch. The procedure is repeated until a hard splitting is encountered, or no more jet constituents remain (in which case the entire event is vetoed). The two branches at the hard splitting are our two subjets, representing our assignment of the jet's radiation to the two quarks in the stop decay. 2 The reconstructed stop candidate is the four-vector-sum of these two subjets.
Most of our ability to discriminate stop pair production from ordinary QCD stems from the fact that the stop events contain two subjet-pair resonances of equal mass. Therefore, instead of looking for a bump within the distribution of the individual stop candidate masses m 1 and m 2 , we look for a bump in the joint distribution of (m 1 , m 2 ). Practically, we can turn this into a 1D bump-hunt by first focusing on the region of small mass asymmetry A ≡ |m 1 − m 2 |/(m 1 + m 2 ), and then constructing the spectrum of the averaged mass m avg ≡ (m 1 + m 2 )/2. We pick a nominal mass asymmetry threshold of 10%.
We have found that S/B can be further purified by a handful of additional cuts. In particular, we can place a cut on the stop-pair CM-frame production angle (as is done in [10,11]), formed by actively boosting the entire 4-subjet system to rest in the lab frame and measuring the angle of either stop candidate with respect to the beamline. We call this angle θ * , and place a cut | cos θ * | < 0.3. Finally, we exploit the fact that energy splittings within QCD jets tend to be very asymmetric, whereas in stop decays they tend to be more democratic. Within each stop pair, we demand that the p T 's of subjets a and b relative to one another satisfy min[p T (a), p The final ingredient needed to define an analysis is to establish a trigger for the events. We choose the total jet-H T trigger of CMS, which sums up the p T 's of all ordinary jets. 3 In Ref. [55], based on 2011 data, the H T trigger was used for a classic jets+/ E T style SUSY search. The trigger was found to be fully efficient for final reconstructed events with H T > 750 GeV, summing over R = 0.5 anti-k T jets with p T > 50 GeV. Anticipating slightly harsher triggers in 2012, we conservatively set our threshold at 900 GeV. 4 2 The original BDRS method uses a very similar procedure [41], relying instead on a mass-drop criterion and a somewhat different momentum-asymmetry criterion. While BDRS can be adapted for use in the boosted stop search, there are important caveats which we discuss in more detail in Appendix B. 3 ATLAS also has a dedicated single-fat-jet trigger, though this is defined with R = 1.0. Nonetheless, this trigger could still be used as the basis of a substructure analysis very similar to the one that we explore here. 4 A looser offline H T cut may in fact be feasible with 2012 data, and would only improve our results. We thank Keith Ulmer for bringing this point to our attention.
We test our substructure methods and cuts on signal and QCD background simulation samples for the 8 TeV LHC. Full details of their generation can be found in Appendix A. We note here that the signal samples are matched up to one jet emission at the production stage, which gives a more accurate H T spectrum. The QCD continuum simulations are matched up to four hard partons (including b-quarks), using the CKKW-L [56] prescription implemented in Pythia8 [57]. We take a quark-flavor-conscious Durham k T distance as our merging measure, and a merging threshold of 50 GeV. Other background samples, such as tt, are generated in MadGraph5 interfaced with Pythia8, or self-contained within Pythia8. All samples are processed through a simple, perfect 0.1×0.1 grid "calorimeter" in η-φ space. To prove the robustness of our methods in a high pileup environment, we also introduce pileup events into the simulations and then apply a form of event-wide trimming [58] tailored to pileup removal before the fat-jet clustering stage. (This procedure also heavily reduces the impact of the underlying event and soft ISR.) After declustering the fat-jets, individual subjets are energy-smeared similar to normal LHC jets (see Appendix A for further details). Figure 2 shows the effect of the cut-flow on the m avg spectra of continuum QCD and some example signal mass points. There are several features that are worth noting. First, the ubiquitous turn-on peak is present, but has been pushed down to the O(10 GeV) scale, far away from our signals. This can easily be understood from the interplay between our H T cut, the 10% relative-p T requirement in the declustering, and the size of our calorimeter cells. (Note that the lowest-p T subjets that we work with are roughly 40 GeV.) For the signal, before the A < 0.1 cut, it is possible to see three distinct mass features: the turn-on, the true stop-mass peak, and an intermediate peak near half of the stop mass. The last feature arises from events where one stop is correctly reconstructed, but one of the quarks was lost for the other stop. After the A cut, the stop peaks become very clear. The additional cuts both reduce the overall size of the background and further tighten the signal peaks. Throughout the entire set of cuts, the QCD background in the vicinity of the signals stays "featureless." The final signal efficiencies relative to the inclusive pair production rates are 5 × 10 −5 for 100 GeV and 4 × 10 −3 for 300 GeV. For comparison, stop pair acceptances for the standard 4-jet searches are usually at the 10 −3 -10 −2 level.
We show how the signal peaks appear on top of the continuum QCD background in linear scale in the left panel of Fig. 3, and give an indication of S/B relative to the QCD background's bin-by-bin statistical errors. It is clear already from this simple plot that the 100 GeV stop should be visible with high significance, and that exclusion reach should extend beyond 200 GeV. We also show on this plot the largest subdominant background contributions, namely tt and W +jets. (Other backgrounds such as diboson and singletop can be shown to be far smaller.) We will not directly address the impact of these backgrounds on a search, but they are clearly important. This is especially true for tt, which contributes a broad bump peaked near 175 GeV, and at a size only O(1) smaller than our stop signal. 5 However, the multibody structure of this background is under much better 5 The fact that tt is not a larger contribution is perhaps somewhat surprising, given that for mt ≃ m t , the inclusive tt cross section is about six times larger thantt * . About half of this factor comes from the tt all-hadronic branching fraction, since only all-hadronic events are efficient at passing the H T cut and subsequent substructure cuts. It is also important to realize that for high-p T central production, the difference in cross sections is not as big. (Asymptotically, the factor of six reduces to a factor of two.) Finally, the large fraction of partial reconstructions with two-body substructure significantly broadens the theoretical control than pure QCD, and its normalization could be extracted in the highly orthogonal semileptonic channel. We therefore anticipate that it could be systematically subtracted or accounted for in a constrained fit. Indeed, it can even serve as a useful calibration peak. If it is necessary to further suppress tt, it might be possible to do so with supplementary substructure cuts that can pick out and reject 3-body features, without highly rescultpting the continuum QCD. (E.g., N-subjettiness [48] observables or the dimensionless variables of the HEPTopTagger [59] would be appropriate to study.) Regardless, some degradation of sensitivity in the vicinity of m t should be expected in reality.
If the RPV coupling obeys MFV, then almost every stop decay will contain a b-quark. It therefore becomes possible to exploit a b-tagged analysis. We show in the right panel of Fig. 3 the m avg spectra after demanding that at least one of the four subjets is tagged, assuming flat (b, c, q/g) tag rates of (60%, 10%, 2%). The S/B (and S/ √ B) improves top peak shape.
dramatically, as does the relative contribution of tt to the background budget. Exclusion reach to above 300 GeV already appears highly likely. These distributions set the stage for our data analysis in the next section.

III. SEARCH STRATEGIES AND SENSITIVITY ESTIMATES
Extraction of a bump on top of a smoothly-falling background spectrum is a classic problem that has appeared many times already at the LHC. Probably the most famous recent application is the observation of the Higgs resonance in the continuum diphoton spectrum [60,61], but this strategy has also been applied by CMS in its paired dijet resonance searches [8,9]. Relying solely on the assumption that the background is "featureless," the observed spectrum can be fit to a parametrized function with or without the addition of a signal bump. The parametrization used by CMS, which we also use here, is a four-parameter function of the form [9] dP dm avg where p 0 , p 1 , p 2 , p 3 are free parameters and √ s is the proton-proton center-of-mass energy (8 TeV for 2012). 6 There are also many other ways to directly estimate the QCD m avg spectrum from the data, using control regions. The common ABCD method was used by ATLAS in its own paired dijet resonance searches [10,11]. This method requires defining sideband cuts in two variables. The nominal cuts define region A, and the other three choices of 2D cuts (nomi-nal+sideband, sideband+nominal, and sideband+sideband) define three signal-depleted regions B, C, and D. Assuming small correlations between the two variables (an assumption that must be justified in Monte Carlo studies), the region-A background spectrum can be derived by taking the bin-by-bin ratios of counts B×C/D. For its ABCD-based search, AT-LAS uses the variables A and | cos θ * | to define its four regions. We run our own version of this search, taking sideband cuts A = [0.1, 0.4] and | cos θ * | = [0.3, 0.8].
In addition to the shape-fit and ABCD methods, we also explore two supplemental techniques which may further improve statistical power.
In the first of these, we run a simultaneous shape fit with Eq. (1) over our nominal signal region and a nearby asymmetry sideband region A = [0.1, 0.2]. This increases the data statistics available to the fit by O(1), as well as folding in more discriminating characteristics between the stop signal and QCD. The method is based on the assumption that the QCD spectrum changes very slowly as a function of A when A ≪ 0, whereas the signal is strongly 6 It may be possible to develop parametrizations that more directly incorporate the analytic structure of QCD (see [62,63]), though we have not explored this possibility. peaked in A for m avg ≃ mt. We observe exactly this behavior in our simulations. To implement the simultaneous fit, we set the QCD fit parameters identical between the two A regions. We can then add into the fit the signal shapes and relative normalizations appropriate to each individual region. This simple two-region fit might also be improved by more finely subdividing in A, and thereby more fully exploiting the nontrivial shape of the signal over this variable. Equivalently, this can be viewed as a full 2D fit over the small-A region of the (m 1 , m 2 ) plane.
The final method that we study (inspired in part by [64]) is based on the assumption that the two fat-jet masses can be considered approximately uncorrelated in background events. If this assumption holds, then the m avg spectrum in the A → 0 limit can be directly predicted from the spectra of individual fat-jets. Given a single-jet probability distribution dP/dm jet (where m jet represents the subjet-pair mass), we get The major advantage of this method is that the statistics available for measuring dP/dm jet are enormous, since no tight cuts should be placed on the event-by-event jet-mass asymmetry. However, the uncorrelated assumption needs to be carefully studied in simulation. (It may also be tested to some extent in data, by inverting some of our final cuts.) To keep the kinematics similar to our final signal region, we measure the dP/dm jet spectrum in a control region with only the H T and | cos θ * | cuts in place. We pick a random jet amongst the leading two, and further demand the min[p T (a), p T (b)]/ max[p T (a), p T (b)] cut for only that jet. (The signal contamination in this single-jet control sample is negligibly small.) Once the dP/dm avg spectrum estimate is derived, we use it as a template for a one-parameter normalization fit to our signal-region spectrum, with or without a signal bump added. In practice, we approximate dP/dm jet with a very finely-binned spectrum measured from our "data," apply the above transformation, and then integrate back to our nominal m avg binning (10 GeV). Since our four search methods (shape fit, ABCD, asymmetry-sideband, and single-jet template) all rely on assumptions about the behavior of the QCD background, we can first get some sense of the validity of those assumptions. To do so, we compare the signalregion QCD spectrum to the predictions of the four methods, without signal. Here and below, we focus on the mass range m avg = [60, 500] GeV. Also, due to the finite statistics of our QCD Monte Carlo sample, which is actually comparable to a 2012-like data set, we smooth out the spectra in the signal and control regions by fitting to the exponential of a fifth-degree polynomial. These fits generally have high goodness-of-fit probability. 7 (We do not apply smoothing to the single-jet spectrum, which has much higher statistics.) We show the results of the comparison in Fig. 4, including the original Monte Carlo data with its statistical errors. For m avg < ∼ 300 GeV, the agreement is generally better than 5% (10% with the smaller-statistics b-tagged sample). At higher masses, more pronounced disagreements develop, but their significance is likely not large given the growing error bars. This is also entering the m avg region where the subjets are nearly separated by ∆R = 1.5 and the spectrum is turning off, in which case we might not be surprised to find more sensitivity to the different control region cuts used in the estimation procedures. Needless to say, a higher-statistics simulation would be useful here to better gauge the level of agreement at both higher and lower masses. Still, given the encouraging agreement over the mass range in which we will shortly find our best sensitivity, we can proceed with our analyses without that the latter does indeed furnish a high-probability fit. However, we have explicitly chosen different functions for smoothing the raw Monte Carlo and for fitting the derived pseudodata, in order to help prevent spuriously good pseudodata fits. fear that their underlying assumptions are grossly invalid, at least for matched QCD. 8 We now apply the four data-driven methods to search for the stop signal in our simulations, using the matched QCD background. For all four methods, we search for a signal using χ 2 differences as a discriminator. In the case of the fits, we take the χ 2 difference between background-only and signal+background fits, using √ N error bars for the observed bin counts. For the ABCD method, which directly predicts the background spectrum with no free parameters, we construct χ 2 bin-by-bin by combining in quadrature the statistical errors of the observed spectrum and the predicted spectrum. (The latter is itself derived by simple propagation-of-errors from the B, C, and D bin counts.) In all cases, the signal strength and shape are fixed, and systematic errors are not assessed. Our results should therefore be indicative of what can be accomplished in the limit of small systematics. 9 We use our (smoothed) background simulations and signal+background simulations as the basis of a large number of pseudoexperiments corresponding to 20 fb −1 at 8 TeV. We run our various search strategies on these pseudoexperiments to build up ∆χ 2 distributions. We find these distributions to be highly Gaussian, at least out to about three standard deviations. To better parametrize their separation, we therefore fit the ∆χ 2 distributions to Gaussians. From these we can derive a "median discovery significance" and a "median exclusion significance." These are the distance between the medians of the two ∆χ 2 distributions, as measured in units of the background-only pseudoexperiments' σ (for discovery) or the signal+background pseudoexperiments' σ (for exclusion). We can set a benchmark for discovery at 5σ B separation, not accounting for the look-elsewhere effect. We can set a benchmark for exclusion at 2σ S+B , which is practically equivalent to the usual 95% CL S criterion. For exclusion, we also consider fluctuations at ±1σ B about the background-only median, and recompute the exclusion level in σ S+B units.
Figs. 5 and 6 show the final results for untagged and b-tagged analyses, respectively. It is clear that all four search methods perform comparably, but that the single-jet template 8 We have also rerun the untagged comparison on a larger background sample generated wholly within Pythia8, based on showered 2 → 2 production without matching. The shape fit and ABCD methods continue to improve across the full mass range, generally agreeing with the signal region spectrum to better than 2% for m avg < 350 GeV (and still consistent within MC errors). The asymmetry-sideband agreement exhibits a -20% dip above 300 GeV, similar to the matched case, but with much higher significance. Interestingly, the single-jet template also develops a broad +20% discrepancy above 300 GeV, which appears to be more severe than for the matched case. Again, these features occur in regions where our sensitivity is in any case becoming poor, but they would need to be more carefully investigated in a fully realistic search. 9 Systematic errors affecting the signal shape and normalization will include the jet/subjet energy scales and resolution uncertainties, and probably PDF and scale errors for the high-p T cross section. Possible systematics affecting our different data-driven estimates of the QCD spectrum would need to be investigated with higher-resolution Monte Carlo simulations with a full detector mock-up and/or in control data. method tends to edge out the other three, and that the asymmetry-sideband method offers a small but consistent improvement over the simple shape fit. (In fact, for exclusion significance, the single-jet template method gives results very close to what would be inferred with a naive S/ √ B analysis with optimized mass windows.) The similarity of the results is encouraging, and suggests that experimentalists will have many alternative choices for performing cross-checks of a tentative signal, or as fall-back options if any of these datadriven methods turns out to be unreliable. From Fig. 5, which shows the untagged analysis, we see that stops less than about 175 GeV could be discovered, and stops less than about 320 GeV could be excluded. For the b-tagged analysis in Fig. 6, masses below 250 GeV are discoverable, and exclusion sensitivity extends to nearly 400 GeV. We note that this analysis was run without re-optimization of our cuts, so it might be possible to construct an even more sensitive search. It may also be possible to make even further gains by considering a double-b-tagged search.
Looking ahead, we have also run versions of these analyses on 14 TeV simulations, assuming 300 fb −1 luminosity, and for simplicity neglecting pileup. Here, we have used a summed-jet H T cut of 1600 GeV, which keeps the rate approximately the same as the 900 GeV threshold under 2012 conditions (assuming quadrupled instantaneous luminosity at 14 TeV). The 100 GeV untagged signal remains visible, with statistical significance slightly better than our 2012 estimate, though with approximately 2-3 times smaller S/B. The discoverable range expands up to about 500 GeV, and masses of 200-300 GeV would be visible at the 10σ-level. Exclusion should extend up to 650 GeV. This last finding is comparable to that of the recent Snowmass 2013 report [32], which uses traditional jet reconstruction methods and a highly approximate background estimate. However, that search assumes 2012-like jet p T cuts, and even then is limited to the mass range above 300 GeV. By contrast, in our jet substructure version of the search there is practically no low-mass cutoff on the search range, with masses from 100 GeV to O(TeV) covered by a single analysis strategy.

IV. CONCLUSIONS
In this paper, we have addressed what has been believed to be one of the most difficult supersymmetry signatures at hadron colliders, and demonstrated that it may nonetheless be made highly visible using the tools of jet substructure. Besides serving as a crucial supplement to the LHC's broad-based program for testing naturalness, this result, if reproducible in a realistic analysis on actual LHC data, will serve as a benchmark for fully jetty searches. The implications extend well beyond just RPV supersymmetry. Thus far, multijet searches at the LHC have successfully placed constraints on the pair-production a variety of colored objects, whose prompt decays contain neither leptons, neutrinos, nor other invisible particles [8][9][10][11][26][27][28][29]. However, compared to the LSP stop with baryonic RPV, which is a color-triplet scalar undergoing a two-body decay, these searches have always relied on more color, more spin, more flavor, and/or more final-state partons. I.e., they can exploit much higher cross sections and/or more complicated event topologies. With the search that we are proposing, we have finally "hit bottom" on colored particle pair-production and decay. Any model with a color-triplet "diquark" scalar can be searched for, whether it is connected to naturalness or not, and independently of the flavor structure of its decay (see e.g. [65][66][67][68]). 10 In terms of concrete performance, we have found that jet substructure at the LHC can conclusively push beyond the 100 GeV threshold set by previous limits, leaving no gaps. In fact, it should be possible to achieve discovery-level sensitivity at 100 GeV using 2012 data, demonstrating the LHC's far superior production cross sections and luminosity compared to earlier experiments. Exclusion-level sensitivity extending up to 300 GeV or higher seems achievable. These results can be obtained using a number of promising data-driven search techniques, and should be realistic if at least one of these techniques exhibits managable systematic errors. Returning to the usual MFV-inspired assumption that most decays contain b-quarks, the option of an extremely powerful b-tagged search opens up. This may provide roughly 10σ discovery sensitivity at 100 GeV, and exclusions extending to nearly 400 GeV, without accounting for possible re-optimization of the analysis cuts or further gains from applying a double-b-tag. The upcoming Run II of the LHC should push exclusion to masses about twice as high, even while further tightening sensitivity at lower masses.
These findings have yet again illustrated the utility of viewing complicated hadronic activity through the lens of jet substructure. This basic change in philosophy from canonical jet-based reconstruction (almost) frees us of the notion of a "minimum distance" between reconstructable hard partons, namely the jet radius, which was in fact a major botteleneck in the traditional 4-jet versions of this search. Clearly with a more flexible viewpoint, the extremely high energies available at the LHC can be made to work for us, not against us, when attempting to search for O(100 GeV) or even O(10's of GeV) objects decaying into jets. Indeed, we have found here that we are in a good position to test whether supersymmetry has exploited the limitations of our conventional analyses, and has been hiding in plain sight this whole time.
of events without light partons in the final state, and can be made to yield a good fit to the matrix element predictions by using shower damping (SpaceShower:pTdampMatch=1). (By contrast, we have found that the "wimpy shower," as well as the virtuality-ordered and p Tordered showers of PYTHIA6, all tend to underestimate the amount of production radiation.) Using self-contained Pythia8 simulations with shower damping, and running through our complete analysis chain described below, we obtain very good agreement with the more complicated matched simulations. Running matching within Pythia8 is also possible, but we have not explored this option for our signal generation.
Our nominal QCD continuum simulations are based on MadGraph5 matched up to four partons within Pythia8, using the more theoretically-rigorous CKKW-L prescription [56]. 13 For our merging measure, we use the Durham-k T distance used internally by MadGraph5, 14 and exploit the program's ability to produce multi-parton simulations with kinematic cuts defined in that space with threshold xqcut. 15 The measure only applies to partons that can realistically be viewed as merging within a QCD diagram according to quark flavor. For example, two antiquarks would never be compared, nor would a u-quark and a c-antiquark. We also forbid mergings with the beam if they would violate flavor, though this occurrence is very rare. Our merging scale is set to 50 GeV for 8 TeV simulations, and 100 GeV for 14 TeV simulations. (For the former, the m avg range over which we run our stop searches is highly dominated by the hard 4-parton events.) For comparison, we have also generated a set of unmatched QCD simulations wholly within Pythia8, based on showered 2 → 2 production with the default (wimpy) shower. These display a qualitatively similar m avg spectrum to the matched sample, though with a less steep falloff, and O(1) higher rate at m avg > ∼ 200 GeV. These differences are much less pronounced when comparing samples processed through ordinary jet reconstruction.
For our other backgrounds (tt, W +jets, t/s-channel single-top, tW , diboson), we mainly relied on Pythia8, though our tt simulation starts as a 2 → 6 decay chain processes in MadGraph5 in order to capture spin correlation effects. We have also damped the shower for tt, which is treated as a power shower by default. Similar to stop pair production, comparison to the kinematics of undecayed matched samples shows good agreement.
All of our simulations are leading-order, and can be normalized to higher-precision calculations using K-factors. For our matched stop samples, we find that a flat K-factor of 1.5 corrects us to the NLO+NLL predictions [76] for almost any mass. We have also verified, by comparing with p T (t) spectra predicted by Prospino 1.0, that the matched spectrum is in excellent agreement with the NLO spectrum. Therefore we do not anticipate a substan-13 v1.7.6 of Pythia8 is affected by a bug that causes crashes for this matching. We thank Stefan Prestel for providing us with a private fix. 14 Beam-distances are computed as k 2 T,iB = p 2 T,i + m 2 i . Inter-parton distances are computed as k 2 T,ij = max(m 2 i , m 2 j ) + 2 min(p 2 T,i , p 2 T,j ) × (cosh ∆η ij − cos ∆φ ij ). 15 To avoid double-counting of α s reweightings, we have commented-out the relevant code in MadGraph5.
tially different K-factor for boosted stops. For tt, we use a K-factor of 1.8. For all of our other simulations, including matched QCD, we coarsely assume K-factors of 1.5. Our own comparisons of matched and unmatched QCD simulations to the data of [9] suggest that this choice of K-factor may be conservatively large. Downstream of the hadron-level simulation, we apply a simple detector model in the form of a 0.1 × 0.1 calorimeter grid in η-φ space. We form our final H T trigger by first clustering the cells into R = 0.5 anti-k T jets with FastJet3 [54], and demanding that the sum of p T s of jets above 50 GeV exceeds 900 GeV (1600 GeV for 14 TeV simulations). At this stage we neglect both jet energy measurement fluctuations and pileup effects. The former should be only a few percent for such a large sum-over-energies, and the latter would be systematically accounted for in a realistic jet measurement.
While we do not concern ourselves with the effects of pileup on ordinary jet energy measurements, we do wish to make sure that our subsequent jet substructure methods are not adversely affected. Therefore, on top of each hard event surviving the H T trigger, we superimpose an average of 20 min-bias events from Pythia8, including both charged and neutral activity. We then actively remove this pileup using a slightly modified form of trimming [58]. We cluster all calorimeter cells in the event into R = 0.2 anti-k T jets, and discard the contents of these jets if their total p T falls below an absolute cutoff of 5 GeV. Unlike canonical trimming, which works jet-by-jet and operates with a relative p T measure, we have targeted this method to subtract the contaminating energy density of pileup, at least in regions that do not overlap with hard activity. We have found that this procedure successfully preserves our 100 GeV stop lineshape, which is otherwise shifted by 5 GeV and broadened due to the pileup. The stop reconstruction rate is unchanged. The final impact on the background is generally modest, though there is an O(1) reduction of the high-m avg tail, e.g. in the vicinity of 300 GeV.
Following this pre-trimming stage, we apply our main clustering into R = 1.5 C/A fatjets. Only the leading two fat-jets in p T are considered, and these must both fall in the region |η| < 2.5. The fat-jets are declustered into subjets as described in Section II. To obtain more realistic stop mass peaks, we smear the energies of the subjets as 16 To cover scenarios where stops have large branching fractions to b-quarks, we also run a b-tagged analysis. We assume that subjets can be tagged similar to ordinary jets, and that the tag/mistag efficiencies are fairly flat in our analysis range. (Our subjet p T 's typically vary between 100 GeV and 400 GeV. See, e.g., [77].) To perform flavor tagging, we keep track of bottom-hadrons and prompt charm-hadrons from the event record, and match them to the closest subjet within ∆R < 0.2. Each subjet's "true" flavor is then determined by the heaviest associated hadron. We apply flat b-tagging efficiencies of 60%, 10%, and 2% for bottom-flavored, charm-flavored, and unflavored subjets, respectively.

Appendix B: Supplementary Results
This appendix contains three supplementary sets of results: the ∆R distributions of subjets for signal events, a comparison of our nominal R = 1.5 jet radius to R = 0.8, and comparisons with the more standard BDRS declustering procedure. Fig. 7 shows the ∆R distributions of subjets within stop-jets, for events passing our complete set of analysis cuts. This plot makes it clear that for mt = 100 GeV, a large fraction of stop decays would comfortably sit inside of a normal-sized LHC jet of R = 0.4 or R = 0.5. It is also notable that, even though we choose a much larger fat-jet radius, very few stop decays are reconstructed with unphysically-large ∆R. In other words, our substructure procedures and analysis cuts adaptively find the "correct" ∆R scale for the signal. For larger stop masses, the separation becomes large enough that an ordinary jet radius could resolve the decays. But in our treatment this regime is continuously connected to the scenarios with ∆R < 0.4, with no artificial threshold. Finally, we can see that with our absolute and relative energy cuts, mt = 300 GeV is about the largest mass that displays complete containment within R = 1.5 fat-jets. Still, a large fraction of mt = 400 GeV decays remain contained, a signal which is important for the b-tagged version of the analysis.
In Fig. 8, we compare the QCD continuum's m avg spectrum with our nominal R = 1.5 to an identical analysis with R = 0.8. It can be seen that, in the vicinity of m avg = 100 GeV, the background increases both in absolute rate and in steepness. Essentially, the entire spectrum has been "squashed" by a factor of 2, since the overall mass scale is set by R × H T . Performing the same analysis with the mt = 100 GeV signal, the lineshape is practically unaltered, but the overall acceptance increases by 30%. This is because, with a narrower fat-jet, there are fewer cases where the declustering picks up a spurious ISR jet. Still, the gain in S/ √ B is marginal, and comes at a cost of slightly reduced S/B in addition to a more difficult background shape. Higher stop masses display significantly reduced efficiencies due to incomplete containment.
Our nominal declustering procedure judges splittings based on the p T 's of subjet candidates relative to the original fat-jet, and on their individual m/p T ratios. This procedure is a direct descendant of the BDRS procedure of [41], which uses a somewhat different set of declustering criteria, and also applies an additional filtering step by reclustering the subjet constituents. Within the kinematic regime of our analysis, the declustering stage of BDRS acts almost identically to our procedure without the m/p T requirement. 17 With 17 The BDRS mass-drop criterion is mostly redundant here, and the declustering is driven mainly by the filtering, the two subjets are further refined into three, using the C/A algorithm with momentum-asymmetry criterion. See [78] for a detailed related study. R = min[∆R(subjets)/2, 0.3]. We reform these back into two subjets, by clustering together the two that are closest in ∆R. This allows us to apply our final cut on the ratio of subjet p T 's, which assumes 2-body substructure. We also consider a form of BDRS without the filtering step. Figs 9 and 10 show a comparison of our nominal procedure with both filtered and unfiltered BDRS. (This comparison is made without pileup or trimming.) It is clear that the 100 GeV stop signal is fairly insensitive to the detailed procedure, but that the 300 GeV stop and QCD spectra can be highly reshaped, and that the overall rates in the vicinity of m avg = 300 GeV are increased by O(1). The filtered QCD spectrum is also flatter in the region between 100 GeV and 200 GeV, which could help a signal bump stand out more clearly there. These last two points can be considered advantages. However, filtering also accidentally introduces new mass scales from its minimum reclustering radius of 0.3 for well-separated subjets, resulting in a bimodal QCD distribution with a local minimum near 50 GeV and a broad local maximum near 175 GeV. And while the 300 GeV stop signal is enhanced, it sits on top of a background that is starting to sharply change shape. Therefore, these advantages must be treated with caution. Removing filtering, we still see some enhancement of the 300 GeV signal and its background, but the background shape becomes less biased. The differences between unfiltered BDRS and our nominal declustering are dominated by the introduction of high-m/p T subjets for the former, and these are more likely to be contaminated by additional radiation. We have also found that the detailed shape of the spectrum at high mass develops more sensitivity to changing analysis cuts to define control regions. However, these differences relative our nominal procedure might be reduced with the use of trimming.