Probing Light Stops with Stoponium

We derive new limits on light stops from diboson resonance searches in the $\gamma\gamma$, $Z \gamma$, $ZZ$, $WW$ and $hh$ channels from the first run of the LHC. If the two-body decays of the light stop are mildly suppressed or kinematically forbidden, stoponium bound states will form in $pp$ collisions and subsequently decay via the pair annihilation of the constituent stops to diboson final states, yielding striking resonance signatures. Remarkably, we find that stoponium searches are highly complementary to direct collider searches and indirect probes of light stops such as Higgs coupling measurements. Using an empirical quarkonia potential model and including the first two $S$-wave stoponium states, we find that in the decoupling limit $m_{\widetilde t_1} \lesssim 130$ GeV is excluded for any value of the stop mixing angle and heavy stop mass by the combination of the latest resonance searches and the indirect constraints. The $\gamma \gamma$ searches are the most complementary to the indirect constraints, probing the stop"blind spot"parameter region in which the $h^0 \tilde t_1 \tilde t_1^*$ trilinear coupling is small. Interestingly, we also find that the $Z\gamma$ searches give a stronger constraint, $m_{\widetilde t_1} \lesssim 170$ GeV, if the stop is primarily left-handed. For a scenario with a bino LSP and stop NLSP, several gaps in the direct collider searches for stops can unambiguously be filled with the next run of the LHC. For a stop LSP decaying through an R-parity violating $UDD$ coupling, the stoponium searches can fill the gap 100 GeV $\lesssim m_{\tilde t_1} \lesssim 200$ GeV in the direct searches for couplings $\lambda"\lesssim 10^{-2}$.

Light stops are a quintessential feature of a natural supersymmetric theory [1][2][3][4], being responsible for the cancellation of the dominant Higgs mass quadratic divergence coming from the top quark. However, after Run 1 of the LHC the possible existence of light stops has been strongly constrained by a suite of dedicated searches by ATLAS and CMS [5][6][7][8].
For instance, in a simplified scenario containing a neutralino as the lightest supersymmetric particle (LSP) and the stop as the next-to-lightest supersymmetric particle (NLSP), with a few gaps and caveats, stop masses as high as ∼ 700 GeV and LSP masses as high as ∼ 250 GeV have been excluded [6,8]. Still, as is always the case with direct searches, the limits depend strongly on the spectrum and decays of the stops, and there remain open windows in which a light stop can hide from LHC searches. Examples include compressed [9] or stealth [10] stops and R-parity violating stop decays [11,12]. It is therefore critical to continue to develop new strategies to directly search for light stops [13].
Another opportunity to probe light stops is presented when the stop is long lived, which naturally occurs in a number of motivated scenarios. Indeed, provided the stop does not have an unsuppressed 2-body decay, stop pairs produced through gluon fusion can form a stoponium bound state due to the Coulombic attraction mediated by the strong force [14][15][16][17][18][19][20][21]. Once the bound state is formed, it can decay via annihilation to gg, γγ, Zγ, W W, ZZ and hh, leaving a distinctive resonance signature. As we will emphasize below, the branching ratios to these final states depend essentially on one parameter, the trilinear Higgs-stop-stop (h 0t 1t * constraints on light stops, including vacuum stability and electroweak precision data.
We also discuss the consequences of our results for two motivated scenarios in which stoponium annihilation decays are relevant. The first is the canonical R-parity conserving case with a bino LSP and a stop NLSP. For mt 1 m t + m χ 0 the stoponium annihilation decays are visible, providing a clean and unambiguous probe of light stops which is independent of possible kinematic degeneracies or model dependencies of stop branching ratios. Second, we discuss the case of a stop LSP decaying via the R-parity violating U DD operator to two jets, which is challenging to probe directly due to the large QCD backgrounds. Again, stoponium can provide a clear-cut test of light stops in this scenario.
It is worth emphasizing that there is a potentially large source of theoretical uncertainty present in the stoponium production cross section coming from the assumed potential model and contribution of the excited states to the signal. The limits we present here are based on the empirical quarkonia potential model of Hagiwara et al. [24]. Furthermore, as discussed in Ref. [20], it is likely that higher S-wave stoponium states contribute to the resonance signals, and we therefore include the first two states in our estimate. As we will argue below, we believe our assumptions in this regard are conservative, but we will also discuss in detail how different choices affect the stoponium production cross section and the resulting limits.
We begin in Section II with a brief review on the production and decay of stoponium at hadron colliders. Our new limits and future projections from LHC resonance searches on several benchmark scenarios for each diboson channel are presented in Sec. III. The current limits from stoponium and other indirect constraints on the full stop sector parameter space are combined in Sec. IV. In Sec. V we discuss the implications of our results on two motivated scenarios containing lights stops. Section VI presents our conclusions and outlook. Finally we include several Appendices A-E which collect and summarize our conventions for the stop sector, the indirect constraints on light stops, the stop decay widths, our procedure for limit extrapolation, and the LHC diboson resonance searches used in our analysis.

II. STOPONIUM PRELIMINARIES
In this work we will investigate the annihilation decay signatures of the S-wave (J P C = 0 ++ ) stoponium bound state, ηt, at the LHC. These signatures result from the following processes: pp → ηt + X, ηt → γγ, Zγ, W W, ZZ, hh . . .
Searches for stoponium resonances provide a highly complementary probe of light stops, testing different assumptions about SUSY spectra and parameters than both direct searches for stops and indirect probes such as Higgs coupling measurements, precision electroweak data, and vacuum stability constraints. In particular, direct searches for stops (and other new particles) are inherently model dependent, requiring assumptions about the superpartner spectrum and stop decay chains. In contrast, the limits and projections we will derive for stoponium will depend on two different assumptions: 1) the stop does not have an unsuppressed 2-body decay, and 2) the value of the trilinear Higgs-stop-stop coupling. As we will discuss shortly, the first assumption of a narrow stop width is actually satisfied in a number of interesting scenarios. The second assumption on the value of the h 0t 1t1 coupling, which is determined by the stop sector masses and mixing angle, will govern the stoponium branching ratios for the various final states in Eq. (1).
In this section we will collect the relevant ingredients entering into our study of stoponium.
Our analysis relies on a number of previous works [14][15][16][17], notably those of Martin and collaborators [18][19][20][21]. Next we describe in more detail the basic requirements for observable stoponium annihilation decays. Following this we will review the production of stoponium at the LHC and the stoponium decay branching ratios.
Before beginning, let us note here that we will work exclusively in this paper with the physical stop parameters {mt 1 , mt 2 , θ t }, in whicht 1 (t 2 ) denotes the lightest (heaviest) physical stop state, and a mixing angle θ = 0 (π/2) corresponds to a purely left-handed (righthanded) stop. We will also restrict to the decoupling limit. We refer the reader to Appendix A for our stop sector conventions.

