A survey for low stau yields in the MSSM

We study the implications of LHC results for the abundance of long-lived staus after freeze-out from thermal equilibrium in a super-WIMP dark matter scenario. We classify regions in the MSSM parameter space according to the stau yield, considering all possible co-annihilation effects as well as the effects of resonances and large Higgs-sfermion couplings. Afterwards, we examine the viability of these regions after imposing experimental and theoretical constraints, in particular a Higgs mass around 125 GeV and null-searches for heavy stable charged particles (HSCP) at the LHC. We work in a pMSSM framework and perform a Monte Carlo scan over the parameter space. To interpret the HSCP searches in our scenario, we consider all potentially important superparticle production processes, developing a fast estimator for NLO cross sections for electroweak and strong production at the LHC. After applying all constraints, we find that stau yields below 10^-14 occur only for resonant annihilation via a heavy Higgs in combination with either co-annihilation or large left-right stau mixing. We encounter allowed points with yields as low as 2x10^-16, thus satisfying limits from big bang nucleosynthesis even for large stau lifetimes.


Introduction
One way to tackle the cosmological gravitino problem [1,2] caused by late gravitino decays in R-parity-conserving supersymmetry (SUSY) is to make the gravitino the lightest supersymmetric particle (LSP) and thus stable [3,4]. This leads to an attractive scenario where the gravitino accounts for the dark matter [5,6] whose density can match the observed one for a relatively high reheating temperature after inflation for gravitino masses in the GeV range [3,4]. However, the next-to-LSP (NLSP) tends to be long-lived due to the very weak coupling of the gravitino. In this case, late decays of the NLSP [7] and catalysis effects [8] can endanger the success of big bang nucleosynthesis (BBN); this is sometimes called the NLSP decay problem. The standard models of particle physics and cosmology successfully describe BBN. In the scenario with gravitino dark matter and a stau NLSP, τ 1 , the preservation of this success translates into bounds on the stau lifetime τ τ 1 and the stau yield Y after freeze-out from thermal equilibrium. For example, if τ τ 1 10 5 s, which corresponds to gravitino masses m G 300 GeV for a 1 TeV stau (or m G 20 GeV for m τ 1 = 300 GeV), the yield is required to be smaller than roughly 10 −15 [9,10]. While the stau is in thermal equilibrium in the hot early universe, its abundance decreases rapidly with time once the temperature falls below its mass. It freezes out from thermal equilibrium at a time determined by the cross section for the annihilation of staus and possibly other superparticles into Standard Model (SM) particles. As a consequence, smaller stau yields correspond to larger annihilation cross sections.
In this paper, we provide a classification of parameter space regions according to the stau yield, with emphasis on a thorough survey for regions where the yield is exceptionally low, i.e., much smaller than 10 −13 , which is the order of magnitude generically expected for a stau with mass around 100 GeV [11]. We work in the framework of the phenomenological Minimal Supersymmetric SM (pMSSM), whose parameters are defined at low energies and thus immediately applicable to the calculation of the annihilation cross section, without first calculating the running from some high scale. Systematically varying these parameters, we also allow all other superparticles to become nearly massdegenerate with the stau. Then co-annihilation effects can lead to large annihilation cross sections and correspondingly small stau yields.
We first consider cases without significant left-right sfermion mixing. Then the order of magnitude of the annihilation cross section is set by gauge couplings. The results can therefore equally be adopted for the case of a smuon or selectron NLSP. In a second part we allow significant left-right mixing for the third-generation squarks and finally for the staus themselves. In such cases, sfermion-Higgs couplings can become much larger than the gauge couplings, leading to a strong enhancement of the annihilation cross section. In our survey we will encounter the known exceptional regions with resonance effects [12] and enhanced Higgs-stau couplings [13,12], as well as regions with interesting co-annihilation effects that had not been studied in detail yet. We also investigate which combinations of effects produce interesting results.
Afterwards, we determine which parts of the parameter space are currently allowed. To this end, we perform a Monte Carlo scan and apply the relevant experimental and theoretical constraints. We require a CP -even Higgs with a mass around 125 GeV. We also consider the results of searches for the heavy Higgs particles of the MSSM. In order to take into account the LHC searches for heavy stable charged particles (HSCP) and R-hadrons, we compute the next-to-leading-order-corrected and next-to-leading-logresummed cross sections for the electroweak and strong production of superparticles and reinterpret the cross section upper limits for the 7 and 8 TeV runs reported by CMS [14] within our scenario. We consider the flavor and electroweak precision observables that are most relevant for the scenario, namely the W mass as well as the branching ratios of b → sγ and B 0 s → µ + µ − . As to theoretical constraints, we impose the absence of chargeand color-breaking minima in the scalar potential. Combining all results will finally allow us to determine the range of stau yields that is allowed by present constraints.
Our discussion reveals the most important parameters the stau yield depends upon and is thus a first step towards answering the question whether one could infer the stau yield from measurements at colliders. Together with a measurement of the stau lifetime, a determination of the yield could indicate whether the standard picture of the early universe is consistent or whether it must be extended. A number of such extensions have been proposed, for example, dilution of the NLSP density by late-time entropy production [15].
Although the gravitino-LSP scenario serves as our primary motivation to consider long-lived stau NLSPs, we would like to stress that this is not the only possibility. A long stau lifetime can occur in a number of scenarios with different super-weakly interacting LSPs. As our study is independent of the exact properties of the LSP, it can be applied to all these scenarios. However, the precise constraints from BBN do vary for different LSPs, which is why we will not consider them here.
The paper is organized as follows. Based on considerations about the freeze-out of MSSM sparticles, we will work out several general insights into the dependence of the stau yield on the model parameters in section 2. This discussion will set the stage for the deductive classification of the pMSSM parameter space in section 3. Afterwards, we will describe a scan in the 17-dimensional pMSSM parameter space with a stau NLSP, applying current theoretical and experimental constraints. Their effects on the yield will then be presented in section 5. We will conclude in section 6.

The freeze-out abundance of staus
In this section we will briefly review the physics of sparticle freeze-out and mention the specific assumptions for our setup. Suitably rewriting the solution of the Boltzmann equation we deduct guidelines for the systematic survey of the parameter space which we perform in section 3.
The abundance of staus around the time of BBN is considered to arise from the freeze-out [16,17,18,19,20,21,22,23] of the stau NLSP, the lightest sparticle of the MSSM. We assume that all MSSM particles have been in thermal equilibrium at some point during a hot phase in the early universe and that R-parity is exactly conserved. When the temperature of the universe decreases below the mass of the stau, the stau number density decreases exponentially. This exponential decrease is maintained as long as pair annihilation of staus are efficient enough to keep their number density close to the equilibrium number density. At the freeze-out temperature, which is typically of the order T f ∼ m τ 1 /25 [20], the stau decouples from the thermal bath and its number density freezes out. Consequently, the number density changes only because of the expansion of the universe, so the stau yield Y ≡ n τ 1 s remains (approximately) constant, where s is the entropy density. Afterwards, for the considered case of a metastable stau NLSP, it decays into the LSP once the Hubble parameter becomes comparable to the decay rate. This simple picture changes slightly when other MSSM sparticles are close in mass with the stau giving rise to co-annihilation effects [17,20]. In this case the annihilation of these sparticles competes with their decay into the stau, and a simultaneous freeze-out of several sparticles can occur. We assume that all heavier MSSM sparticles eventually decay into the stau NLSP and not directly into the LSP. Moreover, we will require that all other MSSM sparticles have a lifetime smaller than ∼ 10 −2 s (correspondingly, Γ 10 −22 GeV) ensuring that none of these decays take place during or after BBN. Consequently, all changes on BBN solely depend on the stau yield. With these assumptions the desired number density of staus is simply the sum of the number densities of all relic sparticles that survive freeze-out. Although our requirement on the lifetime of the other sparticles will set a lower limit on the mass degeneracy of co-annihilating sparticles with the stau, in this section we will blithely consider exact mass degeneracy as a limiting case, keeping in mind that a certain separation is in fact required. The effects of taking all current limits into account are discussed in section 5.
An approximate solution of the Boltzmann equations can be written as [21,23] 1 x 0 x f dx πM Pl 8ḡ 45 where x = m τ 1 /T and x f , x 0 are the corresponding quantities at the point of freeze-out and at the desired point of observation, respectively. Furthermore, σ eff v Møl is the thermally averaged cross section times Møller velocity andḡ is a degrees of freedom parameter as defined in appendix A. Besides, M Pl = (8πG N ) −1/2 is the reduced Planck mass. The stau yield at freeze-out, Y (x f ), can only be computed by solving the Boltzmann equation numerically. For this, we will make use of micrOMEGAs 2.4.5 [24]. Note that this program automatically sets x 0 = ∞, which is a good approximation as long as the stau is not extremely short-lived, τ τ 1 10 −4 s. An approximate solution can be found by neglecting the term 1/Y (x f ) [21]. For the following discussion it is instructive to rewrite eq. (2), where we closely follow [25]. We can express the yield as Y ∝ m τ 1 x 0 where we have introduced the dimensionless thermally averaged cross section [23] Here, i, j run over all supersymmetric particles (including the stau and the antistau) involved in the (co-)annihilation with mass m i = x i T and internal degrees of freedom g i . Besides, K n are the modified Bessel functions of the second kind of order n and z = √ s/T .
Note that σ eff v x = x −2 m 2 τ 1 σ eff v Møl . The rescaled cross sectionσ is connected to the (usual) annihilation cross section bỹ it is a function of dimensionless quantities only, where {a SUSY } denotes a set of SUSY parameters each normalized by m τ 1 , and {a SM } is a set of SM parameters, normalized in the same way.
Considering eqs. (3) and (4) there are three interesting observations that can be made and that are important for the subsequent discussion.
1. For fixed a SUSY , and ifσ ij is asymptotically independent of a SM for a SM → 0 (which corresponds to the limit of large SUSY masses) 1 ,σ ij only depends on the ratio x/z. From this, it follows that the yield is simply proportional to the stau mass, up to effects induced by the dependence of the choice of x f and the correction 1/Y (x f ) on the stau mass, as well as effects of order a SM .
2. For a fixed initial state i, j, opening up additional final state channels can only enhance the cross sectionσ ij and thus lower the yield. Here, m A plays an important role, as it determines the masses of the heavy Higgses, that could be either below or above their pair-production threshold. All other R-parity even particles are considerably lighter than the stau, if we take current direct search limits into account which require m τ 1 340 GeV [14]. 3. In contrast to opening up new final state channels, introducing additional initial state channels in the presence of co-annihilation effects can either raise or lower the yield depending on the involved cross sections and the additionally introduced degrees of freedom. For simplicity, we consider the limiting case of an exact mass degeneracy where the co-annihilation effects are maximal. In this case where we have introduced Additional initial states can only lower the yield if they introduce large cross sections σ ij v x in the numerator that are capable of overcompensating the introduction of additional terms in the sum over the degrees of freedom in the denominator. For instance, the mere introduction of more sparticles with similar interactions as the stau cannot decrease the yield further 2 -in contrast, if there are combinations i, j that lead to smaller cross sections σ ij v x than σ τ 1 τ * 1 v x , the numerator in eq. (9) increases less than the denominator and we obtain a net increase of the yield. This is the case in the slepton co-annihilation region, as we will discuss below. Considering co-annihilating sparticles i, j that introduce cross sections much larger than the stau-stau annihilation cross section, σ ij v x ≫ σ τ 1 τ * 1 v x , we can approximate eq. (9) by The introduction of more and more sparticles i, j of the same kind could only lead to an asymptotical increase of σ eff v x towards the value we would obtain by neglecting the stau degrees of freedom in the denominator altogether. However, this saturation would only be achieved if all cross sections σ ij v x were equally large. This is usually not the case. For instance, introducing an additional squark generation in a squark co-annihilation scenario effectively reduces σ eff v x due to the smaller inter-generational interactions, i.e., the smaller cross sections σ ij v x for i, j belonging to different generations. This is an important fact which severely restricts the possibilities for exceptionally small stau yields as the result of co-annihilation effects and enables an economical discussion of the various combinations of coannihilating sparticles that could occur in the general pMSSM. In particular, it allows us to study the different cases in an isolated way.
After we have discussed the general scaling behavior of the yield with the initial state sparticle masses as well as the behavior under the introduction of additional initial and final states, in the following we will comment on the parameters that govern the size of the cross sectionσ ij for fixed initial and final states.
There are basically two ways how the free SUSY parameters a SUSY could affect the cross sectionσ ij . One is the strength of the involved couplings. Besides the known SM gauge couplings the MSSM contains the couplings of the sfermions f to the Higgses which involve the trilinear soft terms A f , tan β and the higgsino mass parameter µ as a priori free parameters of the theory. Although subject to constraints (see section 4) these couplings can be very large [13,12] and the resulting cross sections can even be larger than the ones for processes dominated by the strong interaction. All other couplings in the theory are given by SM gauge couplings (multiplied by possible suppression factors ≤ 1 due to mixings) or are proportional to the mass of the involved particle and thus do not introduce additional free parameters. The second is the appearance of non-SM particles in the intermediate states of the annihilation processes. On the one hand SUSY particles can appear in the t-channel of the annihilation processes. On the other hand the heavy Higgses whose masses are determined by m A can appear in the s-channel. Especially the latter effect can lead to a drastic enhancement of the cross section close to the resonant pole m A ≃ 2m τ 1 [12] (or m A ≃ 2m co-ann for the case of a co-annihilating particle).
With these general remarks in mind, we will now systematically explore the different distinct regions in the pMSSM parameter space.

