Magnifying the ATLAS Stealth Stop Splinter: Impact of Spin Correlations and Finite Widths

In this paper, we recast a"stealth stop"search in the notoriously difficult region of the stop-neutralino Simplified Model parameter space for which $m(\tilde{t}) - m(\tilde{\chi}) \simeq m_t$. The properties of the final state are nearly identical for tops and stops, while the rate for stop pair production is $\mathcal{O}(10\%)$ of that for $t\bar{t}$. Stop searches away from this stealth region have left behind a"splinter"of open parameter space when $m(\tilde{t}) \simeq m_t$. Removing this splinter requires surgical precision: the ATLAS constraint on stop pair production reinterpreted here treats the signal as a contaminant to the measurement of the top pair production cross section using data from $\sqrt{s} = 7 \text{ TeV}$ and $8 \text{ TeV}$ in a correlated way to control for some systematic errors. ATLAS fixed $m(\tilde{t}) \simeq m_t$ and $m(\tilde{\chi})= 1 \text{ GeV}$, implying that a careful recasting of these results into the full $m(\tilde{t}) - m(\tilde{\chi})$ plane is warranted. We find that the parameter space with $m(\tilde{\chi})\lesssim 55 \text{ GeV}$ is excluded for $m(\tilde{t}) \simeq m_t$ --- although this search does cover new parameter space, it is unable to fully pull the splinter. Along the way, we review a variety of interesting physical issues in detail: (i) when the two-body width is a good approximation; (ii) what the impact on the total rate from taking the narrow width is a good approximation; (iii) how the production rate is affected when the wrong widths are used; (iv) what role the spin correlations play in the limits. In addition, we provide a guide to using MadGraph for implementing the full production including finite width and spin correlation effects, and we survey a variety of pitfalls one might encounter.


Introduction
One of the main drivers of the Large Hadron Collider (LHC) program is to understand the origins of the electroweak scale. In the Standard Model, the Higgs boson mass parameter m 2 H is quadratically divergent, implying that it is not calculable within the Standard Model effective theory since it is sensitive to physics at arbitrarily high scales, e.g., the Planck scale where gravity becomes strong. Fine-tuning away these quadratic contributions leads to the so-called hierarchy or naturalness problem. One compelling way to ameliorate this problem is to assume that nature is Supersymmetric, thereby yielding a calculable Higgs boson mass parameter. Then softly breaking this symmetry maintains the calculability of the theory, at the expense of introducing physical quadratic corrections to m 2 H . Furthermore, every particle now gets a partner (whose spin differs by a half integer), implying a rich phenomenological program at the LHC to hunt down these new states. See [1] for the classic review of Supersymmetry and [2] for some perspective on the status of the naturalness problem.
While the masses of the superpartners are free parameters, the couplings are fixed by the symmetries of the theory. This implies that, given limits on a mass spectrum, one can compute the fine-tuning required to reproduce the measured electroweak scale. Unsurprisingly, the two superpartners with the largest couplings play the most important role: the superpartners of the top -the stop -and the gluonthe gluino. The mass of the partner of the Higgs is also relevant, but the limits on its mass are sufficiently weak due to a low production cross section such that it is not as relevant as the stop and gluino. Taking the connection between the spectrum and the fine-tuning problem seriously motivates that the stop, gluino, and Higgsinos are the lightest new physics states that would appear at the LHC [3,4]. This spectrum has become even more compelling in the light of the comprehensive set of constraints on new physics published by both ATLAS and CMS, as emphasized in [5][6][7]. This has additionally motivated many studies recasting ATLAS and CMS results to comprehensively understand the connection between the allowed parameter space for Higgsinos, stops, and gluinos, and the required tuning, e.g. [8][9][10].
The focus of the paper here will be on the stop-neutralino (t 1 -χ 0 1 ) Simplified Model, with a focus on the status of the parameter space currently constrained by ATLAS. In particular, we want to explore what is still allowed in the very difficult "stealth stop" region where m(t 1 )−m(χ 0 1 ) m t , where m t is the top quark mass. The compressed kinematics make this particularly difficult. Although the final state is ttχ 0 1χ 0 1 , the neutralinos are nearly at rest so that the resulting missing energy is very small and thus the observable signature is a pair of top quarks. Furthermore, the production cross section for stop pair production is O(10%) of the top pair production cross section. Given these challenges, ATLAS has done a remarkable job constraining the stealth parameter space, and at this point only a "splinter" remains [11] as shown in the gray region of Fig. 1.
In more detail, ATLAS constraints on stealth stop production are derived through  The gray region summarizes the ATLAS constraints, not including the tt spin correlation constraint, or the σ tt contamination constraint (which is the subject of this paper) [11]. several independent analysis approaches: searches for direct stop pair production where m(t 1 ) m t ,t 1 → t +χ 0 1 and the top decay is on-shell [12,13]; searches for direct stop pair production where m(t 1 ) − m(χ 0 1 ) m t and the stops each decay viã t 1 → W bχ 0 1 [14]; analysis of tt spin correlations through an angular analysis of the lepton decay products [15]; and a simultaneous measurement of the tt production cross section at √ s = 7 TeV and 8 TeV [16]. The latter two analyses constrain stop pair production in the stealth region through detailed measurements of tt processes. CMS has performed similar searches, e.g. [17][18][19][20][21][22]. There have additionally been some relevant phenomenology studies providing complementary ideas for targeting stealth stops [23][24][25][26][27][28][29][30].
Our focus here will be to recast the limit derived from σ tt [16] into the full m(χ 0 1 )m(t 1 ) plane -ATLAS only provides results for fixed m(χ 0 1 ) = 1 GeV. A careful theoretical treatment is needed in order to properly simulate events in the region where m(t 1 ) − m(χ 0 1 ) m t . In particular, a naive factorizing of the production and decay for the stop production misses important physical effects from the finite size of the widths involved, and from spin correlations. Our goal here is to review these issues in detail, with an emphasis on their impact for recasting the results of [16]. We also provide a guide for the reader for simulating events near threshold using modern Monte Carlo event generation tools MadGraph and Pythia. These lessons are more generic than the application presented here, and our hope is that the reader will find them useful when simulating events in a region where mass thresholds are relevant.
This paper is organized as follows. A review of the theoretical subtleties that arise when simulating stealth stop production/decay (or more generally for generating events near a mass threshold) through the splinter is presented in Sec. 2; additional technical details can be found in App. A and some common pitfalls one might encounter are discussed in App. B. Then we provide a review of the ATLAS top cross section measurement, and stealth stop limit in Sec. 3. Emphasis is given to the CL s limit setting procedure utilized here, and a validation of our framework is presented. Finally, Sec. 4 provides the results of our recast and shows how much of the splinter can be excluded. We further emphasize the role of the systematic error on the top cross section by showing how the limits can be strengthened if this is reduced. A brief closing discussion is given in Sec. 5.