A. Conditions for stoponium formation and annihilation decays
What conditions are required for the distinctive resonance signatures of stoponium to be observable? 1 Provided they are kinematically accessible, stop pairs will be copiously produced in pp collisions. However, whether or not the bound state forms and subsequently decays via the annihilation processes in Eq. (1) depends on the natural width of the stop. The bound state formation time is on the order of the inverse binding energy which, assuming a Coulombic potential, is given by )mt] −1 is the "Bohr radius" of the system. Therefore, stoponium will form provided the natural width of the stop is smaller than the binding energy. Subsequently, the decay of the bound state may proceed either through the prompt decay of the constituent stop or through the annihilation decays. For instance, if the h 0t decay is to a pair of gluons, with partial width (see also Sec. II C below) where R(0) = √ 4πψ(0) is the radial wavefunction at the origin, for which we employ the parameterization in Ref. [18,24]. For mt 1 ∼ O(100) GeV, the factor |R(0)| 2 /m 2 ηt ∼ 0.2 GeV, leading to a partial width Γ(ηt → gg) ∼ O(1) MeV in the light stop mass range of interest.
If the annihilation decay width dominates over the natural stop width we can potentially observe the stoponium through its resonance signature.
In summary we require, In practice, since the annihilation decay width is typically smaller than the binding energy, it sets the upper bound in Eq. (3). In fact, these conditions are satisfied in a number of interesting of scenarios and over a wide range of model parameters, particularly when the stops are light. We will examine two such scenarios in Section V.

B. Production
At leading order (LO) the production of stoponium in pp collisions at the LHC proceeds via gluon-gluon fusion, with a cross section given by [18] where P gg is the gluon parton luminosity as a function of τ ≡ m ηt /s, with s the squared center-of-mass energy, defined in Eq. (D.2) in Appendix D. The partial decay width for Γ(ηt → gg) is given in Eq. (2). Next-to-leading order (NLO) perturbative QCD corrections to stoponium production have been computed in Ref. [20] and with the threshold resummation in Ref. [17]. To account for these corrections we will use the results of Ref. [20]. In particular, we extract the cross section prediction from their Fig. 9 and furthermore estimate a conservative ∼ 10% scale uncertainty on this prediction from their Fig. 7 [20]. We can also reasonably estimate a ∼ 5 − 10% PDF uncertainty in the mass range of interest, similarly to that for Higgs production via gluon-gluon fusion [26].
There are likely larger sources of theoretical uncertainties than the QCD scale uncertainty, which are, however, difficult to quantify. Notice in particular that the cross section in Eq. (4) depends on the value of the stoponium wavefunction at the origin squared, |R(0)| 2 , through Eq. (2), and our imprecise understanding of this matrix element potentially provides a significant source of uncertainty. Ref. [20] employs the Λ (4) MS = 300 MeV parameterization of the wavefunction at the origin from the study of Hagiwara et al. [24]. This parameterization is based on an empirical quarkonia potential model, which deviates significantly from the pure Coulombic form at large mt 1 due to the effect of higher order QCD corrections. For instance, the prediction for |R(0)| 2 for the ground state is smaller by roughly a factor of 2 in comparison to that obtained from a Coulomb potential, suggesting that our estimates of the production rate are conservative. Ref. [24] also present Λ MeV parameterizations, which lead to a ∼ 10% reduction or enhancement in the stoponium production rate, respectively, giving some sense of the uncertainties involved. We note that the study of Ref. [24] is 25 years old, and given the significant experimental and theoretical progress in the studies of quarkonia and QCD in the interim, it would be extremely useful to perform a modern analysis in the spirit of [24], with the particular aim of assessing the theoretical uncertainties.
Another source of uncertainty in the prediction is the contribution of higher S-wave stoponium states to the annihilation decay signals. Ref. [20] includes the first two states in their calculation, but reasonably argues that higher states will likely add to the signal, and in this way our estimates of the production rate are likely conservative. In Figure 1 we display the stoponium production cross section at the LHC under several different assumptions regarding the potential and contribution of excited states. Needless to say, it would be extremely interesting to study the stoponium system in more detail to both quantify the uncertainty due to the modelling of the bound state potential as well as understand the phenomenological consequences of the higher radial states.