A systematic survey in the pMSSM
In this section we will give an overview of different regions in the pMSSM parameter space characterized by the physical processes governing the stau yield and the resulting ranges of values.
• In section 3.1 we will consider the case of no sfermion mixings and no co-annihilation effects. In the case of no sfermion mixings, all couplings are determined by the known gauge couplings and masses of the involved particles, as well as the field decomposition in the EW gaugino sector. Then the only parameters governing the yield are the stau mass and the masses and mixings of EWinos appearing in the t-channel diagrams. We consider a purely right-handed stau as well as a purely left-handed stau and investigate the dependence of the yield on the variation of both stau masses and the EWino mass parameters M 1 , M 2 and µ.
• In section 3.2 we allow for co-annihilation effects by approaching the region a i 1.1, where we successively consider co-annihilation with sleptons, gauginos, and squarks. Naturally, the yield in this case additionally depends on the respective value(s) of a i and thereby the mass of the additional co-annihilating sparticle(s), as well as the handedness of the NLSP. Furthermore, for the case of EWino co-annihilation we vary m A and thus examine the effect of the opening up of additional Higgs final state channels and Higgs resonances.
• In section 3.3 we finally allow for significant left-right mixing of the sfermions, thereby enabling large sfermion-Higgs interactions either in the third-generation squark sector (in a co-annihilating setup) or in the stau sector. In both cases, by varying m A , effects of additional Higgs final state channels and Higgs resonances are studied. Besides the respective mixing angles and masses, in addition there is a strong dependence for diagrams involving light or heavy CP -even Higgses in the final or intermediate state on the A-parameters of the third-generation squarks as well as µ and tan β.
If not stated otherwise, we set all rescaled mass parameters a i ≡ m i /m τ 1 that are not under consideration to a i = 4, which we consider sufficiently large to ensure that processes containing these particles are significantly suppressed and generically do not contribute to the stau yield. All spectra in this section are calculated using SuSpect 2.41 [26] at leading order. 3

Stau pair annihilation in the absence of left-right mixing
In this subsection we consider stau annihilation in the absence of large left-right mixing and co-annihilation effects with other sparticles. 4 If not stated otherwise, in order to achieve purely right-or left-handed mass eigenstates we choose a low value for tan β (tan β = 2) and enforce the cancelation A τ = µ tan β, so that X τ ≡ A τ − µ tan β is zero and thus the stau mass matrix is diagonal, cf. appendix B.
In this case, the stau yield depends only on the stau mass and on the masses and mixings of EWinos appearing in the t-channel of the annihilation processes. Figure 1 shows the stau yield as a function of the stau mass for a purely right-handed and purely left-handed lighter stau τ 1 . For this plot, the tau sneutrino mass is set to m τ 1 by hand for τ 1 = τ L . 5 As we expect from the discussion above, the stau yield has an almost linear dependence on the stau mass. In fact, the expressions 3 For the study of idealized cases in this section we switched off the higher-order corrections in the spectrum generation by setting ICHOICE(7)=0. The computation of the yield in the Monte Carlo scan in section 5 contains the full radiative corrections provided by SuSpect. 4 For a left-handed lighter stau we have to take sneutrino co-annihilation effects into account. 5 Strictly speaking, this choice should be considered as a limiting case which is not a valid point in the MSSM, since m ντ < m τ 1 for a purely left-handed lighter stau. However, it approximates nearby valid points with an almost left-handed lighter stau and a slightly heavier sneutrino.  Figure 1: Left panel: Stau yield Y as a function of the stau mass m τ 1 for a right-handed lighter stau as well as for a left-handed lighter stau mass-degenerate with the tau sneutrino (to be considered as a limiting case for realistic spectra). All other SUSY mass parameters are set to 4m τ 1 . Right panel: Effects of EWinos in the t-channel. The curves show the yield as a function of Mi/m τ 1 for the bino mass parameter, i = 1, with τ1 = τR (blue solid curve) and τ1 = τL (black dot-dashed curve) as well as for the wino mass parameter, i = 2, with τ1 = τL (green dotted curve). All curves are normalized to their respective value at Mi/m τ 1 = 4.
for a right-handed lighter stau and for a left-handed lighter stau describe the results in the given range at a percent level fit accuracy.
In eq. (8), we have argued that we generically expect a constant scaling of the yield with respect to the stau mass if all other SUSY parameters are fixed and the effective cross section is independent of SM-like scales, which are parametrized by a SM . Therefore, the deviation from this expected scaling behavior in the above expressions calls for a more thorough investigation. Indeed, if the cross section defined by eq. (5) exists and is finite in the limit a SM → 0, eq. (3) is expected to hold independent of the x-dependence of σ eff v x . Consequently, these deviations must come from the approximations employed in deriving eq. (3), i.e., from the freeze-out approximation with constant x f and from neglecting 1/Y (x f ). Indeed, the value of x f chosen by micrOMEGAs varies with m τ 1 .
Effects of varying the mass of the bino and wino appearing in the t-channel diagrams are shown in the right panel of figure 1. The additional t-channel diagrams have the effect of lowering the stau yield. For the right-handed stau, when lowering M 1 the channel τ 1 τ 1 → τ τ becomes more important reaching a maximum of around 65% for an approximate degeneracy of the bino-like neutralino and the stau, i.e., M 1 /m τ 1 ≃ 1.1 (where co-annihilation effects are still small). Since the other channels (listed above) are not affected by the variation of M 1 , their absolute contribution remains unchanged and the yield drops by a factor of roughly 0.5, accordingly. For M 1 /m τ 1 = 1.1, the pre-factor in eq. (12) would become 8.55 × 10 −13 , in rough agreement with the results given in [27,11]. For a completely decoupled bino the yield pre-factor in eq. (12) would change to 2.24 × 10 −12 according to the missing channel τ 1 τ 1 → τ τ . For the left-handed stau the t-channel diagrams are less important in comparison. For the case of small M 2 (but again above the region where co-annihilation is efficient) the t-channel processes τ 1 ν τ → τ ν τ , ν τ ν τ → ν τ ν τ and τ 1 τ 1 → τ τ (each of which contributes around 13%) become the most important processes followed by τ 1 τ 1 → W W (10%). In contrast, lowering the bino mass does not have a large effect on the yield of the left-handed stau; in particular, the respective t-channel processes do not become the leading contributions. The (absolute) yield for right-handed staus even becomes smaller than the yield for left-handed staus for M 1 /m τ 1 ≃ 1.1 [12]. In all cases the exponent of the m τ 1 -dependence of the yield stays approximately constant when varying M 1 /m τ 1 or M 2 /m τ 1 .

Co-annihilation regions
Co-annihilation effects can be important whenever the mass splitting ∆m between the stau NLSP and the next-heavier sparticle(s) is of the order of the freeze-out temperature, ∆m/m τ 1 ≃ x −1 f . Given that the typical freeze-out temperature corresponds to x f ≃ 25, co-annihilation effects are expected to be significant for relative mass degeneracies of around 5-10% [20].
We will now systematically investigate how the stau yield changes if additional coannihilating sparticles are introduced. Further, exemplarily we show how simple estimates including the exact consideration of all degrees of freedom can successfully predict the relative change in the yields. In the following we consider Y τ 1 = τ R and Y τ 1 ≃ τ L from eqs. (12) and (13) as reference yields that we normalize our results to. We choose m τ 1 = 1000 GeV and we systematically vary the mass ratios m i /m τ 1 for the different sparticle species i in order to study co-annihilation effects of the sparticles in the MSSM in an isolated way. If not stated otherwise we vary the corresponding soft masses and plot the physical sparticle mass. If several sparticle masses are governed by one parameter that is subject to variation, we plot the smallest among these sparticle masses, if not stated otherwise. For example, if we vary the soft mass of the left-handed sleptons, we plot the sneutrino mass. 1

Co-annihilation with sleptons of the first and second generation
The upper panels of figure 2 show the co-annihilation of staus with right-and lefthanded sleptons of the first and second generation. The stau yield increases with an increasing importance of co-annihilation effects. This can be understood as follows. As an example, let us consider the case of a right-handed lighter stau which is mass-degenerate with the right-handed selectron and smuon. In the limit of complete degeneracy the denominator in eq. (9) is enhanced by a factor of 9. On the other hand, not all cross sections σ ij v x (i, j = 1, 2, 3 denote the generations) are equally large, as for i = j only t-channel diagrams contribute. Accordingly, in the limit of a completely decoupled bino, the numerator in eq. (9) increases only by a factor of 3 and hence, since Y ∝ σ eff v −1 x , the yield would increase by a factor of around 3 with respect to the non-degenerate case. However, for the considered case of M 1 ≃ 4m τ 1 , the t-channel neutralino contribution is important. The six channels ℓ i ℓ j → ℓ i ℓ j contribute around 7.7% each. For the choice m B ≃ 1.1m τ 1 , 7 each channel contributes around 13%. Thus, the net increase of the stau yield is milder in the presence of the t-channel bino contributions and turns out to be 2.3 and 1.8 for the former and latter choice for m B , respectively. This agrees with the findings in [11,28].

Co-annihilation with gauginos
In this paragraph, we discuss the effects of co-annihilation of staus and gauginos from the electroweak as well as strong gauge groups. While the yield is generically lowered for sizeable co-annihilation cross sections, a special case is given for the co-annihilation with bino eigenstates, where we obtain an increase in the overall yield.
The lower panels of figure 2 show the co-annihilation effects for gauginos. We vary M 1 , M 2 , µ and M 3 and plot the mass of the lightest EWino (which is the lightest neutralino) and the gluino, respectively. Hence, we consider (almost) pure gauge eigenstates in the EWino sector. The effects of large mixing in the EWino sector will be discussed in section 3.2.4. Due to the smaller annihilation cross section of the bino the net effect of bino co-annihilation increases the yield. This is, again, due to the increase in the degrees of freedom by a factor of 2 (for the right-handed lighter stau) or 3/2 (for the left-handed lighter stau, accompanied by the tau sneutrino). At the same time the pair-annihilation of the binos is negligible and the associated co-annihilation of neutralinos with staus (or with tau sneutrinos) is sub-leading-in the limit of complete degeneracy the corresponding contributions add up to less than 25% in the case of a right-handed light stau and 20% in the case of a left-handed light stau. Consequently, for the right-handed stau we expect from eqs. (3) and (9) (In this and the following estimates we make use of the fact that for annihilation processes 7 m B denotes the mass of a bino-like lightest neutralino. without resonant or threshold effects the thermally averaged cross section can be expanded in 1/x where the leading contribution is independent of x [21].) At m B ≃ 1.1m τ 1 , where co-annihilation is already inefficient but the t-channel neutralino contributions are (still) maximal, Using g τ R = 1 and g B = 2, we obtain in good agreement with figure 2 (lower right panel, dot-dashed curve). For τ 1 = τ L , a similar estimate yields a net increase in the yield by a factor of around 2 with respect to the yield at m B ≃ 1.1m τ 1 . In contrast, for wino co-annihilation the annihilation of the wino-like neutralino and chargino among themselves is the dominant contribution. For a right-handed lighter stau these contribute almost 100% to the annihilation cross section while for a left-handed lighter stau they contribute more than 50% followed by associated co-annihilation processes of neutralino and chargino with the lighter stau and the tau sneutrino amounting to a contribution around 30%. Hence, due to the larger annihilation cross sections of wino-like EWinos the stau yield is significantly reduced despite the 6 additional degrees of freedom introduced by the mass-degeneracy of one neutralino and chargino.
For higgsino co-annihilation the relative importance of annihilation and co-annihilation processes is not vastly different. However, due to the 8 additional degrees of freedom the net reduction of the stau yield is less. In the case of a left-handed lighter stau the yield even increases slightly at around m H /m τ 1 ≃ 1.05.
For the case of gluino co-annihilation the situation is even more pronounced than for the wino case, since σ g g v x ≫ σ τ 1 τ 1 v x and σ g τ = 0. In the case of a mass-degenerate gluino, eqs. (3) and (9) yield Y ∝ (g g +2g τ ) 2 /g 2 g . Accordingly, the yields for a left-handed and a right-handed lighter stau differ only due to the extra degrees of freedom of the tau sneutrino, Y τ 1 = τ L /Y τ 1 = τ R ≃ (g g + 2g τ 1 + 2g ντ ) 2 /(g g + 2g τ 1 ) 2 ≃ 1.2. (The relative yields in figure 2, which are normalized to the respective yield without co-annihilation, show a larger difference, due to the difference in the reference yields.) Gluino pair annihilation processes are dominant up to a relative mass difference to the stau of 7% and 6% for a right-handed and left-handed lighter stau, respectively.