Calculating Stealth Stop Pair Production
A self-consistent framework for generating events is required to properly reinterpret the experimental constraints in the stealth stop region where m(t 1 ) − m(χ 0 1 ) m t . We address four effects that could impact the recasted limit: carefully computing the decay width in the two-to-three-body transition region, the error introduced by the narrow width approximation (NWA), the importance of including width in the Breit-Wigner propagator, and how spin correlations can change the efficiency to pass the signal region cuts. These are each detailed in Secs. 2.1 to 2.4.
Three approaches are compared. All three use MadGraph5_aMC@NLO 2.6.1 [31] for the hard matrix element calculation and convolution with the NNPDF3.0 leading order parton distribution functions [32], and Pythia 8.2 [33] is always used for the parton shower and hadronization. The stop pair production cross section is calculated using NLLFast [34][35][36]; we additionally correct for the numerical impact of the NWA (see Sec. 2.2) when using the NLLFast cross section to compute our final constraints. Detector effects are approximated using Delphes 3.4.1 [37]. We use the default delphes_card_ATLAS.tcl card, modified so that the jet radius is 0.4, in accordance with ATLAS. Jets are clustered with the anti-k t algorithm [38,39] within the FastJet 3.2.1 [40] framework. More information regarding the event generation is given in App. A.
The differences come into what approximation is assumed when decaying the stops. The relevant Feynman diagrams for stop production 1 and decay are given 1 Here we will follow the ATLAS conventions and express stop pair production as p p →t 1t1 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 Y W a e n w q u w L I + j M r K 6 t r 6 R n W z t r W 9 s 7 t X 3 z + 4 U 0 k m K X N o I h I 5 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 4 7 1 5 N L 2 v q V g N x 9 / X Q Y 7 T o C 5 m I c = " > A A A B 6 X i c d V B N S 8 N A E N 3 U r 1 q / q h 6 9 L B a h e g i b U r W 9 F b 1 4 r G h s o Q 1 l s 9 2 0 S z e b s L s R S u h P 8 O J B x a v / y J v / x k 0 b Q U U f D D z e m 2 F m n h 9 z p j R C H 1 Z h a X l l d a 2 4 X t r Y 3 N r e K e / u 3 a k o k Y S 6 J O K R 7 P p Y U c 4 E d T X T n H Z j S X H o c 9 r x J 5 e Z 3 7 m n U r F I 3 O p p T L 0 Q j w Q L G M H a S D f V k + N B u Y L s Z u O 0 j h B E N p o j I 7 W z Z q 0 O n V y p g B z t Q f m 9 P 4 x I E l K h C c d K 9 R w U a y / F U j P C 6 a z U T x S N M Z n g E e 0 Z K n B I l Z

z Q K c A h H U A E P L q A G t 1 A H H x j 0 4 R l e 4 c 2 R z o v z 7 n z M W p e c f O Y A / s D 5 / A G j a o z n < / l a t e x i t > (⇤)
< l a t e x i t s h a 1 _ b a s e 6 4 = " C P 1 Q p j g W P A U U A 9 6 i F e T e y S 6 x g 7 Q = " < l a t e x i t s h a 1 _ b a s e 6 4 = " C P 1 Q p j g W P A U U A 9 6 i F e T e y S 6 x g 7 Q = "

z Q K c A h H U A E P L q A G t 1 A H H x j 0 4 R l e 4 c 2 R z o v z 7 n z M W p e c f O Y A / s D 5 / A G j a o z n < / l a t e x i t >
< l a t e x i t s h a 1 _ b a s e 6 4 = " C P 1 Q p j g W P A U U A 9 6 i F e T e y S 6 x g 7 Q = "

z Q K c A h H U A E P L q A G t 1 A H H x j 0 4 R l e 4 c 2 R z o v z 7 n z M W p e c f O Y A / s D 5 / A G j a o z n < / l a t e x i t >
< l a t e x i t s h a 1 _ b a s e 6 4 = " C P 1 Q p j g W P A U U A 9 6 i F e T e y S 6 x g 7 Q = "