C. Decay branching ratios
For the stoponium decay branching ratios, we use the LO results presented in Ref. [18] for all decay channels except for the ηt → γγ channel. For the diphoton channel we obtain the partial width Γ(ηt → γγ) by multiplying the LO partial width Γ LO (ηt → gg) [18] by the NLO ratio Γ NLO (ηt → γγ)/Γ NLO (ηt → gg) computed in Ref. [19]. This approximation reproduces the correct NLO BR(γγ) when BR(gg) dominates, and it is this case in which the diphoton channel is most important 2 . When other decay modes such as hh or W W dominate instead, the diphoton branching ratio is anyway too small to induce an observable signal. We use the two-loop running α S (µ) and the one-loop running α(µ) renormalized [��] FIG. 1: Stoponium production cross sections at the LHC. Our baseline cross section for √ s = 8 TeV (solid blue) and √ s = 14 TeV (solid green) from Ref. [20], which assume the Λ (4) MS = 300 MeV parameterization of the stoponium wavefunction at the origin from Ref. [24] and include the contribution of the first two S-wave states. These cross sections are used to derive our limits and projections. For √ s = 8 TeV we also display for comparison the cross section derived with the quarkonium potential [24] including instead only the ground state (blue dashed) and the first 10 S-wave states (blue dotted), as well as the cross section assuming a Coulomb potential (red) (including only the ground state) with α S evaluated at µ = 1/ r 1S = 2 α S (µ) m ηt /9 . from α S (m Z ) = 0.118 and α(m Z ) = 1/128.0.
We will work in the decoupling limit. The minimal low energy spectrum will contain the stops and possibly the left-handed sbottom. It is also possible that the spectrum contains a bino-like neutralino LSP, since the condition (3) is readily satisfied when mt 1 m t +m χ 0 (see Sec. V for a detailed discussion of this scenario), although the stoponium annihilation decays to neutralinos in this case are always negligible. We will not consider Higgsinos or Winos that are lighter thant 1 , as in this case we expect the stop to have an unsuppressed two body decay to a lighter chargino and a bottom quark, generically violating the condition (3) 3 All other heavy Higgs scalars and superpartners are taken heavy and decoupled and will not influence the stoponium decays. The mass of the left-handed sbottom is fixed by assuming vanishing sbottom mixing (see Eq. (B.8)). We emphasize that left-handed sbottoms can have nonnegligible effects on the stoponium branching ratios when the lightest stop is primarily lefthanded. As typically happens, unitarity and gauge invariance induce a (partial) cancellation /v (upper) and stoponium branching ratios (lower) for m t 1 = 120 GeV (left) and 160 GeV (right) as functions of the stop mixing angle. A strong correlation between λ h 0t 1t * 1 and the stoponium branching ratios is observed. The shaded regions are constrained indirectly by Higgs signal strength data, precision electroweak measurements, and vacuum stability as discussed in Appendix B. Here we have fixed m t 2 = 1 TeV. Note that in the left panel we have terminated the plot at small θ t where the sbottom becomes tachyonic.
in the W W partial width between t-channel stop and sbottom exchanges if the stop is mostly left-handed [15]. Finally, the stoponium decays to bb and tt are also subdominant and will not be considered further.
The main decay modes, ηt → gg, γγ, Zγ, W W, ZZ and hh, are described well by the three physical stop parameters, m t 1 , m t 2 , θ t , and t β ≡ tan β. These parameters are defined in Appendix A. The dependence on t β is weak in most cases, so we fix t β = 10 throughout in this paper. The trilinear coupling of the 125 GeV Higgs to the light stops, λ h 0t 1t * 1 , plays a crucial role in determining the stoponium decay pattern. This coupling is written in terms of the physical stop parameters as Due to the last term in Eq. (5), this coupling can become much larger than the electroweak scale, v = 174 GeV, when the stop mixing is large, inducing sizable partial decay widths for stoponium to W W, ZZ, hh via the s-channel exchange of h 0 . This is illustrated in Fig. 2, where we display the mixing angle dependence of the coupling λ h 0t 1t * 1 and the stoponium branching ratios. A strong correlation is observed between the value of λ h 0t 1t * 1 and the pattern of branching ratios. For large values of mixing, the coupling λ h 0t 1t * 1 grows and the stoponium dominantly decays to W W , ZZ, and, if kinematically allowed, hh. Instead, when λ h 0t 1t * 1 is small, the gg decay dominates and, importantly, the γγ and Zγ branching ratios are maximized (see Sec. III A and Sec. III B for further discussion).
Large values of λ h 0t 1t * 1 are constrained, particularly for light stops, by indirect tests including Higgs signal strength measurements, electroweak precision data, and vacuum stability. These indirect constraints are reviewed in Appendix B and shown as the shaded regions in Fig. 2. In particular, we observe that only a small region around λ h 0t 1t * 1 = 0 is consistent with these indirect constraints for very light stops. It is therefore very interesting to consider the ηt → γγ (and Zγ) channel as it can provide an independent probe of this unconstrained region of parameter space.

III. NEW LIMITS ON LIGHT STOPS AND FUTURE PROSPECTS
We now derive the current limits on light stops from stoponium annihilation decays and estimate the future reach of each diboson channel in Eq. (1): γγ, Zγ, W W, ZZ and hh.
In this section, we simply assume that the stoponium forms and decays via annihilation, but we will discuss particular models and relevant parameter space in Sec. V. For each diboson channel, we define and consider the ideal or benchmark branching ratio that best represents the optimal or characteristic reach, respectively, over a wide range of the stop masses. Based on these benchmarks, we present the limits as a function of the stoponium mass. Next, in Sec. IV we exhibit the exclusion limits in the general stop parameter space, i.e., as a function of the stop masses and mixing angles, accounting for the realistic branching ratios and comparing with the indirect constraints described in Appendix B.
Our limits and projections are obtained by recasting the results of the latest heavy resonance searches at LHC 7+8 TeV, which are collected in Appendix E. We will present the strongest limit in each diboson channel. For the future sensitivities, we extrapolate the current expected limits to LHC 14 TeV with 30 fb −1 and 3 ab −1 . Our extrapolation method and assumptions on statistical and systematic uncertainties are described in Appendix D.