Co-annihilation with squarks
In the upper left panel of figure 3 we show the co-annihilation effects of the first two generation squarks for a right-handed stau. We vary the soft masses m Q 1,2 , m u 1,2 and m d 1,2 and plot the mass of the lightest among the squarks whose mass is dictated by the respective parameter. Although the involved strong interactions lead to relatively large cross sections (and in particular σ(qq → X) ≫ σ( τ 1 τ 1 → X ′ )) the decrease in the   , wino (green dotted line) and gluino (blue solid line) mass parameter in a co-annihilation scenario with third-generation squarks. We adjusted m u 3 and m Q 3 such that the corresponding lighter sparticle (the stop and sbottom, respectively) is exactly mass-degenerate with the stau. Lower right panel: Effects of the presence of squarks in the t-channel of gluino co-annihilation diagrams. We adjust the gluino to be exactly mass-degenerate with the stau and varied m u 1,2 (green dotted curve), m Q 1,2 (blue solid curve) and all soft parameters of the three squark generations, namely m u 1,2,3 , m d 1,2,3 and m Q 1,2,3 simultaneously (red dashed curve). yield is significantly less pronounced than in the case of gluino co-annihilation. We can understand this as follows, considering the case of a full degeneracy of the stau with the right-handed up-type squarks of the first two generations. The dominant annihilation channels in this case are u R u R , c R c R → gg and contribute 72% to the annihilation cross section. We compare this contribution with the case when there is no co-annihilation with squarks. In the latter case τ 1 τ 1 → γγ contributes 38% to the total annihilation cross section. From these numbers we can estimate the expected reduction of the yield. From eqs. (3) and (9), for the case of co-annihilation and for the case without co-annihilation, and thus We approximated the ratio between the two cross sections by the unaveraged cross sections in the non-relativistic regime as computed by CalcHep [29]. This estimate comes very close to the value that is displayed in the upper left plot of figure 3 (purple dashed line).
Since for a close mass degeneracy the pair annihilation processes of squarks dominate over stau pair annihilation and associated squark-stau annihilation, the absolute stau yield for a left-and right-handed stau are virtually identical. This is why we refrain from showing the corresponding plot for the left-handed stau in figure 3. The main difference in such a plot would arise from the mere difference in the reference yields for left-and right-handed staus.
The difference between the reductions of the yield for up-and down-type right-handed squarks arises solely from the different cross sections, since the number of degrees of freedom is exactly the same. The difference is induced from subdominant channels containing γ g and Z g final states. These contributions are sensitive to the charge of the corresponding squarks, leading to a smaller yield for the up-type squarks.
In the case of degenerate left-handed squarks additional annihilation channels open up, namely the annihilation of up-type-down-type squark pairs arising from diagrams with t-channel squarks or charginos as well as the four vertex contact interactions u L d L → W g. We found that the cross section for these processes containing electroweak interactions are almost as large as those induced by strong interactions. Furthermore, we observed a constructive interference between gluino and wino exchanging t-channel diagrams. In fact, this leads to a significant increase of the effective thermally averaged cross section with respect to the case of right handed squarks which overcompensate the doubling in the degrees of freedom. Finally, for the case m Q 1,2 = m u 1,2 = m d 1,2 we observe a clear increase in the yield as a mere result of an increase in the degrees of freedom relative to the cases considered before.
The upper right panel of figure 3 shows the co-annihilation effects of the thirdgeneration squarks. The relative behavior of the yield for the case of b R , t R and ( b, t) L is comparable to the yield in the corresponding cases for degenerate first two generations. The overall decrease in the yield due to the co-annihilation effects is smaller as the inter-generation initial states, for which the cross section is considerably smaller, are absent. It is interesting to note that channels with Higgs particles in the final states, which contain contact term interactions arising from the F -and D-terms in the scalar potential, are not suppressed in the absence of sfermion mixing. However, they do not play an important role even for the stop although the diagram arising from the F -term is proportional to the Yukawa coupling squared.
Effects of gauginos appearing in the t-channel of the squark annihilation processes are shown in the lower left panel of figure 3. For the blue solid and green dotted curve we fixed m Q 3 such that m b 1 = m τ 1 (m t 1 is slightly larger) and varied M 3 /m τ 1 and M 2 /m τ 1 , respectively. For the black dot-dashed curve we fixed m u 3 such that m t 1 = m τ 1 and varied M 1 /m τ 1 . 8 All curves are normalized to the respective values at a i = 4. As in the case of slepton annihilation, the t-channel contributions increase the effective annihilation cross section for small gaugino masses. However, the relative effect is smaller, cf. right panel of figure 1. For M i /m τ 1 1.1 co-annihilation effects of gauginos become important. For bino and winos these effects increase the yield due their smaller annihilation cross sections relative to those of the squarks on the one hand and the additional degrees of freedom on the other hand. For the gluino, co-annihilation effects lead to a further reduction of the yield despite the additional degrees of freedom. In order to further understand the interplay between squarks and gluinos, in the lower right panel of figure 3 we fixed m g = m τ 1 and varied certain squark masses (according to our convention, all others are kept at 4m τ 1 ). The contributions from t-channel squarks in the gluino-annihilation processes cause a destructive interference. For the green dotted and blue solid curve we vary the soft parameters m u 1,2 (m d 1,2 obviously give the same result) and m Q 1,2 , respectively. For the red dashed curve we varied all (bilinear) soft mass parameters of the first-to third-generation squarks simultaneously. Interestingly, among scenarios with squark and gluino co-annihilation a scenario with a mass-degenerate gluino and decoupled squarks would have the smallest stau yield.

Varying m A in the case of EWino co-annihilation
We now vary the parameter m A in order to investigate the potential for changes in the cross sections of the EWino co-annihilation scenario due to additional intermediate and final states, especially around the resonant pole of an s-channel heavy Higgs and below the threshold for heavy Higgs final states. The EWino couplings to the Higgs are always of the type H W Φ or H BΦ, where W and B denote the wino and bino gauge eigenstates, H the higgsino gauge eigenstate and Φ a Higgs. Hence, pair-annihilation of EWinos into Higgses requires either a substantial higgsino admixture for the lightest EWinos or the presence of a higgsino-like EWino sufficiently light to significantly contribute in the t-channel. However, as we do not observe any change in the yield and the relative contributions when passing the heavy Higgs production threshold EWino pair annihilation into heavy Higgs final states is not an important channel (see figure 4 at m A /m τ 1 1).
Resonance effects occurring for m χ ≃ 2m A are only important in the case of a large higgsino admixture in the lightest EWinos participating in the pair co-annihilation processes. Figure 4 shows the yield in a co-annihilation scenario where the lightest EWino  is a bino-higgsino mixture (M 1 = µ, M 2 = 4m τ 1 ) or a wino-higgsino mixture (M 2 = µ, M 1 = 4m τ 1 ). We vary m A and show the yield as a function of m A /m τ 1 , for a complete degeneracy m χ 0 1 = m τ 1 and for a relative deviation of 5%, m χ 0 1 = 1.05m τ 1 . The resonant EWino co-annihilation can lower the yield by more than two orders of magnitude. This is analogous to what happens in the H/A-funnel region of a neutralino LSP scenario [30,31,32,33,34,35]. These results were obtained for tan β = 2. However, we found very similar results with tan β = 40, although with a slightly shallower dip in the resonance. For M 2 = µ as well as tan β = 2 and tan β = 40 the dominant annihilation channel is χ ± 1 → tt and χ ± 1 → bb in the resonance, respectively. The small dip in the curves for M 1 = µ slightly below the main resonance is caused by resonant annihilation of staus via the CP -even heavy Higgs. 9 For the case M 1 ≃ M 2 ≃ µ, not shown here, the yield tends to be larger again due to the additional degrees of freedom.
The right panel of figure 4 shows the annihilation contributions for M 2 = µ and m χ 0 1 = m τ 1 . The contribution χ χ → HX denotes all channels with EWinos in the initial state and exactly one of the Higgs fields H 0 , H ± , A 0 in the final state. (Channels with two Higgs fields in the final state contribute negligibly.) Independent of m A , the contribution of channels with one light Higgs h in the final states is roughly a fifth of the remaining contributions denoted by 'others' in the plot.

Co-annihilation with mixed stops
Still restricting ourselves to the case of small left-right mixing of the staus, we will now discuss the case of co-annihilation with squarks that acquire substantial left-right mixing. The potentially large couplings of sfermions to the Higgses are proportional to the leftright mixing and proportional to the parameters appearing in the off-diagonal terms in the mass matrix. We assume no particularly large X f in the first two generations and restrict the discussion to the third-generation sfermions, i.e., to the case of a coannihilating sbottom or stop. The couplings of the sbottom and stop to the neutral, CPeven Higgses h, H are summarized in appendix C. In the decoupling limit, M Z ≪ m A , and for enhanced Higgs-sfermion couplings, these couplings can be approximated by We will here exemplarily focus on the stop. We do not encounter potentially larger enhancements of the couplings for sbottom co-annihilation. The smaller Yukawa coupling, in contrast, tends to require larger SUSY parameters in order to obtain the same coupling strength. Moreover, the couplings of the sbottom and a stau are similar in the sense that they can become very large for large tan β. In this concern it is more interesting to study the case of the stop being important in a complementary corner of the parameter space, namely for smaller tan β. Furthermore, a significant left-right mixing of the stops is preferred from the requirement of large radiative corrections to the Higgs mass when interpreting the Higgs discovered at the LHC as the lighter neutral, CP -even Higgs h. Figure 5 shows the stau yield for a co-annihilating stop which is completely massdegenerate with the stau NLSP. At tree-level, the leading contribution of the coupling of the stop to the light Higgs, eq. (22), can be expressed solely by the spectrum parameters by using the analogon of eq. (68) for the stop sector. We chose m t 1 = m τ 1 and varied θ t for different choices of the mass of the second stop, m t 2 . We set θ t by fixing tan β = 5 (as well as µ = 4m τ 1 as usual) and setting A t accordingly. Note that this treatment of the parameters implicitly determines the mass of the lighter sbottom, so further coannihilation effects can take place that potentially increase the yield. The Higgs mass was set to m h = 126 GeV 'by hand'. The result for the yield is, however, not sensitive to the actual value of the Higgs mass. (The implications of the requirement to actually obtain this Higgs mass from the radiative corrections in the stop sector are discussed in section 5.) As discussed before, strong interaction can already lead to a reduction of the yield by an order of magnitude. In the case of large stop-Higgs couplings the corresponding  processes become dominant and lead to a further significant reduction. We obtain stau yields of 5 × 10 −14 and less than 10 −16 for a mass splitting of m t 2 /m t 1 = 1.5 and 3, respectively. The lower panels of figure 5 show the relative contributions to the annihilation. Due to the exact mass-degeneracy of the stop with the stau, the pair-annihilation processes of stops dominate over annihilation process involving the stau. For a relatively small mass gap between the lighter and the heavier stop, m t 2 /m t 1 = 1.5, i.e., for a moderate coupling C[h, t 1 , t 1 ] m t 1 , annihilation into gluino pairs and pairs of vector bosons are the dominant channels. For very large mass gaps, m t 2 /m t 1 = 3, i.e., for large values of C[h, t 1 , t 1 ], the channel t 1 t 1 → hh becomes important. In this regime the leading contribution from t 1 t 1 → hh comes from the pair annihilation of stops via the t-channel diagram. This contribution involves two stop-stop-Higgs vertices. The cross section is therefore proportional to C[h, t 1 , t 1 ] 4 while all contributions with an s-channel h are only mH/m τ 1 rel. contribution proportional to C[h, t 1 , t 1 ] 2 . For this reason the channels t 1 t 1 → V V and t 1 t 1 → tt, bb become less important with larger mass gaps and larger left-right mixing of the stops.
The upper right panel of figure 5 shows the mixing parameter |X t |/M s , where M s = √ m t 1 m t 2 , corresponding to the lines drawn in the left panel. For m t 2 /m t 1 = 1.5 and 3 the mixing parameter is |X t |/M s = 3 and almost 14, respectively. |X t |/M s ≃ √ 6 maximizes the positive radiative corrections to the Higgs mass [36,37] and thus is preferred in the absence of overly large stop masses.