z Q K c A h H U A E P L q A G t 1 A H H x j 0 4 R l e 4 c 2 R z o v z 7 n z M W p e c f O Y A / s D 5 / A G j a o z n < / l a t e x i t >
< l a t e x i t s h a 1 _ b a s e 6 4 = " C P 1 Q p j g W P A U U A 9 6 i F e T e y S 6 x g 7 Q = " < l a t e x i t s h a 1 _ b a s e 6 4 = " C P 1 Q p j g W P A U U A 9 6 i F e T e y S 6 x g 7 Q = "    Figure 2: The left diagram illustrates the MG Full scheme and the right diagram illustrates the P NWA scheme for the parameter space where m(t 1 ) − m(χ 0 1 ) < m t (see Table 1 for an explanation of the naming conventions). The green circles represent the full tree-level stop pair production matrix element. The gray circles represent decays that do not include any matrix element information, i.e., the particles are decayed using phase space alone. The superscripts ( * ) denote particles that can go off shell. This figure was adapted from the diagrams in [14].
in Fig. 2. In both the left and right diagrams, the circle connecting the proton lines represents all of the possible diagrams that contribute to stop pair production, assuming the only particles beyond the Standard Model are the lightest stop and lightest neutralino. In the left diagram, the stops decay to a neutralino and a top quark, which can be off shell. The top quark then decays to a bottom quark and a W boson, which further decays leptonically (so as to populate the signal region). This is the diagram for the most complete calculation (which we use for computing the recasted limits below), implemented using MadGraph to compute the full parton-level final state (before showering), and is referred to as MG Full . It can also be interpreted as the calculation performed by MadSpin [41] (the approach taken by ATLAS), which computes spin correlations taking all intermediate particles exactly on-shell; we refer to this approach as MS. The final approach (P NWA ) is to allow Pythia to decay the stops, which is equivalent to assuming a flat matrix element i.e. neglecting any spin effects in the final state. This is illustrated by the diagram on the right, where the stop squark goes through a three-body decay directly. The naming conventions for these three approaches are summarized in Table 1.
There are many subtleties that must be accounted for to correctly evaluate production for the parameter space near the three-to two-body threshold. Three-body effects on the width can be sizable, even if the two-body final state is open. The narrow width approximation yields percent level errors with regard to the full production cross sections and/or width calculations. If one uses the wrong width, this can have a non-trivial impact on the total rate. Finally, spin correlations affect the distributions of the final state particles, and are particularly important when instead of p p →t 1t1 * . the three-body decay dominates. The next four subsections explain these points in detail.

Transition from Three-body to Two-body Decays
If the channel is open, two-body decays nearly always dominate three-body decays due to phase space suppression. However, this transition must be accounted for continuously, which requires that care must be taken near threshold. GeV and m(χ 0 1 ) = 1 GeV. The solid blue and dashed green lines are computed with MadGraph, by forcing a two-and three-body decay respectively. Clearly, the twobody approach does not capture near threshold effects (which can be important with multiple widths of the top as denoted by the thin vertical lines), due to the phase space factor implicit in the two-body assumption. For reference, we also show a comparison with SDecay in the red dot-dashed line, since it includes loop corrections from QCD. However, clearly SDecay is not reliable near threshold. In Fig. 3, the stop decay width is plotted as a function of the stop mass. The blue and dashed, green lines are computed with MadGraph, and are the two-or three-body decays, respectively. Well above threshold, the two-and three-body decays yield the same width, which implies that the three-body process is dominated by the onshell (as opposed to an off-shell) top quark. However, as the mass of the stop is lowered towards the threshold, the two-body processes start to diverge as there is less available phase space for the two-body decay. The darkest vertical line denotes the threshold, while the lighter ones mark one, two, and three top widths above threshold. We see that the stop mass needs to be many top widths above threshold for the off-shell top to not contribute significantly.
While Fig. 3 shows that off-shell effects are important even above threshold, not all tools will carefully take this into account since their purpose is usually to cover the parameter space where there are no accidental degeneracies. For instance, the red line displays the width of the stops as calculated by SDecay [42]. As soon as the stop is above threshold, the program only computes the two body diagram, leading to a discontinuity in the stop width. Similarly, using the default compute_widths command during event generation in MadGraph also only computes two-body modes anywhere above threshold. The offset between the SDecay and MadGraph results is due to the inclusion of QCD corrections to the width by SDecay. The impact of incorrectly computing the stop width is explored further in Sec. 2.3 and App. B.
The essential physics is straightforward to understand. An unstable particle with mass M and width Γ contributes to a process P via a propagator: where q is the momentum flowing through the propagator of interest, and M(q 2 ) represents the rest of the (integrated and spin averaged) matrix element squared for P. Note that P captures both decay and production. When the width is small compared to the mass of the particle, the propagator is highly peaked near q 2 ∼ M 2 . The rest of the matrix element can then be evaluated at q 2 = M 2 using the integrated propagator: Thus in the NWA, P is given by: As a toy example with which we can investigate the quality of the NWA, we compute scalar production through an s-channel propagator convolved with the gluon pdf. 2 The production is given by the full propagator: where A is a constant that includes the couplings and phase space factors, Y is the rapidity of the final state particles, and √ s is the center of mass energy of the protons. For comparison, the process in the NWA is: Eqs. (2.4) and (2.5) are evaluated using the MSTW 2008 NNLO parton distributions [48][49][50] at a center-of-mass energy of √ s = 8 TeV for a variety of masses and widths.
As the gluon distribution grows at low momentum fractions, we cannot perform the momentum integral in Eq. (2.4) over the full phase space. To mitigate this issue, we fix the limits of integration in terms of a given number of widths around the resonance peak. Figure 4 shows the ratio of the NWA of Eq. (2.5) to the full calculation of Eq. (2.4) as a function of the width-to-mass ratio Γ/M . The integration range is taken to be 5 or 15 widths in the left and right panels, respectively, and the different colors show different choices of resonant mass. Due to the choice for fixing the limits of integration in terms of a number of widths, there is an intrinsic limit to how large the width-to-mass ratio can be; if the width is too large, the integral range implies imaginary momenta. When this occurs (as marked by the vertical dashed lines in the figure), a constant starting lower bound for the integration is used.
For the simple toy process modeled here, at small widths the NWA overestimates the rate by 6% or 2% for integration windows of 5 or 15 widths, respectively. The NWA is no longer effective when the width is a few percent of the mass, and larger masses imply a faster breakdown. This seemingly small effect can be important as the NWA is commonly used, and this discrepancy factor naively contributes for each NWA assumed in a diagram. For instance, producing the particles on-shell with subsequent decays left to other programs implicitly assumes the NWA as it does not integrate over the intermediate on-shell propagator. As shown in the appendix, this yields a 6% effect on the total cross section for stop pair production. The NWA is also usually assumed when there are particle decay chains, such as when computing the width of the top quark -integrating over the intermediate W from a top decay leads to a few percent effect.  The left and right panels show the result when the integration limits are 5 or 15 widths around the resonance peak, respectively. The vertical dashed lines show the value of the width-to-mass ratio which would lead to unphysical imaginary momenta in the integration; values to the right use a constant lower bound for the integration instead. We see that the NWA is no longer a good approximation for Γ/M 10 −2 , and that this breakdown happens earlier as the mass is increased.