A. γγ
The diphoton channel has long been recognized as a promising place to search for stoponium [15]. At LO, the ηt → gg and ηt → γγ partial decay widths do not depend on stop sector parameters other than the lightest stop mass, m t 1 . They are induced solely by the strong and electromagnetic interactions, such that the ratio, is fixed by gauge couplings and independent of m t 1 . As discussed in Sec. II C, in our numerics we include the NLO corrections to the ratio in Eq. (6) from Ref. [19,20]. These corrections, along with the running gauge couplings, induce a mild dependence on the lightest stop mass, with the ratio (6) varying from 0.0025 − 0.0042 as the lightest stop mass mt 1 varies from 100 -400 GeV. The NLO corrections actually reduce the ratio (6) by 35%-25% from the LO result [19].
Since the gg and γγ decay modes are always present and the size of the corresponding partial decay widths are rather model-independent, the diphoton branching ratio is maximized simply when the other possible decay channels (W W , ZZ, hh) are minimized. We will therefore define the ideal γγ scenario by the idealized limit in which the branching ratio BR(ηt → γγ) is given by the ratio in Eq. (6) (including NLO corrections as discussed above).
We note that the idealized branching ratios are 10 2 -10 5 times bigger than that of the SM Higgs boson with the same mass in the range of m η t = 200 − 600 GeV.
In the full stop parameter space, the ideal diphoton scenario is approximately realized in the region where the Higgs-stop-stop coupling λ h 0 t 1 t * 1 is small, corresponding to mixing angles In this region of parameter space the ηt → W W, ZZ, hh diagrams involving s-channel Higgs exchange are of a similar size or smaller than those mediated by the weak gauge couplings, thus explaining the dominance of the ηt → gg mode. This can also be seen from the dependence of λ h 0 t 1 t * 1 and the stoponium branching ratios on the stop mixing angle in Fig. 2, where we observe that the gg branching ratio is nearly 100% when the λ h 0t 1t * 1 is small. For the stop mass mt 1 = 100 − 300 GeV (mt 2 = 1 TeV), the ideal γγ scenario is approximately realized for θ t = 0.10 − 0.22 or θ t = 1.32 − 1.48. The actual maximum diphoton branching ratio in these parameter regions is smaller than the ideal one of Eq. (6) only by about 0.7−1.8% fraction. In the next section, we will translate our limits to the full stop parameter space, accounting for the realistic branching ratios.
We note that the ideal γγ scenario lies within the blind spot region consistent with current indirect constraints, even for very light stop masses of order 100 GeV. In particular, the constraints from the Higgs signal strength data are relaxed in this region because the coupling λ h 0 t 1 t * 1 is small and the modification of the gluon-gluon fusion cross section from stop loops is minimized. This makes the ηt → γγ channel particularly interesting to study as it can provide a complementary probe of this parameter region.
In Fig. 3, we show the current limits from ATLAS and CMS diphoton resonance searches and the future sensitivities for the ideal diphoton scenario. Stoponium masses m ηt ≈ 2mt 1 in the the vertical blue bands are excluded by the observed limits. We note that the ATLAS and CMS observed limits display complementary patterns of statistical fluctuations, excluding somewhat different stop mass ranges. The union of the two searches excludes stop masses up to about m t 1 ∼ 136 GeV. We emphasize that these exclusions are based on the assumption of including the first two S-wave stoponium states in the production rate, as described in Sec. II B. This is a conservative assumption as it is likely that higher modes also contribute to the diphoton signal. Following Ref. [20], for comparison we have also displayed in Fig. 3 the cross section including only first S-wave mode or first 10 S-wave modes, allowing the reader to see the effect of this uncertainty. The reader can also see the difference between the two potential models: the thick blue line indicates the empirical quarkonia potential in Ref. [24] and the dot-dashed red line represents the Coulomb potential.
We also show in Fig. 3

B. Zγ
Another important channel is Zγ since the branching ratio is also maximized when the coupling λ h 0t 1t * 1 is small and the indirect constraints are absent. As for the diphoton channel, we can write the ratio of Zγ to gg partial widths at LO: As in the diphoton case, the BR(Zγ) is maximized when the BR(gg) dominates and the other decay modes (ZZ, W W, hh) are suppressed, corresponding to small λ h 0t 1t * 1 . However, unlike the diphoton case the mixing angle dependence in Eq. (8) implies that the BR(Zγ) is largest when the lightest stopt 1 is mostly left-handed. This asymmetry is also clearly visible in Fig. 2. We therefore define the ideal Zγ scenario by evaluating BR(Zγ) with stop parameters corresponding to λ h 0t 1t * 1 = 0 and a mostly left-handedt 1 . In practice, we choose different mixing angles satisfying this condition for different stop masses. The mixing angle and the BR(Zγ) for the ideal Zγ scenario gradually increase from θ t = 0.169 to 0.184 and from 0.84% to 1.37%, respectively, as the stop mass varies in the range m t 1 = 100 − 400 GeV. We note that these branching ratios are 100 -3000 times larger than those of the SM Higgs boson with the same mass in the range of m η t = 200 − 600 GeV.
We show the current limits on the ideal Zγ scenario from the Zγ → + − γ resonance searches in the upper panels of Fig. 4  can be sensitive to 300 GeV with 3/ab. As we will see in Fig. 8 in the next section, the diphoton channel is most complementary to the indirect constraints, as it is sensitive to the entire small λ h 0t 1t * 1 region including primarily right-handed stops. On the other hand, the Zγ channel provides stronger limits than the diphoton channel for mostly left-handed stops.
It is important to mention that our Zγ results are based on the LO calculation while the γγ result in the previous subsection is from the NLO calculation in Ref. [19,20]. For the latter case, the NLO effects reduce the ideal γγ branching ratio by 35%-25% for the stop mass in the range 100-400 GeV. If the same NLO reduction is applied to the ratio of the Zγ to gg partial widths (hence, to the ideal BR(Zγ)), the expected exclusion on the stop mass is approximately reduced by about 20 GeV. In addition to the NLO perturbative QCD effects, uncertainties from the number of S-wave stoponium states contributing to the signal and from potential models are again significant as shown by multiple theoretical predictions in Fig. 4.