Varying m A in the case of co-annihilation with mixed stops
If we relax our assumption m A ≃ 4m τ 1 we can study the effects of heavy Higgs resonances and of opening up channels with heavy Higgs final states. Figure 6 shows the yield in a co-annihilation scenario where the stop is maximally mixed, θ t = π/4 or 3π/4, and m τ 1 = 1 TeV as a function of m H /m τ 1 . We show the relative yield for an exact degeneracy (blue, solid and red dashed curves) as well as for m t 1 = 1.05m τ 1 (green, dotted and black, dot-dashed curves). We choose two sets of parameters, one with tan β = 2 and A t = 4 TeV (blue, solid and green, dashed curves) and one with tan β = 20 and A t = −4 TeV (red, dashed and black, dot-dashed curves). For both cases we set A τ = A b = 0 and µ = 4 TeV. The soft parameters m Q 3 and m u 3 are determined by tree-level relations from the desired m t 1 and θ t . Again, we use an iterative algorithm in order to control these parameters after spectrum generation.
In the right panel of figure 6 we show the relative contributions of the annihilation channels for the case m t 1 = m τ 1 , tan β = 2 and A t = 4 TeV. The red curve is the contribution of all channels that are not explicitly displayed. Its pronounced peak at around 2250 GeV is caused by the channel t 1 b 1 → tb, which contributes around 38%.
The mass of the sbottom is around 1150 GeV. Similar to the case of EWino co-annihilation, we obtain a strong reduction of the yield in the presence of a resonant pole. In contrast, we also see a decrease of the yield below the threshold for heavy Higgs final states. This effect is only pronounced for small tan β since the annihilation into final state heavy Higgses contributes significantly only for very large stop-Higgs couplings.

Large stau mixing
We will now discuss large mixing in the stau sector itself. Accordingly we will switch off any avoidable effect of co-annihilation. In the decoupling limit, and for enhanced Higgs-stau couplings, these couplings read approximately We first vary the stau mixing angle while keeping m A ≃ 4m τ 1 . Analogous to the case of the stop, we perform the scan for different choices of m τ 2 /m τ 1 . We choose tan β = 20, µ = 4m τ 1 and achieve the required X τ by choosing A τ accordingly. Figure 7 shows that the yield can be reduced by several orders of magnitude for large mass splittings and significant left-right mixing, i.e., large couplings C[h, τ 1 , τ 1 ]. This result was first discussed in [13,12]. In [12] the authors equally scanned over θ τ but chose a fixed value for X τ . The results obtained there are compatible with ours. The upper right panel of figure 7 displays the size of |X τ | which is required in order to provide the fixed ratio m τ 2 /m τ 1 when varying θ τ . This reveals that low stau yields can only be obtained for very large values of |X τ |.
The upper left panel of figure 7 shows that for a moderate mass ratio m τ 2 /m τ 1 , the curves are not symmetric around θ τ = π/2. This effect arises from the interference term of a heavy Higgs in the s-channel of the annihilation processes τ 1 τ 1 → tt, bb, hh. While the coupling C[h, τ 1 , τ 1 ] is completely symmetric around θ τ = π/2, C[H, τ 1 , τ 1 ] is not. The asymmetry could, however, be reduced or removed by another choice of A τ , µ and tan β to achieve the required X τ , or by a stronger decoupling of m A .

Varying m A in the case of large stau mixing
If we relax our assumption m A ≃ 4m τ 1 the contributions with heavy Higgs intermediate or final states can dominate the annihilation cross section. On the one hand, the heavy Higgs can appear in the s-channel leading to a resonant pole in the propagator when m H ≃ 2m τ 1 [12]. On the other hand, heavy Higgses can appear in the final state around or below threshold, i.e., when m h + m H 2m τ 1 or m H m τ 1 . The upper left panel of figure 8 shows the yield for a maximally mixed stau, θ τ = π/4, and m τ 1 = 1 TeV as a function of m H /m τ 1 . For small values of tan β the yield does not significantly deviate from the one for a right-handed stau, except for the resonance where the yield is reduced by up to more than two orders of magnitude. The upper right panel in figure 8 shows the relative contributions to the annihilation. For most of the displayed range of m H /m τ 1 co-annihilation channels contribute the most (red dot-dot-dashed curve). This is caused  Aτ /m τ 1 = 8 tan β = 5 mH/m τ 1 rel. contribution  by the relatively small X τ that requires a small mass splitting of m τ 1 , m τ 2 and m ντ in the presence of maximal mixing. In the resonance, the channel τ 1 τ 1 → tt dominates (black dot-dashed curve). Note that the peak in the contribution of co-annihilation channels slightly above the resonance m H /m τ 1 = 2 stems from the resonant annihilation of τ 2 and ν τ , which are slightly heavier than the lighter stau.
For tan β = 50 we obtain a reduction of the yield by about four orders of magnitude. This result is independent of the chosen sign of A τ and therefore independent of the sign of the coupling C[H, τ 1 , τ 1 ]. Since the couplings of the heavy Higgs to the bottom quark are proportional to tan β in the decoupling limit, the dominant channel for tan β 8 is τ 1 τ 1 → bb. Since the coupling C[H, τ 1 , τ 1 ] is also proportional to tan β, for very small stau yields in the resonance region we typically obtain bb final states.
Another interesting observation can be made in the region m H /m τ 1 1. Below the threshold for two heavy Higgses in the final state the stau yield is significantly reduced in the case of a negative A τ (blue solid curve) while this feature is not present for positive A τ (red dashed curve). (The other parameters are identical.) This asymmetry is due to an interference of the t-channel diagram τ 1 τ 1 → HH with the s-channel diagrams τ 1 τ 1 → h, H → HH. The diagram τ 1 τ 1 → H → HH is sensitive to the sign of C[H, τ 1 , τ 1 ] and introduces a constructive (destructive) interference for C[H, τ 1 , τ 1 ] negative (positive). When decreasing the mass of the heavy Higgs, this diagram is reduced by the increasing denominator of the heavy Higgs propagator.

Differences in the scaling behavior
In all the processes discussed in sections 3.2 and 3.3 we have set the stau mass to m τ 1 = 1 TeV. Point 1 in the list of observations in section 2 implies that for fixed ratios of SUSY masses and for an annihilation cross section that is independent of the masses of SM particles in the limit a SM → 0, the results can be extrapolated to any value of m τ 1 by a simple rescaling of the yield that is approximately linear in the stau mass. Indeed, we explicitly checked the scaling behavior of all limiting cases considered in section 3.2 and found where δ ≃ 0.9 as in eqs. (12) and (13). However, the argument does not apply for non-vanishing left-right mixing in the sfermion sector, since the limit a SM → 0 would imply sin 2θ f → 0, cf. eq. (68) in the limit m τ /m τ 1 → 0. The mixing term introduces an explicit scale dependence, since it is proportional to the fermion mass. However, eq. (26) holds approximately if we keep the ratios of parameters fixed that govern the annihilation processes sensitive to the left-right mixing. As an example, for large m H , i.e, no resonance effects, the leading term for stau pair annihilation in the presence of large stau-Higgs couplings is annihilation into light Higgs final states via t-channel stau exchange, which behaves like [13]σ where we have additionally adopted the limit m h ≪ m τ 1 . In figure 9 we plotted Y /m 0.9 τ 1 as a function of m τ 1 between 300 GeV and 10 TeV for the case of maximal left-right mixing, θ τ = π/4, andĈ h, τ 1 /m τ 1 fixed to four different values. Although the scaling is slightly changed forĈ h, τ 1 /m τ 1 ≥ 1, eq. (26) remains a reasonable approximation. The slight deviations arise from the interplay of different contributions, most importantly stau pair annihilation into h, W and Z. ForĈ h, τ 1 /m τ 1 = 0.3, diagrams involving stau-Higgs couplings that introduce sensitivity to the left-right mixing are negligible. A similar behavior can be found for the other cases considered in section 3.3 when holdingĈ h, f 1 /m f 1 fixed.
KeepingĈ h, τ 1 /m τ 1 and θ τ fixed implies that m τ 2 /m τ 1 and m ντ /m τ 1 vary with m τ 1 . This results in an increasing importance of co-annihilation effects with the heavier stau and the tau sneutrino. This effect is only significant forĈ h, τ 1 /m τ 1 ≤ 2, however. For Ĉ h, τ 1 /m τ 1 = 0.3, the co-annihilation effects in the extremely compressed stau sector lead to a net increase of the yield compared to the case of a single right-handed stau. Note that keeping the ratio m τ 2 /m τ 1 and θ τ constant would requireĈ h, τ 1 to increase proportionally to m 2 τ 1 , which would result in large deviations from eq. (26), in particular a net decrease of the yield with increasing m τ 1 . When raising m τ 1 in this setup the required large values forĈ h, τ 1 would quickly drive the model into phenomenologically unfeasible regions (see sections 4.4 and 5.1).

Summary and classification of regions
One important outcome of the performed survey is the fact that in all regions the scaling behavior of the yield with m τ 1 (and for fixed a i otherwise) is approximately linear. Hence, the desire to achieve low stau yields points to low stau masses and in the same way to low masses for co-annihilating sparticles close in mass to the stau. On the other hand, lower mass limits for the different sparticles can be derived from LHC searches. These searches potentially translate into lower limits on the stau yield for a given region. In order to be able to discuss the impact of these and further experimental and theoretical bounds on the yield in section 5, we will now summarize and classify the most important phenomenologically different regions we found in this section.
Without co-annihilation effects and significant left-right mixing in the stau sector the stau yield is roughly Y ≃ 10 −12 for m τ 1 = 1 TeV. The yield does not change orderof-magnitude-wise when introducing (nearly) mass-degenerate selectrons and smuonsslepton co-annihilation effects lead to a slight increase in the yield. In this region in parameter space the sleptons dominantly annihilate into vector bosons and leptons. We will refer to this region as the bulk region. The possibility of obtaining low stau yields in this region is simply restricted by the lower bound on the stau and slepton mass.
If we allow for EWino co-annihilation the effect on the yield ranges from a slight increase (for bino co-annihilation) to a decrease by almost an order of magnitude (for wino co-annihilation) with respect to the case of no co-annihilation. If we additionally approach the region of a resonant s-channel propagator, m A ≃ 2m χ 0 1 , we found a net decrease of the yield by up to more than two orders of magnitude for maximally mixed EWinos. The main limiting factor for achieving small yields in the EWino co-annihilation region is given by the lower limits on EWino masses in the long-lived stau scenario.
Similarly, lower bounds on the yield in the gluino co-annihilation region and squark co-annihilation region potentially arise from the respective bounds on the gluino and squark masses in the long-lived stau scenario. In both scenarios a reduction of the yield by around one order of magnitude could be achieved.
When considering scenarios where sfermions with large left-right mixing are involved in the annihilation processes, lower limits on the stau yield do not arise solely from the lower limits on sparticle masses. In contrast, the involved Higgs-sfermion couplings depend on a priori free parameters of the theory which can only be restricted by theoretical bounds (from vacuum stability and unitarity) or indirect experimental bounds (e.g., precision and flavor observables, MSSM Higgs searches). We distinguish two characteristic regions for which stau-Higgs couplings are important. The Higgs final state region and the Higgs resonant region are characterized by a dominant annihilation of staus into two final state Higgses (any combination of h, H, A 0 , H ± ) and annihilation into bottom or top quarks, respectively. Both effects are also present in the case of co-annihilation of stops or sbottoms with large left-right mixing which we summarize as the 3rd generation squark co-annihilation region.
In the following section we will describe the various constraints on the model parameter space whose implications for the seven regions defined here will be applied in section 5.

Implications of the first LHC runs
The LHC has brought important insights into the physics of elementary particles constraining possible extensions of the SM. In the first runs of proton-proton collisions at center-of-mass energies of 7 and 8 TeV, searches for the SM Higgs boson as well as for SUSY and other theories beyond the SM have been in the focus. Searches for additional Higgs bosons have imposed severe bounds on the MSSM Higgs sector, especially for the region of low m A . Furthermore, the determination of the couplings of the discovered Higgs boson to the SM particles may lead to indications for new physics that could serve as discriminators between different models beyond the SM.
We will here consider the latest results from the LHC experiments and further constraints and will work out their implications for scenarios with a long-lived stau. A key ingredient of the analysis is the interpretation of the searches for heavy stable charged particles (HSCP) performed at the 7 and 8 TeV LHC in the considered pMSSM parameter space. We will include not only the searches for charged sleptons but also the searches for R-hadrons, which can appear for small mass gaps between the gluino or squarks and the stau. Those regions are of particular interest for us in order to cover scenarios with gluino and squark co-annihilation. To interpret the collider bounds in the pMSSM parameter space, we have to compute the complete SUSY cross sections for each generated point in the parameter space. The enormous computing time for the calculation of the full SUSY cross sections at next-to-leading-order (NLO) precision makes it necessary to find other methods allowing for a fast estimation of the cross sections suitable for a large number of points. We achieved this by developing a fast cross section estimator based on grids and interpolation routines. This is particularly important for EWino production which in principle depends on many parameters.
In the present section, our main goal is to reveal the interplay between the constraints on the Higgs sector, the limits from HSCP searches, and other theoretical or experimental constraints from flavor and precision observables. To this end we will perform a Monte Carlo scan over the pMSSM parameter space.