Inconsistent widths
As shown above, common practice approaches to width calculations near threshold can yield errors near a factor of 10. In this subsection, we examine the numerical impact on the production rate if one uses the wrong width. We will specialize to the stop pair production and decay process of interest: where we do not put any kinematic requirements on the internal top or W .
By inspection of Fig. 2, two stop propagators must be considered. Since these widths are very small, we can reliably estimate how an error in the stop width propagates to the full process using the NWA in Eq. (2.3). If the stop width is computed inconsistently by an amount ∆Γ, then including the effect of both stop propagators implies that the cross section prediction will be wrong by approximately This effect is illustrated numerically in Fig. 5  We see that in order to accurately estimate the cross section and efficiencies using the full matrix element, it is important to consistently calculate the widths. Due to the presence of two stop propagators in the production diagrams, inconsistencies in the widths lead to changes in the rates by ∼ 2 ∆Γ/Γ.

Spin correlations
The last effect explored here, which is especially important for light stops near or below the top threshold, is the spin correlations of the decay products. When events are simulated in the NWA, i.e., as is usually done when decaying particles in Pythia, production is factorized from decay. This implies that the angular distribution of the decay products do not include the correlations that result from the spin of the parent. While it is possible to include matrix elements for decays in Pythia 8, the three-body decay of the stop is not currently implemented. When Pythia decays particles without a matrix element implementation, it assumes a constant matrix element and determines the kinematics strictly from phase space.  Table 1 above for conventions). The data sets are normalized and binned in the plane given by the p T of the second hardest particle versus the p T of the hardest particle for bottom quarks b, charged leptons , the neutralinoχ 0 1 , and the neutrino ν. The color scales show the result of subtracting the P NWA bins from the MG Full bins. P NWA events do not include spin correlations which results in the b quarks (upper left) being softer than they are when the full matrix element is used. This is compensated by harder neutralinos (lower left). There is some effect on the lepton and neutrino distributions, coming from the subsequent decay of the W , but the result is more subtle.
As an example for the importance of including spin correlations, Fig. 6 examines the distributions of the final state particles for pair production of 150 GeV stops which decay to a 1 GeV neutralino, b quark, lepton, and neutrino. We then compare the entire process, including decays, generated in MadGraph (MG Full ) with a calculation where the stops are decayed by Pythia (P NWA ). The events are then binned according to the harder/softer truth-level b quark, lepton, neutralino, and neutrino, normalized such that the sum of the bin content multiplied by the bin area is unity. The figure shows the result of taking the content of the MG Full bins and subtracting the content of the P NWA bins. Bins with more MG Full content are shaded blue, while the red bins contain more P NWA events.
Physically, the results of Fig. 6 are due to the fact that the the outgoing particles for the P NWA events, are distributed according to phase space. This neglects the fact that the b and theχ 0 1 are connected to the same scalar vertex and should therefore be correlated. This results in b quarks from the P NWA events being softer than from the MG Full events. In contrast, the neutralinos are softer when the full matrix element is considered. Now that we have reviewed a variety of physical effects that can impact the production calculation, we will move on to discuss how the ATLAS top cross section measurement is performed and how the limit on stealth stop production is derived.