C. ZZ and W W
The branching ratios BR(ZZ) and BR(W W ) are intimately related with the coupling λ h 0t 1t * 1 as discussed in Sec. II C. The indirect probes of light stops discussed in Sec. IV place strong constraints on large values of this coupling and in turn determine the maximum phenomenologically allowed branching ratios BR(ZZ) and BR(W W ), as is depicted in Fig. 2.
Unlike the diphoton channel discussed above, there are several reasonable ways to define the benchmark ZZ, W W scenarios. For instance, one could simply take the largest theoretically allowed branching ratio, or alternatively the largest branching ratio consistent with indirect constraints; in both cases one is led to vary the mixing angle as a function of the stop masses. One could also choose a fixed mixing angle for a range of the stop mass. For simplicity, we will instead choose a prescription with a fixed branching ratio that is broadly consistent with the current indirect constraints. Referring to Fig. 2 of m t 1 = 160 GeV case, we define the benchmark ZZ, W W scenarios by choosing BR(ZZ)=6% and BR(W W )=12%. In fact, for the range of θ t allowed for the 160 GeV stop we obtain maximum BR(ZZ) 4 -8% and BR(W W ) 11 -16% for the stop mass m t 1 = 100 − 300 GeV (with m t 2 = 1 TeV). Thus, the maximum branching ratios in this range of θ t do not vary much with the stop mass. The current indirect constraints are somewhat weaker for heavier stops and a wider range of θ t is viable, but we envisage that indirect constraints will also improve with more data, probing a broader range of mixing angles and heavier stop mass.
We refer to the mt 1 = 160 GeV case as a representative benchmark scenario in view of the interplay between stoponium and indirect constraints. Again, we provide limits on the full stop parameter space accounting for true branching ratios in the next section. Using

D. hh
Similarly to the ZZ and W W branching ratios, the maximum branching ratio Br(hh) is limited by indirect constraints on the coupling λ h 0t 1t * 1 . We define the benchmark hh scenario by a fixed BR(hh)=0.25. The actual maximal BR(hh) for the range of θ t currently allowed for the 160 GeV stop gradually increases from 16% to 41% for the 130 GeV to 300 GeV stop.
In Fig. 7, we display the current sensitivity and 14 TeV projections on the benchmark hh scenario from the hh → bbγγ resonance search. No new limits can be obtained with this channel from the Run 1 data. However, in the long term this channel is expected to probe stops as heavy as 215 GeV in the benchmark hh scenario at √ s = 14 TeV with 3/ab. Assuming the stoponium forms and annihilate decays, we observe that light stops with masses mt 1 130 GeV are excluded by the combination of diboson resonance searches and the indirect constraints. Stoponium searches alone probe stops lighter than about 125 GeV. Remarkably, a strong complementarity is seen between the γγ resonance searches and the indirect constraints (particularly the Higgs coupling constraints), with the γγ channel probing the remaining open stop blind spot region corresponding to small or vanishing coupling λ h 0t 1t * 1 . Moreover, we also see that primarily left-handed stops with masses mt 1 170 GeV are constrained by Zγ resonance searches.
Again, we emphasize that the limits we present here rely on the Λ (4) MS = 300 MeV parameterization of the stoponium wavefunction at the origin from Ref. [24] and the inclusion of the first two S-wave states [20]. While it is likely that this estimate is conservative, particularly in light of the possible additional signal coming from the higher excited states, it is worthwhile to examine how the limits change if different assumptions are made. For example, if we instead include only the contribution of the 1S state, the expected limit on the stoponium mass will be weaker by about 20 GeV, as can be seen by examining Figs. 3  and 4. Alternatively, if one assumes the pure Coulomb potential, the expected limit on the stoponium mass will be stronger by about 70 GeV. In any case, it is clear that stoponium is now probing the hypothesis of light stops beyond the LEP limits and will explore uncharted territory with the next run of the LHC.

V. EXAMPLE SCENARIOS OF INTEREST FOR STOPONIUM
We now discuss two motivated examples scenarios in which stoponium provides a complementary probe to direct stop searches and indirect tests. The condition (3) will generically hold when the stop does not possess an unsupressed 2-body decay. This can happen due to kinematics or small couplings, and to demonstrate this we will focus two examples -1) a "canonical" case of R-parity conservation with a bino LSP and a stop NLSP, and 2) a R-partiy violating case with a stop LSP which decays to a pair of jets through the U DD operator.