Monte Carlo scan in the 17-dimensional parameter space
The pMSSM is based on the following assumptions on the general MSSM: (i) R-parity is conserved, (ii) all complex phases in the soft breaking potential vanish, so that no new sources of CP violation are introduced beyond contributions from the CKM matrix, (iii) sfermion mass matrices are diagonal in flavor space and the trilinear couplings are proportional to Yukawa couplings, so that no new sources of flavor violation are introduced, (iv) universality and vanishing trilinear couplings for the first and second generation sfermions are assumed. After imposing the electroweak symmetry breaking conditions, this leads to 19 free parameters. 10 To simplify the estimation of cross sections for collider bounds we further reduce the number of parameters by imposing which does not affect the qualitative discussion in the work. This way we are left with a 17-parameter pMSSM, with all parameters defined at the TeV scale. We impose the following hard restrictions on the generated points. First, the lighter stau is taken to be the NLSP (and thus the lightest sparticle of the MSSM), Second, at least one of the neutral CP -even Higgses lies within the LHC Higgs discovery window m h or/and m H ∈ [123; 128] GeV, which we will discuss in section 4.1.3 in more detail. The analysis performed here is independent of the nature of the LSP. We merely assume that a very weakly interacting non-MSSM sparticle is the LSP and that the lifetime of the stau is larger than O(10 −7 s), i.e., it is long-lived and thus leaves the LHC detectors before decaying. Note that this does not restrict our 17-dimensional parameter space since the LSP mass is an independent parameter which can easily ensure this constraint and does not have any further consequences for our analysis at this point. For the example of a gravitino LSP, τ τ 10 −7 s implies m G 0.4 MeV for m τ 1 = 1 TeV. Under these assumptions we perform a numerical random scan over pMSSM parameter space and generate points according to the following procedure.
1. After a random selection of the parameters at the low scale, we generate the physical masses as well as mixing angles using the spectrum generator SuSpect 2.41 [26]. The input parameters and scan ranges are described in section 4.1.1. Minimal requirements on the scan points are imposed on this stage-we only proceed with points obeying eq. (29) and the accepted output intervals of m τ 1 , m t 1 , m b 1 (see section 4.1.2). In the following, we describe 1-6 in detail.

Input parameters and scan ranges
As the 17 independent input parameters at the TeV scale we choose [max(m τ 1 , 1000); 5000] Table 1: Parameter ranges for the 17-dimensional pMSSM scan. The second column shows the intervals of the randomly generated input parameters. In the third generation sfermion sector we choose masses and mixing angles as input parameters and determine the corresponding soft masses from these input parameters at tree-level. The third column displays the accepted intervals for these masses and mixing angles after computing the full spectrum including higher order corrections. All dimensionful parameters are given in GeV. ⋆ The interval [0; π/2] is mapped onto [0; π/2] or [π; π/2] according to the sign of Xτ = Aτ − µ tan β, see section 4.1.2 for details. In order to avoid numerical instabilities we choose 10 −4 as a lower limit on scan range of the mixing angles.
sfermions, and treat them as independent input parameters. This has two reasons. First, this way we achieve a better control over the third-generation sfermion masses in the presence of large mixings and thus by choosing appropriate scan ranges we avoid scanning over regions where the stau is not the NLSP or which are already forbidden by conservative model-independent collider bounds. Second, spectra with large mixings are equally strongly represented as those with small mixings. This has an important impact on our considerations of stau yield which is potentially sensitive to the stau mixing angle.
If not stated otherwise, for all input parameters we choose linearly flat priors in the scan. The scan ranges are summarized in table 1. The ranges are motivated by the requirement of a τ NLSP as well as conservative collider bounds on individual particles (see section 4.1.2). First generation squarks and sleptons are kept degenerate. In addition to this 'blind' scan we performed dedicated scans accumulating more points in certain sub-ranges which are of particular interest according to the results of section 3. Those dedicated regions are summarized in appendix E. If not stated otherwise we refer to the complete set of scan points including the dedicated scans. In total we generated 5 × 10 5 points. Note that we do not attach any physical meaning to the density of generated points in parameter space.

Spectrum generation
After the random generation of the input parameters given in (31) we determine the soft masses of the third-generation sleptons and squarks, m 2 , from the respective free parameters in (31), using the tree-level relations (70) and (71) (see appendix B) and analogous expressions for stops and sbottoms. Points with negative mass squares are rejected at this point. From these input parameters the SUSY spectrum is computed with SuSpect 2.41. Points which do not fulfill (29) and (30) are rejected as well as points that do not lie within the accepted output intervals for m τ 1 , m t 1 and m b 1 listed in table 1. The lower limits of these intervals are motivated by conservative collider bounds on individual sparticle masses in the long-lived stau scenario we derived earlier (see [53,54]). However, we will see that they are well below the limits we will finally infer from the interpretation of the HSCP searches at 7 and 8 TeV. Hence, these lower limits only serve to gain efficiency in generating valid points and have no impact on the physical results.
SuSpect computes up to 2-loop corrections for the sparticle masses. For illustration, figure 10 shows the relative correction to the input parameters as a function of the output parameters computed by SuSpect. For the stau mass and mixing angle, loop corrections that are taken into account in the computation via SuSpect are relatively small. The bulk of points acquire corrections well below 10%. However, deviations up to 30% are present in the stau sector. For the stop and sbottom mass higher order corrections are much more significant. Especially in the case of the stop the output value for m t 1 turns out to overshoot the intended value by several 100%. However, as we only use the output values for all further discussions, within the limitations of SuSpect we achieve self-consistent spectra which we will use in the following discussion.
We re-compute the Higgs sector of the spectrum as well as the Higgs decay table with FeynHiggs 2.9.2. The value for the Higgs mass m h computed by FeynHiggs is smaller than the value computed by SuSpect for most of the parameter points. Since larger m h tend to be more challenging to achieve, we consider the lower limit on the Higgs mass to be more important. Thus, to be conservative we use the FeynHiggs value. The resulting spectrum is used for the further analysis. All points that fulfill (30) are recorded and count as generated points.

Meeting the LHC Higgs window
Within this work we interpret the discovered Higgs boson with a mass of 125.5 ± 0.2 (stat.) +0.5 −0.6 (syst.) GeV at ATLAS [55] and 125.7 ± 0.3 (stat.) ± 0.3 (syst.) GeV at CMS [56] as either the light or heavy neutral CP -even Higgs of the MSSM (or even as both contributing to the signal). Accordingly, taking into account the theoretical uncertainty in the prediction of the Higgs mass (see, e.g., [57]), we demand (30). The ranges of the input parameters have an effect on the distribution of the resulting Higgs masses. Through loop corrections, the sparticle masses, especially the stops, are intimately related to m h .
In a neutralino LSP scenario, a pMSSM scan with flat priors and input parameter ranges just above the current collider bounds, the distribution for m h typically peaks at values below the interval (30) and falls off over the interval towards large values (see, e.g., [58]). This implies that for the case of m h ∈ [123; 128] GeV this window would mostly be populated towards its lower end reflecting the preference of the MSSM for a lighter m h . In this work we aim to avoid an asymmetric distribution of m h around the experimental value since the allowed window is to account for the theoretical uncertainty in the computed Higgs mass. Instead, we choose to aim for a flat distribution in m h in our scan. Hence, we allow for relatively large A t in this scan. Remarkably, with the scan ranges given in table 10 we achieve an almost flat distribution in m h over the interval (30). This is partly due to the fact that in the long-lived stau scenario stronger modelindependent bounds on the sparticle masses exist which shifted our scan ranges towards higher masses (see section 4.1.2). The blue line in the upper panel of figure 11 shows the distribution of the Higgs mass m h for the blind scan (the distribution for the complete set of points is virtually identical).
A second effect on the Higgs sector is induced by the allowed range for m A . For the range chosen here 11 most parameter points end up in the decoupling limit avoiding to cover the region where m H could make up the discovered Higgs. In other words the ratio between the number of points with m h versus m H in the interval (30) depends strongly on the chosen scan range for m A . In order to have control over this arbitrary bias we require m A < 140 GeV for half of the generated points, i.e., half of the points in our scan lie explicitly not in the decoupling limit. This way, around 65% (35%) of the generated points feature m h (m H ) to lie in the interval (30). For around 0.7% of the points both Higgs bosons lie in this interval. (30) To obtain m h in the window (30) demands the presence of large radiative corrections on m h requiring an interplay of several parameters that govern these radiative corrections, namely the masses and the mixing in the stop sector and furthermore-in descending order of importance-in the sbottom and stau sector. While the stop contributions to the Higgs mass are large (as demanded) and positive, the sbottom and stau contributions typically diminish the Higgs mass and can be significant for negative µM 3 and large tan β, making it harder to satisfy (30) [59].

Selection effects induced by
These features induce a selection effect resulting in a non-flat distribution in some of the input parameters that we initially scanned over with flat priors. Although we do not assign any physical meaning to the absolute point density in the parameter space in our later results, it is, however, interesting to see in which way the flat priors are 'bent' by the additional requirement (30). This shall be subject of a brief discussion in this subsection. For this discussion we consider the 'blind' scan only.
The largest effect can be observed for A t , which shows a clear preference for large absolute values, |A t | > 3 TeV, according to the large mixing required in order to obtain high m h , see blue line in the right panel of figure 11. This effect is much less pronounced for those points where m H lies in the window (30) (blue curve). Further, the distributions in m b 1 and m t 1 are bent towards disfavoring the upper and lower part of the allowed scan range, respectively. Interestingly, if we restrict A t to a smaller range (e.g., |A t | < 3 TeV or less), the m t 1 distribution changes to favor the lower part of the scan range. This is due to the large (relative) mixing required and shows that this mixing is in fact more important than the overall stop mass scale. The maximal radiative correction is present Other parameters that are affected by the requirement of the Higgs mass and by the accepted output intervals for m τ 1 , m t 1 and m b 1 listed in table 1 are tan β, disfavoring values below ∼ 10, µ, peaking around ±2 TeV and the stop mixing angle, slightly disfavoring maximal left-right mixing, i.e., θ t ≃ π/4 or 3π/4. All

Interpretation of the HSCP searches in the pMSSM
Long-lived staus show up as heavy stable charged particles (HSCP) in the detectors at the LHC, i.e., they are recognized as muons but with two features that potentially allow for a discrimination against real muons: an anomalous time-of-flight (ToF) and an anomalous ionization loss (dE/dx). Both are accessible at the LHC experiments. So far, HSCP searches have been performed at ATLAS [60] (based on 4.7 fb −1 at 7 TeV) and CMS [14] (based on 5.0 fb −1 at 7 TeV and 18.8 fb −1 at 8 TeV) and no significant excess over background has been reported. The null searches have been interpreted in a few long-lived stau scenarios: for a GMSB scenario (ATLAS and CMS) as well for direct production of mass-degenerate sleptons (ATLAS) and the direct production of staus only (CMS). The latter analysis provides an almost model-independent lower bound on the stau mass of 339 GeV. 12 We will here interpret the recent search of CMS in the framework of the 17-dimensional pMSSM. To do so, we determine the cross sections for all relevant SUSY production processes for each scan point, as described in section 4.2.1. The estimation of the cross section upper limit extracted from the search [14] will be described in section 4.2.2.

Fast estimation of SUSY cross sections
For each pMSSM point we determine the cross sections for various production channels at the 7, 8 and 14 TeV LHC in order to estimate the viability of each point after the HSCP null-searches and to discuss the prospects for the LHC long-term run. The computation of all potentially relevant SUSY cross sections at NLO precision is time-consuming 13 and especially not convenient for the use in Monte Carlo scans containing a large number of points. In order to achieve a sufficiently fast determination of the cross sections for each generated pMSSM point we develop a fast cross section estimation tool based on grids and interpolation routines. However, some production processes in principle involve many parameters requiring high-dimensional grids, which would mean to shift the problem of large computing time to the generation of the grids. Therefore, we exploit the potential for approximations wherever suitable. By factorizing the dependence on certain combinations of parameters we describe all channels approximately with a set of up to maximally three-dimensional grids. In the following we will list the respective parameterizations and approximations chosen in the different sectors.

Slepton sector
In the slepton sector we build up one-and two-dimensional grids in the corresponding sparticle masses for the processes e R e R , e L e L , ν e ν e , τ 1 τ 1 , τ 2 τ 2 , ν τ ν τ and e L ν e , τ 1 τ 2 , τ 1 ν τ , τ 2 ν τ , respectively. For this purpose we compute the cross section for Drell-Yan (DY) production (via an s-channel γ/Z or W ± ) with Prospino [46] at NLO. SUSY QCD contributions have been kept small by setting the mass of all colored sparticles to 5 TeV. For the third-generation sleptons the left-right mixing introduces an additional variable the cross section depends upon. Here, we make use of the fact that the dependence on the stau mixing angle θ τ factorizes once the center-of-mass energy √ŝ of the production process is well above M Z [62]. This limit is easily reached for the rather heavy stau masses we are considering. To be concrete, the cross section for the process pp → ij, i, j = τ 1 , τ 2 , ν τ , can be written in the form We choose where m ref = 500 GeV. For the evaluation of the cross section in the scan we interpolate logarithmically over the cross section A and linearly in the correction factor B.