Recast of the ATLAS Measurement and Search
This section begins with a review of the ATLAS measurement [16] of the tt cross section that is used to probe the stealth stop region. The non-trivial statistical procedure is explained, and a validation of our simulation framework is provided. For details of the simulation framework, see Sec. 2 above.
The ATLAS constraint on stealth stops [16] is derived by placing a limit on the allowed contamination when measuring the tt cross section σ tt . The approach is to extract σ tt using both √ s = 7 TeV and 8 TeV data in the e ± µ ∓ channel. An additional selection requiring either exactly one or two b-tagged jets is applied. A simultaneous fit (which additionally allows the top mass to vary) utilizing all four bins is then performed to extract both the cross section and the b-tagging efficiency. The benefit of using a simultaneous fit is to minimize the effect of systematic uncertainties common to the cross-section measurements at the two masses. The result is σ tt (8 TeV)/σ tt (7 TeV) = 1.326 ± 0.057, which is consistent with the Standard Model prediction [16,[51][52][53][54][55][56].
If one is willing to fix the top quark mass using independent observables, the measurements performed in [16] can be interpreted as limiting the number of events coming from the pair production of stop squarks in the stealth stop region; see [25,26] for related pheno studies. ATLAS uses this method to exclude a right-handed stop with masses between the top quark mass and 177 GeV, but for fixed neutralino mass of 1 GeV. We reinterpret their result to constrain regions of the m(χ 0 1 ) − m(t 1 ) plane, including stop masses below the top quark mass. 3 We will present a statistical approach to handle the suite of correlated and un-correlated systematic errors, yielding a 95% CL exclusion region using the CL s method. As emphasized in the previous section, we also carefully include off-shell effects and spin correlations.

Summary of the ATLAS tt Cross Section Measurement
This search relies on electrons, muons, and b-tagged jets with the following kinematic requirements. An electron candidate must have a transverse momentum of p T (e ± ) > 25 GeV and pseudorapidity |η(e ± )| < 2.47. In addition, electrons near the transition region between the barrel and endcap (1.37 < |η(e ± )| < 1.52) are removed. Muons are required to satisfy p T (µ ± ) > 25 GeV and |η(µ ± )| < 2.5. Jet candidates must have p T (j) > 25 GeV and |η(j)| < 2.5. We utilize the b-tagging efficiencies provided by Delphes for our stop signals; there is a different treatment of b-tagging for the tt prediction, as detailed in what follows. The events are triggered on either a single electron or muon, and the efficiencies of each of these triggers plateaus by a transverse momenta of 25 GeV.
Each event considered in the measurement must have exactly one electron and one muon of opposite sign. Events are further classified as having exactly one or two b-jets (with no requirement on the number of non-b-tagged jets) resulting in N b and N bb number of events respectively. The expected number of events in these channels is expressed as where L is the integrated luminosity, eµ is the efficiency for a tt event to pass the e ± µ ∓ preselection (including geometry and detector effects); ATLAS provides eµ = 7.7 × 10 −3 . Additionally, b is the efficiency to tag a b-jet, specifically in events coming from top quark decays. Then naively the probability of tagging the b-quarks coming from each of the two tops in an event would be 2 b . However, kinematics and detector responses can lead to correlations for the two b-jets. This is taken into account by defining C b ≡ bb / 2 b , where bb is the probability of tagging both b-quarks in a tt event. Finally, there are additional Standard Model processes that contribute to both N b and N bb ; these are denoted by N bkg b,bb . This is dominated by single top production, which mainly contributes to N bkg b . The expressions for N b,bb in Eq. (3.1) can be solved for the tt cross section and the b-tagging efficiency, To extract σ tt and b , the values of C b and eµ are estimated from simulations and given in Table 2. The cross section measurement is affected to a greater extent on the uncertainties from N bkg b than from N bkg bb . Conversely, the uncertainty in N bkg bb has a bigger impact on b . The measurement yields σ tt = 182.9 ± 7.1 pb at √ s = 7 TeV and σ tt = 242.4 ± 10.3 pb at √ s = 8 TeV [16]. The result of the two measured cross sections can then be used to infer the mass of the top quark; ATLAS gives the best fit as m t = 172.9 +2.5 −2.6 GeV [16].

Validation of Constraint on Stealth Stop Production
The tt cross-section measurement is sensitive to new physics that decays to top quarks. As opposed to the cross section measurement, where m t is left as a free parameter, taking m t from other measurements allows the use of the predicted value of σ tt to constrain the possibility of contamination from new physics. The ATLAS limits for stop pair production are set by simultaneously fitting the 7 and 8 TeV datasets and using profile likelihood ratios in the asymptotic limit [58]. Correlated uncertainties are accounted for through the use of nuisance parameters. For instance, the dominant top quark pair production cross section uncertainties have a common origin, e.g. the top mass uncertainty, and are thus treated as fully correlated. We take a simplified approach by modeling these uncertainties as Gaussians with widths of the averages of the upper and lower uncertainties for each distribution. To enforce that the 7 and 8 TeV predictions shift together, we introduce the nuisance parameter δ tt as: Given the cross section, the number of expected events in a given bin coming from top quark pair production can be computed using the first terms in Eq. (3.1). The additional nuisance parameters are given in Table 2.
In a similar manner, we treat the non-tt Standard Model contributions to the one and two b-tagged bins as correlated between the two energies. ATLAS provides the number of expected events and size of the uncertainties in each of the bins. We translate this to cross sections, such that the uncertainty in the number of events is the result of the cross section and luminosity uncertainties added in quadrature. The resulting cross sections are then given by where δ N b and δ N bb are the nuisance parameters.   [16]. Each parameter is chosen from a Gaussian centered at the mean with a width given by the uncertainty. The luminosity L has units of fb −1 and all other parameters are dimensionless.
Finally, we can compute the contribution from stop pair production to the signal regions. The central value of the stop production cross section is taken from NLLFast [34]. The efficiencies of the selection criteria described in the previous section are estimated with the simulated events and are denoted by ξ b or ξ bb for the one and two b-tagged bins, respectively. The number of expected events from stop pair production in each of the four bins can then be given by Nt 1t1 (7,8) (1,2) = L (7,8) where the superscripts denote the center of mass energy in TeV and the subscripts denote the number of b-tags. The 7 and 8 TeV predictions are correlated through the cross section and luminosity, but each efficiency is treated as independent. As shown in Table 2, we include a 10% uncertainty on the stop production cross section. We have tested varying this between 5-20% and find minimal changes in the final exclusions.
Our framework uses likelihood ratios L(µ, θ) to set limits on stop production, where µ is the signal strength parameter which is ultimately what gets constrained (µ = 0 is defined as Standard Model and µ = 1 is Standard Model + full stop signal), and θ are the nuisance parameters. The likelihood is defined as the product of the Poisson probabilities of the observed number of events in each of the signal regions N obs i multiplied by the product of the Gaussian probabilities for yielding the specific value of the nuisance parameters, where Nt 1t1 is the predicted number of stop pair production events, N tt is the predicted number of top pair production events, N bkg is the predicted number of non-tt background events, and (mean j , var j ) are given in Table 2. The number of events observed by ATLAS are The likelihood is converted into a test statistic by taking the ratio of the constrained and unconstrained maximum likelihoods. In the unconstrained maximum, both the signal strength and the nuisance parameters are free to vary in order to maximize the likelihood. The values of the signal strength and nuisance parameters at this maximum are denoted byμ andθ, respectively. The constrained maximum likelihood fixes the signal strength to a chosen value, and then the nuisance parameters roam to maximize the likelihood for that value of µ; the values of the nuisance parameters in this case are given byθ. Putting this together allows us to define (3.8) Large value of q represent greater incompatibility with the data. Next, we define a similar quantity under the assumption that only the Standard Model contributes to the data where θ takes the central values for all of the nuisance parameters.
Finally, the CL s method combines q µ and q µ,A . In the asymptotic limit of large number of events, the 95% CL limit is set by solving for the value of µ which yields where Φ is the cumulative distribution of the standard Gaussian with unit width. For more information on this statistical procedure, see App. A of [59].
ATLAS placed limits on the stops assuming a neutralino mass of 1 GeV and a top mass of 172.5 GeV [16,57]. In order to trust our extension of the analysis, we must first ensure that our framework can reproduce the ATLAS results to a good approximation. To this end, Fig. 7 shows the limit obtained by our analysis, i.e., the value of µ which satisfies Eq. (3.10), in comparison with the ATLAS results as a function of the stop mass. The yellow uncertainty band is taken from ATLAS. Values of µ = 1.0 or less imply the stop pair production rate is excluded for that value of the mass. Over most of the region, our analysis follows very closely with the observed limit.
We do note that our results start to deviate from the ATLAS observed limit (staying within the expected uncertainty band) below a stop mass of around 180 GeV. This is not surprising as our statistical methods are different than those used by ATLAS for the signal. ATLAS assumes that the b-tagging efficiency, b , is the same for both tt andt 1t1 events. However, as noted in [57], for low stop masses, the "fitted b is different from what is expected for tt events alone". Trying to keep the tt and thet 1t1 efficiencies correlated in this regime (which is not done in our approach) is likely the reason for the difference. The best fit value of b in our scenario (for tt alone) is shown in Fig. 12.  [16,57] of the signal strength µ vs. the stop mass; our result matches the observed limit over most of the parameter space. As the stop mass drops below 180 GeV, three-body effects and spin correlations start to impact the b-tagging efficiency, and our limit drops to the lower range of the expected uncertainty band as discussed in the text.

Impact of Spin Correlations and Finite Widths
The impact of non-trivial spin correlations are investigated by examining each of the different event generation methods presented in Table 1. The number of expected events coming from stop pair production in each of the four signal regions is shown in Fig. 8 as a function of the stop mass. We sample the stop masses between 150-200 GeV with step sizes of 5 GeV. The blue, green-dashed, and red-dot-dashed lines represent whether the generation is preformed in the MG Full , P NWA , and MS approximations respectively.
The upper panels show the one b-tag regions. In these, the number of expected events grows as the stop mass is decreased, due to the increase in the production cross section. The MG Full and MS predictions, which both include spin correlations, are not particularly sensitive to the position of the threshold as expected from Fig. 3. However, the P NWA line shows a dip at the threshold as it does not account for the three-body to two-body transition correctly. The two b-tag regions are displayed in the lower panel. The effect of the threshold is more dramatic, and the number of expected events is nearly constant below threshold for each method of decay. However, the P NWA approach yields ∼ 20% less events. This is due to the lack of spin correlations for three-body stop decays in Pythia, which results in softer p T on average for the second b-jet.
The derived limits from the three different generation method are shown in Fig. 9. The MS line is slightly lower across the entire region. This is because there are slightly  Table 1.
more events in each panel of Fig. 8 due to the overestimated production cross section coming from the NWA, see Sec. 2.2. Otherwise, since this approach models spin correlations, the results are very similar to the MG Full exclusion.
The P NWA exclusion line is less excluded at large masses, and then flips at the threshold to being more excluded than the approaches that include spin correlations. There are two reasons for this behavior. Even though the cross section is over estimated due to the NWA, the lower right panel of Fig. 10 shows that the softer lepton is more often too soft in comparison to the full calculation. This accounts for there being less events in each signal region above threshold in Fig. 8 pp →t 1t1 → e ± µ ∓ bbχ 1 0χ 1 0 νν MG Full P NWA MS Figure 9: We show the 95% CL limit on the signal strength when the stops are decayed using the three approximations given in Table 1. Beyond the discontinuity around the threshold, there is little difference in the excluded values.
in the weaker exclusion. The strong exclusion for the P NWA events below threshold is unintuitive. As the upper left panel of Fig. 10 shows, below threshold the second b-jet tends to be softer, leading to fewer events, as seen in Fig. 8. It would seem that with fewer events, the exclusion should be weaker than the MG Full line. However, the larger shift in the ratio of the one-and-two b-tagged regions is harder to accommodate in the fit, resulting in the stronger exclusion.
While all of the generation methods result in limits near a stop mass of 180 GeV, Fig. 8 shows that there are thousands of extra events expected from the stops for masses larger than 180 GeV. The uncertainties impacting these regions are explored more in Sec. 4.1. It is surprising that the excluded signal strengths are so similar below threshold, even when the number of events are so different.

Recasted Stealth Stop Results
The results of our recast extends the ATLAS analysis to allow for a variety of neutralino masses. In addition, we interpret the results in terms of both right-and lefthanded stop pair production.
To cover the stealth stop splinter region, events are generate in a grid of the stop and neutralino masses given by: 155,160,165,170,175,180,185,190,195,200 GeV, m(χ 0 1) ∈ 0, 1 × 10 −4 , 1 × 10 −3 , 1, 5, 10,15,20,25,30,35,40,45,50,55,60 GeV . For each parameter point, we generate 250,000 events at each center of mass energy for the full matrix element, and 1,000,000 events when using Pythia for the decays (because we do not specify the decay mode in this case). This is done both for rightand left-handed stops, which decay through t L and t R , respectively. Note that for the left-handed stops, we do not include a light left-handed sbottom in the spectrum.
At each point in parameter space, the value of the signal strength which is excluded at 95% CL s is computed, and a linear interpolation is used to extrapolate between the parameter points. The shaded regions in Fig. 11 show which points are excluded for the full matrix element using the MG Full approximation. The left panel is for stops which decay through a right-handed top, and in the right panel, the stops decay through a left-handed top. The red line shows how the boundary shifts in the NWA without including the spin correlations, the P NWA approximation. These results show that neutralinos in the range 0 GeV < m(χ 0 1 ) 55 GeV for m(t 1 ) 180 GeV are excluded. There are only minor differences in the left-and right-handed stop exclusions when using the full matrix element with spin correlations. The limit extends smoothly between m(χ 0 1 ) = 1 GeV, through the sub-GeV parameter space (which is a region of interest for new ideas in direct detection [60]), to m(χ 0 1 ) = 0 GeV. Spin correlations were not critical for the validation plot Fig. 9 where m(χ 0 1 ) = 1 GeV; the exclusions were very similar whether or not spin correlations were included. However, we see that in the full plane, spin correlations play a role. If they are not included, exclusions are too aggressive at large stop masses, pushing the would-be-excluded region to m(t 1 ) 190 GeV. In contrast, for lighter stops, the bounds are too conservative and do not cover the full area that is excluded when including more physics in the event generation.

Limits with Reduced Systematics
Our recast does not cover the entire stealth stop region; some of the splinter remains embedded. In fact, it turns out that there can be thousands of events coming from stop production contributing to the signal region in the parameter space that we were not able to exclude. In this section, we examine the dominant systematic responsible for causing our limits to saturate and explore how the limits can improve if the uncertainty in the top cross section prediction were reduced. In order to understand its quantitative impact, we have provided Fig. 12 which gives a sense of the relative size of the σ tt uncertainty. Rather than plotting the number of expected events, the We see that in order to compensate for extra events coming from stops (as m(t 1 ) becomes smaller), the tt best fit value can be driven below the error band. In particular, this is true for the best fit for the number of tt events in the gray excluded region. As the mass of the stop increases and the stop production cross section falls, the number of tt events trends toward the central value. The right panel helps expose what is driving these changes by showing the number of standard deviations away from the mean for the best fit values of a subset of the nuisance parameters. In order to accommodate the extra events from stops, the tt cross section is lowered, while the b-tagging efficiencies are simultaneously being driven to larger values. At low stop masses, even the number of events coming from Standard Model backgrounds (which is a small contribution to the total number of events) are pushed outside their uncertainty band.
As shown in Eq. stop events are still not excluded. Measuring the mass of the top quark to better accuracy using unrelated kinematic measurements, as well as improving the pdf and scale uncertainties, could reduce this systematic, thereby yielding a stronger limit. To illustrate the potential impact this could have, the left panel of Fig. 13 shows the value of the signal strength that could be excluded at the 95% CL for a few different mass points in our parameter scan as a function of the improvement factor for the tt cross-section uncertainty. For the massless neutralino parameter points (blue and green), the value of the signal strength which is excluded drops very quickly as the tt cross-section uncertainty is reduced. The green line, with m(t 1 ) = 190 GeV, goes from being allowed to excluded by reducing the uncertainty by less than a factor of 2. However, the exclusion on the signal strength for the m(t 1 ) = 230 GeV parameter point never gets below 1, so it cannot be ruled out by reducing the tt cross-section uncertainty alone simply due to the lack of raw signal events.
The red and cyan lines are for a heavier neutralino m(χ 0 1 ) = 60 GeV. The improvement on the excluded value of the signal strength does not change as dramatically when the tt cross-section uncertainty is reduced. This occurs because the b-jets are much softer when the neutralinos are heavier, and as such do not pass the selection cuts as often. Intuitively, the size of the tt cross-section uncertainty does not have a large impact when so few stop events pass the selection cuts.
The right panel of Fig. 13 shows the exclusion computed here overlaid on the stop splinter region. The blue shaded region is the same as in Fig. 11. The dashed line shows how the excluded region would change if the uncertainty on the tt cross section were reduced by a factor of 2, while keeping the rest of the analysis fixed. If this uncertainty could be reduced, a large portion of the lower sliver of space not yet excluded by ATLAS could be covered. The dot-dashed line shows the projected exclusion with a factor of 10 reduction of the tt cross-section uncertainty, demonstrating that it is in principle possible to cover the entire lower sliver. While it is unlikely that this approach could ever be used to constrain the upper part of the open region, a reasonable expectation is that complimentary direct searches for stops will play an important role as more data is analyzed.

Discussion
In this work, we have examined the light stop splinter region not yet excluded by ATLAS. We showed the impact of including both finite-width effects and spin correlations. Our work shows that the splinter region can be excluded for m(t 1 ) 180 GeV. We also demonstrated that a reduction in the uncertainty on the top quark production cross section by a factor of 2 would raise the limit on m(t 1 ) by around 20 GeV. This would go a log way toward closing the lower splinter region, leaving the upper part to be closed by dedicated stop searches. In addition, ATLAS typically uses m(χ 0 1 ) = 1 GeV for the lightest point in their exclusion scans. We verified that the limits smoothly extend to m(χ 0 1 ) = 0 GeV. When performing this study, careful attention was needed due to some of the assumptions build into the event generation tools. It is very common when calculating the decay width to assume that the three-body effects are not important if two-body modes are open. While this is often a very good approximation, it is simply not true near the two-body to three-body threshold. Furthermore, the decay widths of the stops in this region are on the order of 0.1-10 MeV; MadGraph recognizes that these widths are smaller than the QCD scale and as a result treats the stops as stable by setting their width to zero. The result is that the stop cannot decay and event generation fails. It does do this for good reason since colored particles, such as the stop, with widths this small could form bound states like stoponium. However, they also decay a large fraction of the time before this happens; see [61][62][63][64][65] for more information and limits on stoponium. A detailed workaround is presented in App. A.
Everywhere in this paper, we assumed the LSP is a bino-like neutralino. It would also be interesting to reinterpret this search assuming the LSP is a gravitino. For much of the parameter space, the stops have small enough widths that they would live long enough to leave the detector. However, if one imagines that the stop width were larger, e.g. due to the presence of another decay mode, it would be interesting to study the impact of the spin correlations for this scenario. Finally, generalizing the exclusion to the m t -m(t 1 )-m(χ 0 1 ) parameter space (as was done in [26]) is not possible within our framework since we do not simulate the top and Standard Model backgrounds. We leave such explorations to future work.
consistently defined compared to the model parameters and not necessarily equal to the PDG values.

B Stealth Stop Event Generation Pitfalls
In Sec. 2, we detailed the need for using correct widths, and showed that common calculators actually get the width wrong near threshold. This forces the user to compute the width and production in completely independent processes. In order to get correct results, the parameters used between these two processes needs to be consistent. Needing to do these separately introduces the extra possibility of user error. This appendix details some of the pitfalls we ran into along the way and should serve as a guide to correct calculations.
In Fig. 14, the importance of using consistent widths is demonstrated. The left panel shows the cross section for stop pair production in the dileptonic channel at √ s = 8 TeV, and the right panel displays the ratio of different methods to our procedure (denoted in black). The results of using the width computed by SDecay are shown in green. As SDecay includes NLO QCD corrections, the width in the denominator of the propagators is larger than expected from the couplings, which leads to a production cross section which is too small. Just above threshold, SDecay only uses the 2-body channel, yielding too small of stop widths and extra large cross sections.
The need to account for more than just the two-body width, along with having to compute the width and production separately, opens the possibility of extra (user) inconsistencies. For instance, in generating the decay width of the stops through to the final state particles, it is tempting to use generate t1 > b n1 all all / susy, expecting that the two "all"s account for the W decay. However, in addition to the W decay products, this also produces diagrams with an un-decayed W and radiated gluons. Because of these extra diagrams, the result of the process is not the width of the stop; the width is much too large. This leads to reduced cross sections, as shown in red in Fig. 14. Alternatively, the cyan lines show the effect of using different widths for the top and W between calculating the stop decay and the production. For this, the compute_widths 6 24 command is used when getting the stop width, but the user "forgets" to update the top and W width from the (slightly larger) default values in the production process. In order to get valid results, the matrix element for the stop width needs to be consistent with that used in the production and subsequent decay.
As a last consistency check, we compare our leptonic cross section with the full stop pair production cross section (p p > t1 t1˜). Using the narrow width approximation (NWA) we take the total cross section and multiply by the leptonic branching ratio (of the W s), and show this in the blue lines of Fig. 14. All of the collider level cuts in the run_card.dat are removed for the full (production and decay) process to get the most accurate comparison. With this, the NWA is larger , while the right panel shows the ratio compared to the method used in this paper. The black lines show the value when the stop width is computed to the final state particles. The blue lines are the narrow width approximation, wherein the total stop production cross section is multiplied by branching ratio for each W to decay leptonically. The green lines correspond to using the width obtained from SDecay, which includes NLO QCD corrections. The red and light blue lines come from common user errors. than our result by ∼ 5 − 6% across all stop masses, consistent with our observations in Sec. 2.2.
To get a further sense of what causes the differences, we show the values of the widths and corresponding cross section for the point m(t 1 ) = 200 GeV in Table 3. In the first block, the widths of the W , t, andt 1 are given, underscored by the method for which they are calculated. The t widths marked by three-body implicitly use the given W width. Similarly, thet 1 widths are calculated all the way to the final state leptons, and use the corresponding W and t widths. In doing so, thet 1 width used in production process is consistent and therefore yields the appropriate cross section. This is exemplified in the row where the W and t widths are set to 1.0, and the stop is calculated consistent with that. The resulting cross section remains the same. In the last block of numbers, the stop width is not consistent.  Table 3: Cross section in pb for different widths in GeV of t, W , andt 1 . We set m t = 173.1 GeV, m(t 1 ) = 200 GeV, and m(χ 0 1 ) = 1 GeV.