A. R-Parity conserving bino LSP and stop NLSP
A canonical benchmark scenario with light stops consists of a stop NLSP and a neutralino LSP, under the assumption that R-parity is conserved. In this case the neutralino is stable, leading to the signature of missing energy at colliders. As is well known, in this scenario the stop width is very sensitive to the mass spectrum, and several kinematic regions can be distinguished: I. mt 1 > m t + m χ 0 : The stop decays via the 2-body processt 1 → tχ 0 .
To illustrate these features, in the left panel of Fig. 9 we have plotted the decay width of the stop as a function of its mass for the case of a 10 GeV bino LSP (see the Appendix C for the 3-body partial decay width of the stop). For comparison, in the same plot we show the stoponium binding energy and annihilation decay width. As the stop mass approaches the boundary between 2-and 3-body decays, we observe that Γ η t becomes larger than the natural stop width. Therefore, in region I above, the stop generally decays too quickly for stoponium resonance signatures to be visible, although it may be possible near the edge of region I where phase space suppression is relevant. However, in regions II and III the dominant stop decay channels are 3-and 4-body (or suppressed 2-body), respectively, allowing both sufficient time for the formation of the stoponium and the dominance of the bound state annihilation decays. Therefore, regions II and III are prime targets for resonance searches for stoponium, which is clearly illustrated by the "phase" diagram in the stop -neutralino mass plane shown in the right panel of Fig. 9. We represent the current exclusion limit from stoponium resonances by the blue hatched region; this region is excluded for the given lightest stop mass for all other stop parameters (see Fig. 8). For comparison, we have overlaid the current ATLAS limits from direct stop searches [6], represented by the orange dotted line. While these direct limits naively appear to cover most of light stop region, it is important to emphasize that they rely on the assumption of a 100% branching ratio to the final state under consideration. If this assumption is relaxed, the limits are weakened considerably. For example, in region III if the stop has comparable branching ratios in the 4-body bf f χ 0 and 2-body cχ 0 channels, a significant portion of parameter space opens up [28], and only the orange-shaded region can be excluded regardless of the assumption on the branching ratio. In contrast, constraints from stoponium resonances are already starting to unambiguously probe this region.
It is also worthwhile to note that the "stealth" stop region -the phase space-suppressed part of region I near m t 1 ∼ m t , m χ 0 0 -is also potentially amenable to stoponium searches (within the red-shaded region in the bottom part of Fig. 9) This region is challenging to probe with direct stop searches, and several dedicated observables such as spin correlations [29,30], top pair production rates [31,32] and special kinematic variables [33] have been studied.
However, they primarily rely on precision measurements and calculations, which are often subject to subtle uncertainties. On the other hand, the stoponium resonance searches are cleaner and less ambiguous in most cases. Thus, the stoponium can provide an important . For comparison, we also display a union of ATLAS limits on the stops (orange dotted), which assume 100% branching ratios to the final states under consideration [6] as well as a branching ratio-independent exclusion limit taken from Ref. [28] (orange shaded). The blue-hatched region is excluded from stoponium Fig. 8. complementary probe of the stealth region. In this region, stops with a mostly left-handed mixture (small θ t ) can be probed via stoponium over a wider range of stop masses (see Fig. 9 right) because the individual stop decay width is smaller with a smaller θ t . In any case, stoponium does not yet have sensitivity to stealth stops (the blue-hatched region extends to about 130 GeV), but can potentially be sensitive to a significant portion of this challenging region in the early part of Run 2 with ∼ O(10)/fb of data (see Fig. 3).
Let us finally comment on the case in which the LSP is a very weakly interacting particle such as the gravitino. In this scenario, the stop NLSP is typically long lived due to its weak coupling to the gravitino, and therefore stoponium annihilation decays can be relevant.
For instance, taking a low supersymmetry breaking scale √ F = 10 TeV (corresponding to an essentially massless gravitino), a 200 GeV stop has a decay width Γ( t 1 → t G) ∼ 10 −9 GeV [34], which is much smaller than the typical stoponium annihilation width. In practice, however, this scenario is already strongly constrained by direct searches covering for typical values of stoponium binding energies and annihilation decay widths. We also display a contour of the stop lifetime cτt 1 = 0.1 mm. Stops with lifetimes in this range and longer are subject to stringent LHC searches for long-lived charged particles [35,41]. Stops with prompt RPV decays in the red-shaded region are also constrained from RPV searches [36]. Stoponium can therefore be relevant for a broad range of RPV couplings, λ 10 −2 , and the current exclusion from stoponium resonances is shown as blue-hatched.
both prompt and displaced decays [35], and only a small window around the stealth stop parameter region for prompt decays is still open. The stoponium searches will be able to fill this open region.

B. RPV stop LSP
Another benchmark scenario of interest is a stop LSP that decays to a pair of jets due to a small U DD R-parity violating coupling, which we customarily denote as λ . Such a RPV stop is challenging to probe with direct searches at LHC. Despite its large production cross section, the signature of four jets (paired dijet resonances) is difficult to disentangle from the QCD background, and stops as light as ∼ 100 GeV are still allowed [11,12,36] although some range of the stop masses 200 − 350 GeV has recently been excluded [37].
Stoponium resonance searches provide another means of probing this interesting scenario.
Since the decay of the stop is through a 2-body process, the RPV coupling must be somewhat suppressed in order for the condition (3) to apply. In fact, this may be motivated in explicit models of R-parity violation which aim to protect the lifetime of the proton through other symmetries. For example, in the scenario of MFV SUSY [38,39], the dominant RPV coupling is λ tsb ∼ O(10 −4 ), which is small enough for (3) to hold. Note also that if the RPV coupling is extremely small, λ O(10 −7 ) then the stop will have a macroscopic lifetime. In this case, there are additional strong constraints from searches for displaced dijets [35,40] or long-lived colored/charged particles [35,41]. In Fig. 10 we display the RPV stop LSP phase diagram. We conclude that stoponium resonance searches can provide new complementary probe of this scenario for a broad range of couplings, 10 −7 < λ < 10 −2 .