EWino sector
In the EWino sector we parametrize the cross sections by the underlying SUSY input parameters instead of the physical masses and mixings, namely M 1 , M 2 , µ, tan β and the common first-and second-generation squark soft mass m Q 1,2 . In order to describe the cross section as a function of these five parameters with maximally three-dimensional grids we factorize the dependence on these five parameters as follows. First, we decouple the bino from the spectrum and consider M 1 separately from M 2 and µ. This is motivated by the hierarchy in the respective couplings: for degenerate M 1 and M 2 (or µ) the bino contribution is relatively small. Second, we factorize the dependence on the squark masses. This dependence is introduced by t-channel squark diagrams which can lead to a significant net reduction of the cross section even though taking squark pair and associated squark production into account. This arises from a negative interference between the DY production of EWinos and the t-channel contribution and is relevant for intermediate mass gaps between m Q 1,2 and M 2 where the squarks are still light enough to contribute in the t-channel but already too heavy to (over-) compensate the reduction by squark production.
We found that the complete cross section from neutralino and chargino production (including the associated squark-EWino production) can be well approximated by three functions, each depending on three parameters: where χ i = χ 0 1 , χ 0 2 , χ 0 3 , χ ± 1 , χ ± 2 , and where χ i q denotes the associated EWino-squark production. Where they are not displayed as an argument in the brackets, we set mass parameters to 5 TeV and tan β = 15. We computed the three grids, corresponding to the three functions in eq. (35), with Prospino [46] at NLO precision. To save computing time we only ran the NLO computation for a subset of points and extracted the K-factor from the resulting coarser grid, under the assumption that the K-factor varies more slowly with varying parameters than the cross section itself.
The functions have the weakest dependence on the last argument in each of the brackets in eq. (35). Accordingly, we computed significantly fewer points in the corresponding directions in the grid space. Contributions from associated gluino-EWino production were neglected. For the generation of the spectrum from the SUSY parameters we used SuSpect 2.41, as we did for the generation of the pMSSM points in the Monte Carlo scan. We interpolated logarithmically over the cross sections and linearly in the correction factor R as well as in the K-factors. With this description we found an agreement within a 15% error with the full NLO computation with Prospino for a variety of very different spectra.

Squark and gluino sector
For the production of third-generation squarks, contributions from t-channel gluino diagrams are small due to the small parton densities of the required heavy-flavor quarks.
Furthermore, electroweak production is relatively unimportant. Hence, the relevant production channels are t 1 t 1 , t 2 t 2 , b 1 b 1 , b 2 b 2 via an s-channel gluon diagram, a t-channel squark diagram or the gluon-squark four-vertex. The production cross sections for these processes only depend on the mass of the respective squark alone. For the 7 and 8 TeV LHC cross sections we take the corresponding one-dimensional grids from NLLfast [50] which include NLO and next-to-leading-log (NLL) corrections. For the 14 TeV case we compute the grid with Prospino [48] at NLO.
For the first-and second-generation squark and gluino production, g g, q g, q q and q q * , we interpolate two-dimensional grids in the variables m g and m q ≡ (m u L m u R m d L m d R ) 1/4 which are taken from NLLfast [49] for the case of 7 and 8 TeV LHC cross sections and which we compute with Prospino [45] at NLO precision for the 14 TeV LHC cross sections. We interpolate logarithmically over the cross sections. The error from the interpolation is typically less than 1%.
The total cross section obtained from summing over all the processes described above was compared to the full cross section from Prospino for a variety of different spectra and found to agree within an error of typically 10%. For a few points we found errors up to 15% where we underestimate the cross section computed by Prospino.

Stau production via intermediate Higgs
In addition to the above channels we include the direct production of staus via an schannel Higgs intermediate state [63]. The channel pp → h → τ 1 τ 1 can be important in the presence of large left-right mixing of the stau. Additionally, we take into account the heavy Higgs intermediate state pp → H → τ 1 τ 1 . As mentioned earlier, for the general case (no decoupling limit) these processes depend on a variety of parameters. Accordingly, we compute the production cross section for these channels for each of the generated pMSSM points using the complete spectrum. We perform the computation at the leading order via Whizard 2.1.1 [43] where the effective gluon fusion vertex for the MSSM [64] has been implemented. We consider gluon-fusion and bottom-fusion. For the production via bottom-fusion we reweight the cross section according to the resummed bottom-Higgs coupling (for the leading contributions in tan β), as described in appendix D. For this computation we employed the value for the correction to the bottom mass, ∆ b , from micrOMEGAs.

Estimation of cross section upper limits
As shown in [54], the signal efficiency 14 for the signatures of long-lived stau scenarios at the LHC is much less sensitive to the spectrum than, e.g., in the case for scenarios with neutral stable sparticles escaping the detector, where compressed or widely spread spectra are typically much harder to find. In this reference it has been shown that for the production via colored sparticles the signal efficiency of long-lived staus only drops below roughly 20% for widely spread spectra for which this production mechanism is no longer the dominant channel but is exceeded by the direct production of staus which provides higher signal efficiencies. This way, the signal efficiency for the total SUSY production does not drop below about 20% in the mass ranges of interest for the LHC analysis, provided that there is no long-lived sparticle other than the stau and thus all decay chains terminate in the stau before traversing the sensitive parts of the detector. The significant decrease of the signal efficiency for the production via colored sparticles for widely spread spectra is due to the potentially large boost of the stau developed in the decay of a very heavy colored sparticle. Staus with a velocity close to the speed of light, β ≃ 1, are extremely difficult to discriminate against background muons since the discrimination heavily relies on a deviation from β = 1.
Following this argument, electroweak production mechanisms, e.g., chargino production, offer even less potential to cause a drop in the overall signal efficiency. This is because, due to the smaller electroweak cross sections, the mass gap between the produced sparticle and the stau is smaller if the electroweak production process in question is demanded to give a significant contribution compared to the direct stau production. This fact facilitates the estimation of the signal efficiencies (and for the resulting cross section upper limits) requiring the extrapolation of the results given in [14] to a general pMSSM point. In the following we will describe this procedure in more detail.
If the decay of heavier sparticles into the stau is not prompt, the analysis becomes more complicated. We will examine the case of long-lived colored sparticles which we found to be the most relevant in this study. In particular, gluinos can become long-lived even for relatively large mass gaps m g − m τ 1 100 GeV. The treatment of long-lived colored sparticles is described below.

Application for prompt decays into the stau
We consider a point to be excluded at 95% C.L. if the signal strength, σ limit /σ th , obeys where σ limit is the observed 95 % C.L. upper cross section limit from the experiment and σ th is the theoretical prediction for the total cross section. σ limit is a model-dependent quantity. In the simplest case, for a given spectrum, the upper cross section limit is determined by where S is the required number of expected signal events for the considered spectrum which allows for a 95% C.L. exclusion in the presence of the observed number of (background) events. ε S is the signal efficiency for this spectrum and L is the integrated luminosity. S and ε S both are affected by the applied cuts-the latter directly and the former via its background rejection capability. In HSCP searches the highest sensitivities are typically reached for cuts that supply S = 3 for a 95% C.L. exclusion [54].
In the CMS analysis [14] the observed upper cross section limits are given for the two benchmark models (GMSB model and direct DY production) for the 7 and 8 TeV run as a function of the stau mass σ limit (m τ 1 ). Here, we take the combined Tracker+ToF data. In order to estimate the signal strength for a point in our pMSSM parameter space we assign the upper cross section limits channel-wise: For the direct DY production of the lighter staus we apply the direct DY production cross section limits. For all other slepton production mechanisms, the EWino production and the production of third-generation squarks we applied the cross section limits from the GMSB model as a function of the stau mass. This is done under the assumption that the signal efficiencies and corresponding background rejection for these channels are similar to the GMSB model, which is based on the arguments given above. 15 For an arbitrary stau mass we interpolated linearly between the analysis points given in [14]. For stau masses above 500 GeV we will only be in the vicinity of the exclusion limit if we have a rather degenerate spectrum and thus an important strong production of sparticles. For these production modes the signal efficiency can decrease due to difficulties in the triggering of very slow staus [54]. In order to account for these spectra we extrapolated the upper cross section limits by conservatively assuming σ limit = 3.0 fb for the 7 TeV run and σ limit = 1.0 fb for the 8 TeV run. These values are in accordance with the signal efficiencies that have been reported in [54] in the limit of mass-degenerate spectra where one stau has been required to have a velocity above β = 0.6 in order to ensure an efficient triggering of such events. For the production of staus via first-and second-generation squarks and gluinos as well as for the direct production via an s-channel Higgs we take as a conservative estimate a constant σ limit = 3.0 fb (1.0 fb) for the 7 TeV (8 TeV) run. 16 The signal strength is then obtained by where σ th ik is the computed cross section for the channel i at the LHC energy k and σ limit ik is the corresponding estimated observed cross section upper limit for the respective channel.

Application to delayed decays
For the application of collider limits to the present scenario, it is crucial to know if there are long-lived sparticles other than the stau which play a role in the production and decay at the collider. We therefore compute the width of all sparticles. We used a modified version of SDecay [65,42] which includes all relevant 3-body decays of sleptons into the lighter stau. We compute further 3-and 4-body decays of squarks and gluinos into the stau, relevant if m q < m χ 0 1 and m g < m q , m χ 0 1 , with Whizard 2.1.1 [43]. 15 For the GMSB model considered in [14] (m τ 1 = 308 GeV) the EWino production contributes 53% while the direct DY production of the lighter stau and all other sleptons make up 13% and 33% of the total SUSY cross section, respectively. The contribution from first-and second-generation squarks is negligible. 16 For the direct production of staus via an s-channel neutral, CP -even Higgs (h/H), stau production near threshold is enhanced and so the fraction of very slow staus is large. For this channel the decreasing trigger efficiencies for smaller velocities (below ∼ 0.6) are expected to be the restricting factor of the Figure 12: Scatter plot displaying the dependence of the width of the squarks (blue points) and gluino (red points) on the absolute mass difference to the stau, ∆m = m q −m τ 1 and m g −m τ 1 , respectively. We consider all squarks here, including stops and sbottoms. We only plot a point if the corresponding width of the squark or gluino is the smallest width. The horizontal lines correspond to Γ g, q = 2 × 10 −14 GeV and Γ g, q = 2 × 10 −16 GeV. Figure 12 shows the mass gap between the squarks and the stau (blue points) as well as the gluino and the stau (red points) versus the resulting decay width of the respective sparticles. We only plot the points for which this width is the smallest among all widths of sparticles heavier than the stau. (This also ensures that parameter space points where both a squark and the gluino are long-lived appear only once.) For the gluino, even for mass gaps of up to 300 GeV we encountered points with small gluino widths which imply non-prompt decays into the stau. Note that these situation can only appear in the case that the masses of the squarks and EWinos are well above the gluino mass such that the 4-body decays are suppressed by two off-shell propagators. For other situations the gluino width is typically much larger. We do not take into account loop-induced decay modes of gluino and squarks into staus leaving this for future investigations.
In order for the tracker analysis (dE/dx) to be efficient, the longitudinal and transversal impact parameter of the track candidates, d z and d xy , are required to be smaller than 0.5 cm [14]. Bearing in mind that non-prompt decays typically play a role in the case of rather small relative mass gaps between the heavier mother sparticle and the stau we do not expect a very pronounced kink in the track. We therefore consider a mother sparticle X to be sufficiently short-lived to allow for the daughter stau to pass the tracker requirement, if This corresponds to a decay length of cτ X < 1 cm. 17 For neutralinos and sneutrinos it requires very small mass gaps in order to violate eq. (40). Consequently, these cases appear very rarely in our scan-0.15% of the points signal efficiency. A detailed study of the signal efficiency in this channel is left for future work. 17 The decay length for a relativistic particle X is cβγτX . However, βγ ≃ 1 for β ≃ 0.7. For heavy colored sparticles produced close to threshold, β 0.7 is a typical velocity.
contain metastable neutralinos while 0.6% of the points contain metastable sneutrinos. The determination of the appropriate collider limits for these cases requires a detailed analysis of all branching fraction and the consideration of various missing energy searches. Since these points are not of particular interest for this work we will leave the investigation of these cases for future work and will simply reject the corresponding points from the scan. For metastable charged sleptons other than the stau as well as metastable charginos, we expect the analysis to be virtually identical, regardless whether they decay into the stau or not, assuming that a possible kink in the track will not significant change the sensitivity to the signature.
The case of metastable squarks and gluinos appears more frequently in our scan. We found that 5.8% and 6.7% of the points contain metastable squarks and gluinos, respectively. On the one hand, this relatively large fraction arises from the suppression of the required 3-and 4-body decays, on the other hand, it results from the dedicated scans, specifically accumulating points in the corresponding mass-degenerate regions (see table 2). In the following we describe the treatment of metastable squarks and gluinos in the determination of the cross section upper limits.
If a metastable squark or gluino decays delayed, Γ g, q < 2 × 10 −14 GeV, the stau is assumed not to be recognized in the tracker. Consequently, we only apply the ToF analysis taking into account the data from the muon chambers only. We refer to this data as the 'muon-only' analysis in the following. The cross section upper limits for the muon-only analysis have been reported for stops and gluinos only, where the direct production of these sparticles is taken to be the only production mechanism. If we apply the muon-only analysis on long-lived staus we have to assume that the kinematics of the staus are similar to the strongly interacting mother sparticles that dominate the production. This is indeed the case for the small mass gaps that are required to cause the delayed decay of the stau. Furthermore, the detector response of the drift-tubes in the muon chambers to an R-hadron carrying one unit of electric charge is virtually the same as for long-lived staus. Hence, we estimate the cross section upper limits for staus in the muon-only analysis by the limits derived for stops. Note that the muon-only analysis has only be performed for the 8 TeV LHC run.
If the metastable colored sparticle has an even smaller decay width, Γ g, q < 2 × 10 −16 GeV, corresponding to cτ g, q 1 m, the muon-only analysis might not be applicable anymore. We therefore assume in this case that the strongest sensitivity arises from the R-hadron itself that is recognized in the tracker. 18 Consequently, we apply the cross section upper limits from the corresponding R-hadron search where we conservatively choose the charge suppression model for the gluinos and squarks. To all production processes whose decay chains terminate in late decaying staus seen in the muon-only analysis or in R-hadron searches containing a gluino or squark the respective cross section limits are applied. By doing so, we implicitly assume that production modes of sparticles are only relevant if the mass gap between the produced sparticle to the respective sparticle seen in the detector is small or that the corresponding signal efficiencies do not depend strongly on the mass of the produced sparticles. The final signal strength is then determined by eq. (39). (For those production processes that lead to a prompt decay into the stau we employ the Tracker+ToF analysis as described above.) The interpretation of the HSCP searches leads to very restrictive bounds on the sparticle masses. For example, we did not find any allowed point in our scan with m t 1 , m b 1 850 GeV, m q 1400 GeV and m g 1200 GeV. Regarding the EWino sector, no point with |µ|, M 2 800 GeV survived the bounds.

Further experimental constraints
In this section we will discuss the implications of the most important experimental and theoretical constraints on the considered 17-parameter pMSSM beyond direct SUSY searches considered in section 4.2.

Constraints from Higgs searches at colliders
In addition to the condition (30) we require that the scan points pass a variety of collider bounds from the Higgs searches at LEP, the Tevatron and the LHC imposed at the 95% C.L. For the application of these bounds we use the program package Higgs-Bounds 4.0.0 [66], which tests the compatibility of the predictions for the Higgs sector in a given model against Higgs rates and masses measured in the mentioned experiments. We employed the full set of experimental results supplied by HiggsBounds. For the predictions for the spectrum of the MSSM Higgs sector HiggsBounds is linked to FeynHiggs 2.9.2.
The constraints have a large effect on our parameter space. Most importantly, the bounds depend on m A . Generically, we find that the parameter space is constrained much more strongly for smaller values of m A . Accordingly, in the subset of points where the heavier CP-even Higgs takes the role of the SM-like Higgs, i.e., where 123 GeV < m H < 128 GeV, nearly all points (99.88%) were rejected by the application of HiggsBounds. Most of these points (around 98%) were rejected 19 by the CMS search for MSSM Higgs decays into tau pairs (h, H, A 0 ) → τ τ [67]. The majority of the remaining points were excluded by the search for Higgsstrahlung processes at LEP, where the Higgs is assumed to decay into bb, (h, H, A 0 )Z → (bb)Z [68]. Other processes are less important.
In the subset of points where the lighter CP-even Higgs plays the role of the SM-like Higgs, i.e., where 123 GeV < m h < 128 GeV, around 27% of the points were excluded. Again, for most of these (around 91%) the CMS search for (h, H, A 0 ) → τ τ provides the highest significance. Further analyses of high importance are the search for (h, H, A 0 ) → ZZ → ℓℓℓℓ at CMS [69] and searches for a charged Higgs at CMS [70]. Figure 13 shows the allowed (green) and rejected points (blue, yellow and red) in the m A -tan β plane. Considering the rate of allowed versus rejected points in the different regions, the decoupling limit appears to be strongly favored by the current data.

Constraints from flavor and precision observables
Supersymmetric corrections to the mass of the W boson impose another constraint on the parameter space. Here, we use the experimental value M W = (80.385 ± 0.015) GeV [71]. Following [72,73], we increase the uncertainty by a theory error of 15 MeV, combine the uncertainties linearly and multiply them by a factor of two in order to estimate the allowed range at the 95% C.L. Thus, we apply the limit M W ∈ [80.325; 80.445] GeV (41) to the value calculated by FeynHiggs 2.9.2. The flavor observables BR(B → X s γ) and BR(B 0 s → µ + µ − ) can be directly obtained from micrOMEGAs. We use the world average BR(B → X s γ) = (3.43 ± 0.21 ± 0.07) × 10 −4 [74]. Treating the uncertainties as above we find the allowed range at the 95% C.L.: The rare B 0 s decay has been observed with a branching ratio in the 95% C.L. range [75,76] BR Figure 14 illustrates the impact of these limits on the considered pMSSM parameter space. The limit on M W rejects the largest number of points. The lower panel shows that our choice (41) ensures that the deviation of the ρ-parameter from its SM value, ∆ρ, does not exceed 0.0018. The limit from B → X s γ is particularly restrictive for the subset of points with m H in the LHC Higgs window as given by (30). Both flavor constraints imposed here favor large m A .

Bounds from charge or color breaking minima
For large values of certain parameters, the MSSM scalar potential can acquire minima where U(1) em or SU(3) c is broken (charge or color breaking, CCB). For large tan β, requiring the standard electroweak vacuum to be stable or metastable with a lifetime larger than the age of the universe implies an upper bound on the product µ tan β [77,78,79,80]. We use [80], 0 < − |µ tan β eff | + 56.9 m L 3 m e 3 + 57.1 m L 3 + 1.03 m e 3 − 1.28 × 10 4 GeV where and I(a, b, c) is defined in eq. (85) in appendix D. The quantity ∆ τ describes the higherorder corrections to the tau Yukawa coupling in the limit of large tan β, analogous to the corrections to the bottom Yukawa coupling discussed in appendix D. Depending on m A and A τ , the upper bound (44) can become more stringent by approximately 20% [79].
In order to take into account CCB constraints on the trilinear couplings, we apply the simple conditions [81,82,83,84,85] We caution that the listed analytical constraints are not always reliable [86,87,88,89,90]. We impose them as a conservative first estimate, leaving a detailed numerical analysis employing the recently released program Vevacious [91] for future work. Figure 15 shows the impact of the constraints (44) and (47)(48)(49) on the considered pMSSM parameter space. The bounds on the trilinear couplings are quite restrictive. Furthermore, we see that the chosen range for A τ almost saturates the allowed region.

Stau yields in the Monte Carlo scan
The results of section 3 allowed us to identify all regions that potentially lead to exceptionally small stau yields. In this section we will investigate the limiting factors for low stau yields that could arise from various constraints. This is especially important for regions that contain large Higgs-sfermion couplings which are governed by a priori free parameters of the theory. In the presence of large left-right mixings of the sfermions one can only constrain the possible values of the yield by imposing constraints on the parameters that govern the Higgs-sfermion couplings. Working out the impact of these constraints is the subject of the present section. Furthermore, we will quantify how HSCP searches constrain the possible values of the yield. These searches are especially constraining in the case of co-annihilation with colored sparticles. Therefore, we will utilize the pMSSM Monte Carlo scan introduced in section 4.

Application of constraints
The upper panels of figure 16 show the effect of the constraints discussed in section 4. The blue points are rejected by the HSCP searches performed at the 7 and 8 TeV LHC (see section 4.2 for details). The most obvious result is that the HSCP searches reject all points with m τ 1 340 GeV. This is the most conservative bound on the stau mass in agreement with the bound reported in [14]. 20 For small stau yields the bound on the stau mass tends to become more restrictive-the border between blue and yellow points shows a kink at around Y = 10 −16 . This feature can be understood as follows.
In the region of small stau masses, small yields Y 10 −15 are typically achieved in the Higgs final state region (green points in the middle left panel in figure 16) where the couplings to the Higgs are enhanced. For these points the production of staus via a light or heavy CP -even neutral Higgs at the LHC is typically the dominant contribution to the stau production (see green and red points, respectively, in the lower left corner of the lower left panel in figure 16). This additional production mode raises the stau mass limit and forbids this region. Here we see a first correlation between the observable in the early universe and the measurements at the LHC. A similar effect occurs in the Higgs resonant region. This is best seen in the right panels of figure 16, where we plot the yield against m H /m τ 1 . In the resonance peak, m H /m τ 1 ≃ 2, very small stau yields are obtained. However, the very tip of this peak is excluded by HSCP searches, to a large extent due to the resonant production of staus via the heavy Higgs (see lower right panel of figure 16). For co-annihilation scenarios the bounds on the sparticle masses restrict the possible stau yields according to the scaling of the yield with the stau mass. The yellow points in the middle left panel of figure 16 show the domain of the co-annihilation regions in the m τ 1 -Y plane.
The bounds from MSSM Higgs searches taken from HiggsBounds and the flavor and precision bounds (abbreviated by FP in the following) are particularly restrictive in the region of small m H /m τ 1 . The yellow points in the upper panels of figure 16 are rejected by HiggsBounds and FP constraints. In fact the complete region of m H /m τ 1 0.2 is excluded by these bounds. For smaller yields Y 10 −13 , even higher values of m H /m τ 1 20 The bound in [14] was obtained for an almost completely right-handed lighter stau. As a slightly smaller DY cross section for τ1 production is obtained for θ τ = π/2 [62], we could expect allowed points lying O(10) GeV below the limit of [14]. However, such a value for θ τ requires either the heavier stau and tau sneutrino to be relatively light or the off-diagonal elements of the stau mass matrix, Xτ , to be relatively large leading to an enhanced stau-Higgs coupling. In both cases additional contributions enhance the overall production rate. Hence, we do not find any allowed points below the limit of [14].
Y Dominant production ch. mH/m τ 1 Y Dominant production ch. Figure 16: Distribution of scan points in the m τ 1 -Y plane (left panels) and mH /m τ 1 -Y plane (right panels). Upper panels: Effect of the constraints on the parameter space. The blue, yellow and red points are rejected by the HSCP searches, HiggsBounds+FP constraints and CCB bounds, respectively. The green points pass all the constraints. Middle panels: Dominant annihilation channels. The red, green and yellow points belong to the Higgs resonant region, Higgs final state region and co-annihilation regions, respectively. The blue points do not belong to one of these classes. Lower panels: Production channels that contribute dominantly to the strength of the HSCP signal. For the green and red points direct stau production via a the light and heavy Higgs is dominant, respectively. The yellow points are dominated by other production processes in the stau sector. The blue points are dominated by other processes. Note that the point density is saturated in parts of the plane such that blue points are simply covered by the others, etc.
are rejected by the HiggsBounds and FP constraints. This is partly due to the fact that the regions with smaller yields Y 10 −13 are dominated by the Higgs final state region and Higgs resonant region (see green and red points, respectively, in the middle panels of figure 16), which require large stau-Higgs couplings. These are more easily achieved with large values of tan β for which the constraints on m A from HiggsBounds become even stronger, see figure 13.
Constraints from CCB reject points in all corners of the displayed planes. The constraints on A t and A b can affect points without co-annihilation effects with stops or sbottoms and are therefore not necessarily related to the stau yield. However, a clear correlation is seen in the region of smallest stau yields. The CCB bounds push up the minimal yield allowed by the HSCP bounds by about another order of magnitude, see red points in the upper right panel of figure 16.
Note, finally, that the allowed points (green points in the upper panels of figure 16) all lie within a relatively narrow band in m H /m τ 1 . They span about four orders of magnitude in the yield, 2 × 10 −16 Y 4 × 10 −12 . In figures 17 and 18 we show the effect of the constraints on the parameter space for the above defined regions separately. The red points belong to the respective region, while the blue points belong to the complete set of points. The pure colors denote the allowed points, while the pale points are excluded by one or more of the constraints. Points with yields smaller than 10 −14 occur only for resonant annihilation via a heavy Higgs (see middle panels of figure 17). Among these points, the smallest yields (Y 10 −15 ) are achieved for dominant stau annihilation and no co-annihilation effects. Away from the heavy Higgs resonance we find yields as small as 2 × 10 −14 in the Higgs final state region with large stau-Higgs couplingĈ h, τ 1 ∼ 1 (cp. middle left and upper right panels of figure 17), and slightly below 10 −13 in the gluino and 3rd generation squark co-annihilation regions (see upper left and middle right panels of figure 18).
It is interesting to note that in the EWino co-annihilation region (lower panels of figure 17) and in the 3rd generation squark co-annihilation region (middle and lower panels of figure 18) stau yields down to roughly 5 × 10 −15 are allowed. The smallest yields are again reached for resonant annihilation via a heavy Higgs. In these regions no particular left-right mixing in the stau sector (and for EWino co-annihilation no particular left-right mixing in the sfermion sector at all) is required. Hence, these are the lowest values we found that could equally be realized in scenarios with a selectron or smuon NLSP.
The points with the largest yields almost always belong to the bulk region (see blue points in the middle panels of figure 16). Note that there is a relatively sharp limit of existing points in the high yield end, in contrast to the lower end of the range featuring a few scattered points with very low yields. This is due to the fact that the potential to increase the yield is limited by the number of sparticles that could increase the yield by virtue of co-annihilation effects. In fact, the estimate given in eq. (12) lies approximately in the middle of the band of blue points (bulk region) in the middle left panel of figure 16. Thus, eq. (12) is not too far from the largest yields that can be achieved in the pMSSM.
The percentage of surviving points in the regions is 4.4% in the bulk region, 0.18% in the Higgs final state region, 5.2% in the Higgs resonant region, 5.8% in the EWino co-annihilation region, 1.1% in the gluino co-annihilation region, 3% in the squark coannihilation region and 3.7% in the 3rd generation squark co-annihilation region. We plot the yield against the rescaled Higgs-sfermion couplingĈ Φ, f /m f for the case of the stau as well as for the case of the stop and sbottom in the upper and middle right panels of figure 17 as well as the lower left and right panels of figure 18, respectively. In the latter case we exemplarily plotĈ h, t 1 /m t 1 andĈ H, b 1 /m b 1 ; the couplings to the respective other Higgs behave roughly similarly. Large values are typically excluded mainly by the CCB bounds and precision observables as well as by flavor constraints.
Finally, we note that unitarity of the S-matrix sets further bounds on the involved couplings, see, e.g., [92,13,25]. The minimal yields allowed by unitarity are roughly Y ≃ 7 × 10 −18 (m τ 1 / TeV) [25] for stau-stau annihilation and Y ≃ 4 × 10 −17 (m τ 1 / TeV) for third-generation squark co-annihilation, taking additional degrees of freedom into account. As the minimal yields allowed in the respective regions are more than an order of magnitude larger than these values we assume that the bounds from the requirement of unitarity are significantly weaker than the other bounds considered in this paper, especially those from CCB minima. However, a detailed analysis investigating the particular annihilation processes relevant for our scenario and the effects of relaxing the approximations that were used in [25] appears worthwhile and may lead to more stringent bounds.