VI. OUTLOOK
A natural supersymmetric solution to the hierarchy problem predicts light stops, which can lead to a variety of novel phenomena. Besides the suite of direct searches at the LHC (which depend on the SUSY spectrum) and indirect tests such as Higgs couplings, precision electroweak data and vacuum stability requirements, it is possible that light stops can manifest through a stoponium annihilation decay leading to a resonance signature in the γγ, Zγ, ZZ, W W and hh channels. This will happen if the stop does not have an unsuppressed two body decay and naturally occurs in several motivated scenarios.
In this work we have derived new limits on light stops from ATLAS and CMS diboson resonance searches using Run 1 data. Our limits are derived using the empirical quarkonia potential model of [24] and assuming the first two S-wave states contributing to the resonance signal, which we have argued is conservative. Assuming the stoponium can form and annihilate decay, light stops below about 130 GeV are excluded by the combination of constraints on stoponium resonances and other indirect constraints, such as Higgs couplings measurements. Notably, we have demonstrated that the γγ channel is especially complementary to these indirect probes, as the resonance signature in this channel is maximized in the blind spot region where the Higgs-stop-stop trilinear coupling is small and these constraints are weakest. Furthermore, Zγ resonance searches provide even stronger limits for primarily left handed stops, limiting stop masses below about 170 GeV. In the long term, the LHC experiments will be able to probe stops in the 300-400 GeV range with a high-luminosity run with 3/ab at 14 TeV.
We have also discussed the implications of these searches for two specific SUSY scenarios with lights stops. In the case of a bino LSP and stop NLSP, stoponium annihilation decays can be relevant since the stop width is naturally suppressed over a wide range of parameters.
Stoponium resonance searches provide a relatively clean and unambiguous probe, and can give a handle on the challenging compressed and stealth stop regions of parameter space.
Instead for the case of a R-parity violating stop LSP decaying to two jets, stoponium already provides unique sensitivity to stops lighter than about 130 GeV and RPV couplings in the range 10 −7 < λ < 10 −2 , which is difficult to probe directly due to the large QCD backgrounds.
Looking forward, it is clearly of interest on the theoretical side to bring under better control the uncertainties coming from our imprecise understanding of the stoponium bound state system. A modern study of empirical quarkonia potential models along the lines of Ref. [24] is warranted, and perhaps lattice studies could provide further insight into the non-perturbative matrix elements entering into stoponium production and decays. A more detailed investigation into the decay patterns of the excited states would also be helpful in order to understand their contribution to the resonance signals and perhaps uncover additional signatures of stoponium. On the experimental side, stoponium clearly provides a well motivated target for resonance searches, which should be a high priority during the next run of the LHC. In terms of the gauge eigenstates (t L ,t R ), the stop mass matrix is given at tree level by where m 2 Q 3 , m 2 U 3 are the left-and right-handed squark soft mass parameters, X t ≡ A t − µ/ tan β, with A t the soft trilinear coupling, µ the supersymmetric Higgs mass parameter, and tan β the ratio of up and down type Higgs vacuum expectation values, and D t L = m 2 Z cos 2β 1 2 − 2 3 s 2 W , D t R = m 2 Z cos 2β 2 3 s 2 W . For simplicity, we will assume all parameters are real. The physical stop mass eigenstates are related to the gauge eigenstates through the orthogonal transformation: where the mixing angle θ t satisfies which we take in the range [0, π/2]. The sbottom sector is described in a similar fashion. The sbottom mass matrix can be obtained from Eq. A.1 with the replacements

Appendix B: Indirect constraints on light stops
In this appendix we summarize the existing indirect constraints on light stops coming from Higgs signal strength measurements, precision electroweak data, and vacuum stability. In the MSSM, there will be additional constraints on the stop parameters coming from the requirement of obtaining a 125 GeV Higgs mass, but we will remain open to new contributions to the Higgs quartic coming from physics beyond the MSSM. Furthermore, if Higgsinos or gluinos are light there can be further constraints from flavor physics, which we will not include here. Several summary plots of these constraints are shown in Fig. 11.

Higgs signal strength data
Light stops modify the hgg and hγγ couplings at one loop and are therefore subject to constraints by ATLAS and CMS measurements of the Higgs signal strength parameters [42].
To derive the exclusions, we use the method of Ref. [43], with the recent updates in Ref. [44], which include the latest Run 1 results as of summer 2014. In this approach, the raw signal strength data are combined to derive a ∆χ 2 function for the γγ, V V , and b/τ channels as a function of two combined production mode signal strengths: 1) gluon-gluon fusion and tth (ggF+ttH), and 2) vector boson fusion and associated production (VBF+VH). This yields eight composite signal strengths,μ (ggH+ttH,VBF+VH) (γγ,V V,bb,ττ ) . In terms of the ratios r i ≡ Γ(h → i)/Γ(h → i) SM , the signal strength predictions are given bŷ One loop sbottom and stop exchange lead to the following expressions for r gg , r γγ : where C(r) = 1/2, d(r) = 3, Qt = 2/3, Qb = −1/3, and A SM gg ≈ 1.3, A SM γγ ≈ 6.6. The Higgs-squark trilinear couplings λ i can be found in, e.g., Ref. [45]. The scalar loop function A 0 (m) is defined in Ref. [46] and approaches A 0 (m) → 1/3 in the limit m m h 0 .
With these ingredients we construct the ∆χ 2 function. As r gg and r γγ are correlated, we define a one parameter best-fit region in the physical stop parameter space at 2σ C.L. by demanding ∆χ 2 < 4.

Electroweak precision data
Regarding the precision electroweak data, the largest effect of the stops and sbottoms is to provide a new one-loop contribution to the ρ parameter [47]: √ 2v 2 and we have defined the function We have used the left-handed sbottom mass, m b L , in Eq. (B.8) with the zero sbottom mixing.
We apply the constraint from Ref. [ We have assumed that the effects of H d are negligible, appropriate in the large t β regime, and have substituted s β → 1, c β → 0, and H u → h/ √ 2. To have a proper electroweak vacuum and the measured Higgs mass, we set m 2 H = − 1 2 m 2 h , m h = 125 GeV, and δλ = m 2 h 2v 2 − g 2 8 − g 2 8 . In particular, we assume the presence of some physics, either radiative stop corrections in the MSSM or physics beyond the MSSM, that generates an appropriate correction to the Higgs quartic coupling, δλ. To arrive at the above form of the g S terms, we assume that the stops are aligned or anti-aligned in the SU (3) C color space, which maximizes the vacuum stability constraint [50].
The negative X t -linear term can induce new minima when all three scalar fields obtain vevs. We numerically scan the field space to determine all local minima, and we use the CosmoTransitions program [51] to calculate the tunnelling rate between our electroweak minimum and any deeper charge-color breaking minima. The classical action for this tunneling, S E ≥ 400 is required for the metastability of our vacuum. We find that the parameter space with large values of the Higgs-stop-stop coupling, λ h 0t 1t * 1 , is constrained from the requirement of vacuum (meta)stability.
When the left-handed stop is light, we also check whether the mass of the accompanying left-handed sbottom is positive (and evading LEP bounds). For the given stop parameters and heavy enough right-handed sbottoms, the left-handed sbottom is the heaviest when the sbottom mixing is zero (the equality below holds in this case) Throughout in this work, we use this left-handed sbottom mass with the equality sign by assuming the zero sbottom mixing.
The partial decay width fort → t * χ 0 → W bχ 0 can be written as where we have defined the "off-shell" partial widths: Here, we have defined the couplings, with N ij the neutralino mixing matrix elements, and the kinematic function λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc. Notice that the three-body partial decay width reduces to the two-body one in the appropriate kinematical regime through the use of the narrow-width approximation. FIG. 11: Indirect constraints on light stops: Here we display for light stop masses mt 1 = (100,130,160,200,300) GeV the constraints in the θ t -mt 2 plane coming from Higgs signal strength measurements (gray), precision electroweak data (orange), vacuum stability (yellow), tachyonic sbottoms (blue) and sbottoms below the LEP kinematic reach (light blue). The white regions are currently unconstrained. The red dashed line indicates where the Higgs-stop-stop coupling λ h 0t 1t * 1 vanishes. Here we have fixed tan β = 10 and the sbottom mixing to zero (X b = 0).

Appendix D: Limit Extrapolation to Future Searches
To estimate the future prospects for probing stoponium at the LHC we extrapolate the current expected limits of the diboson resonance searches to √ s = 14 TeV with integrated luminosities L of 30/fb and 3/ab. To accomplish this, we employ several simplifying assumptions. For the uncertainties we consider two extreme cases: (1) statistical uncertainties dominate, and (2) systematic uncertainties dominate and improve in proportion to √ L. The results in the main text of Section III are obtained with the former assumption, but we will offer a comparison between the two cases below.
Statistical uncertainty dominant: At a specified collision energy and integrated luminosity denoted collectively by a superscript i, we have where σ i bound is the upper bound on the signal cross-section, σ i B is the background crosssection, and N i B is the number of background events after all cuts. The latter two quantities are related as N i B = σ i B L i i =σ B P i qq L i i , whereσ B is the partonic background cross-section (which does not depend on i), L i is the integrated luminosity, i is the selection efficiency, and P i ab is the parton luminosity for initial partons ab, defined as In Eq. (D.1) we have taken P i qq since the dominant backgrounds in the most sensitive channels, γγ and ZZ, arise from qq initial states. In our numerics we take the up quark parton luminosity for simplicity, as it will dominate in pp collisions.
We furthermore assume a constant efficiency, i = , which is reasonable since the signalto-background ratio near a resonance is largely determined by the resonance mass. This is similar to gluino pair searches [52] and perhaps more generally useful [53].
Taking ratios, we can extrapolate the current bounds to those for other collision energies and luminosities in a simple way in terms of parton luminosity ratio and luminosity ratio Again taking ratios, we are able to extrapolate the current limits: Comparing Eq. (D.3) and Eq. (D.5), we see that if errors are dominantly from systematics, one obtains a weaker cross-section bound by a factor P i qq /P j qq . Numerically, this equals to 1.45 and 1.73 for m t 1 = 160, 400 GeV at 14 TeV collision for q = u. These results roughly show how sensitively the future bounds may depend on the assumptions on errors.

Appendix E: Resonance Searches
In Table I, we collect the latest experimental results on diboson resonance searches. The searches shown in bold give the strongest constraints in each diboson channel and are used to derive limits on stoponium. The CMS and ATLAS γγ limits are similar but constrain complementary mass ranges due to different statistical fluctuations, and we therefore show both results. The CMS Zγ is found to be somewhat stronger than that of ATLAS, but due to the importance of this channel, we show both results and encourage more careful analysis. All searches use 8 TeV datasets except for the ones with ∼24/fb and ∼10/fb, in which case the combined 7+8 TeV dataset is used. The CMS ZZ → 4 result is stronger than the ATLAS one partly because CMS presents combined results of all production channels and uses the 7+8 TeV dataset. The W W → ν2j result has similar sensitivity to that of the W W → 2 2ν, but we show only the latter result for simplicity. The searches listed in the last four lines give weaker constraints. Lastly, ATLAS γγ and Zγ present limits on the fiducial cross-section, so we assume and unfold a constant fiducial efficiency for each channel as mentioned in the relevant part of text.