Prospects to narrow down the stau yield at the LHC
In the case of a discovery at the upcoming LHC runs it would be desirable to determine the stau yield from the LHC data and conclude on the viability of the underlying cosmological model. This is a difficult task as the yield depends upon various parameters with very different accessibility at the LHC. As a first step, we discuss in this subsection how one might be able to determine the parameter space region the scenario belongs to, which allows to narrow down the allowed range for the stau yield. The discussion remains qualitative and is not intended to be exhaustive.
The points in the scan that are close to the exclusion limit from the HSCP searches at 7 and 8 TeV typically provide a SUSY cross section at the 14 TeV LHC run of σ SUSY 14 TeV ≃ 100 fb. This gives us a rough idea of the prospects for studying long-lived stau scenarios at the LHC. For instance, with 300 fb −1 we obtain a total amount of 3 × 10 4 SUSY events for these points. In fact, due to the prominent signature of staus at the LHC, we could already learn a lot about the spectrum from much fewer events. First, already at the stage of discovering a long-lived stau scenario by the measurement of charged highly ionizing tracks in the detector, we are provided with a good determination of the stau mass with a precision around 15% [61]. In the search for long-lived staus, discovery is expected to take place on the basis of very few observed signal events [54], which translates into a total amount of O(10) produced stau pairs. 21 Second, the cross section for direct stau production differs from that for the production of colored sparticles 21 The discovery reach for an integrated luminosity of 300 fb −1 at the 14 TeV LHC is m τ 1 ≃ 700 GeV for the most conservative case of minimal direct DY production and up to m τ 1 ≃ 3 TeV for the case where the stau, the gluino and the squarks are close in mass [54]. The exclusion reach is similar.     with a similar mass by around five orders of magnitude. Indeed, in our scan the SUSY production cross section for a given stau mass spans four to five orders of magnitude, where the lower edge corresponds to points with dominant direct DY production while the upper edge corresponds to scenarios with a very small mass splitting between the staus and the colored sparticles, in particular the first and second generation squarks. Thus, from the relatively precise determination of the stau mass and a rough idea of the production cross section one might, already at the stage of discovery, be able to decide whether the data is compatible with a gluino or squark co-annihilation scenario or not. Namely, if the stau is relatively light such that the number of observed events is compatible with direct DY production, a co-annihilation scenario that could provide low yields could be excluded. On the other hand, if the stau is relatively heavy with respect to the measured rate of events such that dominant direct DY production is excluded, there are a variety of possibilities that could apply. This is in particular true for the intermediate range of production rates which could be compatible with stop or EWino co-annihilation or resonant stau annihilation via a heavy Higgs. In this case, more data is needed to distinguish between different scenarios. Let us briefly comment on two of them.
The first case concerns closely mass-degenerate staus and gluinos or squarks. As shown in section 4.2.2, here the appearance of delayed decays is a quite common feature, at least in the absence of other nearly mass-degenerate sparticles. Provided a very good understanding of the detector, such a scenario could hence be identified by the appearance of charge flipping tracks or other peculiarities that could occur due to the presence of long-lived or late-decaying R-hadrons in the detector.
Another scenario, which is particularly interesting, is the Higgs resonant region. Due to the appearance of the equally resonant production channel at the LHC this scenario provides a distinct signature [63]. We have seen from the lower right panel in figure 16 that this production channel can indeed be the dominant production channel of staus at the LHC particularly in the region of low stau yields. As discussed in [63], the velocity distribution of staus arising from the s-channel Higgs diagram peaks at significantly lower velocities than, for instance, that for direct DY production. Although challenging for the trigger settings (see, e.g., [54]), this signature can provide a way to distinguish the resonant s-channel Higgs region from other regions. Furthermore, the invariant mass of these events would reveal a distinct peak at twice the stau mass once more data is accumulated. Note that this signal is quite clean with minimal dilution by background. Consequently, such a peak might be visible with a comparatively small number of events.

Conclusions
In this work we presented a thorough survey for possible values for the stau yield in the framework of the MSSM with a long-lived stau NLSP. Focussing on the mass region that might still be accessible to a discovery at a long-term LHC run at 14 TeV, we pinned down the various possibilities for obtaining small stau yields in the pMSSM parameter space. In particular we showed the different possibilities to lower the stau yield by co-annihilation effects, resonance effects, enhanced Higgs-sfermion couplings and combinations thereof. We were able to determine the following configurations with an increasing potential to achieve low stau yields. In the absence of any left-right mixing in the stau sector a light neutralino in the t-channel of the annihilation diagram can lead to a decrease in the yield with respect to the decoupled neutralino case, typically by a factor of about 2.
In contrast, a co-annihilating bino as well as co-annihilating first-and second-generation sleptons increase the yield, again by factors of roughly 2. Scenarios with squark and gluino co-annihilation can lead to a decrease of the yield by a factor of O(10). We found that a decrease of the yield by significantly more than one order of magnitude can only be achieved through annihilation processes which involve large Higgs-sfermion couplings, a resonant Higgs in the s-channel or both.
In order to evaluate the phenomenological viability of the considered parameter space regions we performed a Monte Carlo scan over the 17-dimensional pMSSM with the stau being the lightest among the MSSM sparticles. We interpreted the Higgs boson recently discovered at the LHC as one of the CP -even neutral Higgses of the MSSM. By restricting m A to small values we forced around half of the scan points to explicitly lie outside the decoupling limit in order to cover interesting effects of large mixing in the Higgs sector. However, we found that almost all of these points are rejected by MSSM Higgs searches, most strongly by recent LHC searches. We placed special emphasis on interpreting the current LHC limits for heavy stable charged particles. Data from the 7 and 8 TeV LHC runs were taken into account. Further, we explicitly included the possibility of long-lived colored sparticles appearing due to phase space suppression. We found that long-lived gluinos can appear for mass gaps up to ∆m 300 GeV if all 2-and 3-body decays are kinematically forbidden. Accordingly, we included the R-hadron searches performed by CMS in our analysis. The obtained results imply conservative mass limits on some of the model parameters. These limits most importantly constrain the yield in co-annihilation regions. Furthermore, we showed the effects of the constraints from collider searches for MSSM Higgs signals, from flavor and precision observables as well as from CCB bounds on the allowed values of the stau yield in different regions.
We found that all points with stau yields Y 10 −14 that feature a dominant annihilation into Higgs final states were excluded by these bounds. Points with Y < 10 −14 only survived in the vicinity of the resonant pole of the Higgs propagator at m A ≃ 2m τ 1 . However, we encountered different scenarios with this feature. For staus with a large left-right mixing their annihilation via an s-channel heavy Higgs provides the most effective way to achieve low stau yields, which can reach roughly 2 × 10 −16 . For cases without mixing in the stau sector, we found two other possibilities to obtain small stau yields: co-annihilation with EWinos with a significant higgsino admixture as well as coannihilation with stops or sbottoms with considerable left-right mixing-in both cases annihilation near the resonant pole of an s-channel Higgs is required. We found allowed points down to Y ≃ 5 × 10 −15 and Y ≃ 10 −14 in the former and latter case, respectively.
Thus, our results show that the current constraints on the parameter space of the MSSM with a long-lived stau NLSP still allow for a stau relic abundance small enough to satisfy the strict bounds from big bang nucleosynthesis. The smallness of the corre-sponding region in parameter space suggests distinct features that will be probed in the upcoming LHC run.
The thermally averaged annihilation cross section is defined by where the sum runs over all supersymmetric initial state particles i, j. Further, n eq i,j and n eq are the individual and total equilibrium number densities, respectively. The thermal average σ ij v ij is given by where p i and f i are the three-momentum and the equilibrium phase-space density of particle i, respectively. Further, v ij is the Møller velocity, defined by Yield Y and density fraction Ω The relation between the yield and the density fraction Ω of a relic particle is where ρ 0 would be the current density of the relic if it had not decayed, ρ c is the critical density, and s 0 is the current entropy density of the universe. Inserting the numerical values [93] yields Y = 3.747 × 10 −9 Ωh 2 GeV m .
This expression is used to compute the yield from the output of micrOMEGAs.

B Mixing in the stau sector
Considering real parameters, we denote the stau mass matrix by where with T 3 and Q referring to the weak isospin and the electric charge, respectively. The stau mixing matrix reads R τ = cos θ τ sin θ τ − sin θ τ cos θ τ .

C Sfermion-sfermion-Higgs couplings
In the MSSM, the couplings of the lighter mass eigenstates of the third generation sfermions, τ 1 , b 1 and t 1 , to the CP -even neutral Higgses h and H are given by

D Resummation of the Higgs-bottom couplings
The tree-level Higgs-bottom couplings in the MSSM read (see, e.g., [95]) Radiative corrections to these couplings can be significant [96,97,98,99]. For positive µ and A t , they typically lead to a suppression of the couplings. The leading tan βenhanced terms can be resummed to all orders in perturbation theory [100,101] leading to the approximate relative corrections [102] h hbb h tree h Hbb h tree The leading contributions to ∆(m b ) come from the gluino-sbottom loop and from the charged higgsino-stop, wino-stop and wino-sbottom loops and are given by [100] ∆(m b ) ≃ 2α s 3π m g µ tan β I(m b 1 , m b 2 , m g ) + h 2 t 16π 2 µA t tan β I(m t 1 , m t 2 , µ) where I(a, b, c) = 1 (a 2 − b 2 )(b 2 − c 2 )(a 2 − c 2 ) a 2 b 2 log a 2 b 2 + b 2 c 2 log b 2 c 2 + c 2 a 2 log c 2 a 2 . (85) Note that in the decoupling limit, α ≃ β − π/2, so the hbb coupling remains SM-like even in the presence of large values for ∆(m b ) since tan α tan β ≃ −1. Thus, the correction vanishes. The Hbb coupling reads in the decoupling limit.

E Ranges for the dedicated scans
In table 2 we list all dedicated scan regions for the 17-dimensional pMSSM scan introduced in section 4. The dedicated scan regions are motivated by the results of section 3.
In each region, the table displays those parameters that are constrained to a smaller range than given in   Table 2: Summary of all scan regions, the corresponding percentage of points, and the parameters whose scan ranges deviate from the ones given in table 1. All parameters not listed are scanned over according to table 1. We generated a total amount of 5 × 10 5 points in the 17-dimensional parameter space.