Coscattering in next-to-minimal dark matter and split supersymmetry

In some models of thermal relic dark matter, the relic abundance may be set by inelastic scattering processes (rather than annihilations) becoming inefficient as the universe cools down. This effect has been called coscattering. We present a procedure to numerically solve the full momentum-dependent Boltzmann equations in coscattering, which allows for a precise calculation of the dark matter relic density including the effects of early kinetic decoupling. We apply our method to a simple model, containing a fermionic SU(2) triplet and a fermionic singlet with electroweak-scale masses, at small triplet-singlet mixing. The relic density can be set by either coannihilation or, at values of the mixing angle θ . 10−5, by coscattering. We identify the parameter ranges which give rise to the observed relic abundance. As a special case, we study bino-like dark matter in split supersymmetry at large μ.


Introduction
Weakly interacting massive particles (WIMPs) with electroweak-scale masses are among the bestmotivated candidates for particle dark matter. In WIMP models, the dark matter particle is in equilibrium with the thermal bath of Standard Model particles at early cosmic times, until the rates of WIMP-number changing processes drop below the Hubble expansion rate. Below the temperature of this so-called freeze-out, the dark matter number density becomes effectively constant.
It has long been known that the presence of additional states ψ whose masses are close to the dark matter mass can significantly affect the prediction for the dark matter relic density. This is the case when these states are able to coannihilate with the dark matter particle χ [1]. It is also well known, but perhaps less universally appreciated, that coannhilations can lead to the observed relic density even when the actual χχ → X and χψ → X annihilation cross-sections are several orders of magnitude below the typical electroweak cross-section. All that is needed is for the coannihilation partners ψ to annihilate efficiently among themselves, and for the dark matter particle to remain in equilibrium via inelastic scattering χX → ψX , such that any χ overdensity is rapidly converted into a ψ overdensity which is subsequently washed out.
It is interesting to study scenarios where such inelastic scattering processes start becoming inefficient before ψψ → X annihilations does [2,3]. In that case, the usual coannihilation formalism [1,4,5] fails, since one of its core assumptions is that χ is always in equilibrium with ψ. Having the dark matter relic abundance dictated by the freeze-out of inelastic scattering processes rather than of coannihilations has been dubbed "coscattering" in [2]. This mechanism has since been explored in the context of several different models [6][7][8].
In this paper we will study coscattering in one of the simplest models where it can occur, namely, the singlet-triplet next-to-minimal dark matter model of [9,10]. A special case of this model is split supersymmetry [11,12] with a bino-like lightest neutralino and a large higgsino mass parameter. Besides the Standard Model particles, the model contains only two multiplets with electroweakscale masses, namely, a fermionic singlet χ and a fermionic SU(2) triplet ψ. We will give Majorana masses of the order of the electroweak scale to both these fields, with χ being slightly lighter than ψ, and impose a Z 2 symmetry under which they are odd. This latter symmetry forbids any renormalizable interactions between χ and the Standard Model states, but allows for the dimension-5 operator 1 Λ χψH † H which mixes the neutral components of χ and ψ (here H is the Standard Model Higgs doublet and Λ is a cutoff scale). At large Λ, or equivalently small mixing angles, χχ → X annihilations as well as χψ → X coannihilations quickly become negligible as the temperature drops below the dark matter mass. However, ψψ → X annihilations as well as χX → ψX scattering remain efficient at first, the former because they are not suppressed by the mixing angle and the latter because of its less severe Boltzmann suppression (taking X and X to be relativistic Standard Model particles). Depending on the masses and the mixing angle, either ψψ annihilation or inelastic scattering may be the first process to start decoupling as the temperature drops further, thus giving rise to either coannihilation or coscattering.
A subtlety of the coscattering phase is early kinetic decoupling, as local (or "kinetic") equilibrium is lost through the decoupling of the very same inelastic scattering processes which determine the dark matter relic abundance. Therefore, another key assumption of the usual coannihilation formalism is violated in coscattering scenarios, since the momentum distribution of dark matter at freeze-out is not necessarily an equilibrium (Maxwell-Boltzmann) distribution. This needs to be taken into account properly when predicting the thermal relic density, as emphasized already in [2,3], and may significantly change the result in certain cases [13,14].
The aim of this paper is to establish a framework allowing for an accurate computation of the dark matter relic abundance in coscattering models, including the effects of early kinetic decoupling, building on the formalism proposed in [3]. We will apply our method to analyze the singlet-triplet model at and beyond the transition between coannihilation and coscattering.
After briefly recalling the essential properties of singlet-triplet next-to-minimal dark matter, we will proceed to review and to further develop the formalism for calculating the dark matter relic density. We will then present some numerical results in the singlet-triplet model, and as a special case, we will discuss the parameter space for bino-like dark matter in split supersymmetry. We will conclude with some remarks about present constraints and possible future experimental signatures relevant to these models, and with a brief summary.
2 Singlet-triplet next-to-minimal dark matter We will proceed with a brief review of the singlet-triplet model studied in [9,10], emphasizing its possible origins in split supersymmetry. We add to the Standard Model a fermionic singlet χ as well as its fermionic coannihilation partner ψ with quantum numbers (1, 3) 0 under SU(3) × SU(2) × U(1). The most general Lagrangian compatible with this particle content and a Z 2 symmetry under which χ and ψ are odd is where L 5 contains the dimension-5 operators Here H is the Standard Model Higgs doublet. Apart from the Standard Model Weinberg operator these are the only possible dimension-5 terms. We will assume that dimension-6 and higher operators are negligible in the following. We will restrict our study to real parameters and choose M > 0 without loss of generality. Since the dark matter particle is to be χ-like, we assume that M > |m|.
This model can be UV-completed with a Z 2 -odd Dirac fermion doublet Ψ with hypercharge 1 2 whose mass is of the order of Λ. Identifying χ ∼ B, ψ ∼ W , Ψ ∼ H † u , H d , the particle content is then exactly the low-energy particle content of split supersymmetry, except for a gluino which will play no role in our considerations (assuming that it is heavy enough not to conflict with collider bounds). 1 In split supersymmetry, the cutoff scale Λ is identified with the Higgsino mass scale |µ|, the singlet and triplet masses m and M with the supersymmetry-breaking gaugino mass parameters M 1 and M 2 , and the Wilson coefficients λ, κ and κ in Eq. (2) are given at the renormalization scale |µ| by the matching conditions λ = g ug d +g dg u sign µ , κ =g ugd sign µ , κ =g ug d sign µ .
Hereg u,d andg u,d are the usual gaugino-Higgsino-Higgs Yukawa couplings of split supersymmetry [12]. It is understood that the split supersymmetry model will be embedded in a fully supersymmetric model at an even higher scale M SUSY which is the mass scale of the remaining superpartners. Hence we are assuming the mass hierarchy Since we are not assuming any particular model of supersymmetry breaking, this is a possible choice of parameters. In fact, the origins of the µ parameter being supersymmetric, there is no a priori reason for it to be correlated with any of the supersymmetry-breaking parameters (up to the issue of radiative corrections, which we will comment on in Section 4).
Regardless of the nature of the UV completion, after electroweak symmetry breaking and upon replacing the Higgs field by its vacuum expectation value the first two terms in Eq. (2) will induce a shift in the effective χ and ψ mass parameters respectively, which we will absorb in the definitions of m and M . The third term will cause the neutral component of ψ to mix with χ. The mixing angle is We will study O(1) Wilson coefficients λ and mass differences M − m which are not parametrically smaller than the electroweak scale, such that the mixing angle is roughly of the order v/Λ. The mass eigenstates in the dark matter sector are then a neutral Majorana fermion χ 0 which is mostly χ-like, a neutral Majorana fermion ψ 0 which is mostly ψ-like, and a purely ψ-like charged Dirac fermion ψ ± . (In split supersymmetry, these would correspond to a bino-like lightest neutralino χ 0 1 , a wino-like next-to-lightest neutralino χ 0 2 and a wino-like chargino χ ± 1 respectively, with the higgsino-like states at the higher mass scale |µ|). In the limit of sending the cutoff scale to infinity, χ 0 becomes purely χ-like and completely inert. The mass degeneracy between ψ 0 and ψ ± is lifted by ψ 0 − χ 0 mixing and by electroweak loops, the latter inducing a mass difference of around 160 MeV between the neutral and the charged ψ-like states [15]. The mixing-induced mass difference is even smaller for the mixing angles of interest to us.

Coannihilation and coscattering in the singlet-triplet model
The fermionic singlet-triplet model described in Section 2 is a prime example of a model which can produce the dark matter relic density via either coannihilation or coscattering. While coannihilation in the singlet-triplet (or wino-bino) model has been extensively discussed in the literature (see e.g. [9,[16][17][18][19][20][21]), the coscattering phase has not been studied previously.
We take M and m of the order of a few 100 GeV and degenerate within < 10%. The effective dark matter number density will be depleted by χχ → X, χψ → X and ψψ → X annihilation processes (where ψ is any of the triplet-like states ψ 0 , ψ + or ψ − , and here and in the following we write χ instead of χ 0 to simplify notation -it should be understood that we are dealing with a mass eigenstate which contains a small admixture of ψ 0 ). As the cutoff scale Λ becomes large and the mixing angle becomes small, θ 1, the annihilation cross-sections involving χ become subdominant, while ψψ → X annihilations are independent of Λ. The dark matter particle χ is kept in chemical equilibrium as long as χX → ψX scattering remains efficient; these processes are also suppressed by Λ but, in contrast to annihilations, they are only singly-Boltzmann suppressed if X and X are relativistic Standard Model states.
In the coannihilation phase of the model, the mixing angle is large enough to keep χ in equilibrium until ψψ → X annihilations freeze out. For this phase, the resulting dark matter relic density can be reliably computed using standard public codes. By contrast, in the coscattering phase, the mixing angle is so small that the χX → ψX scattering rate drops below the Hubble rate before the ψψ → X annihilation rate does. The χ number density will then decrease much more slowly and eventually become constant, while the ψ-like states will continue to annihilate until the freeze-out temperature T f ≈ M/25 is reached. At even lower temperatures, the few remaining ψ-like states will ultimately slowly convert or decay into χ, thus increasing the final dark matter relic abundance by a small amount.
We proceed by presenting the details of the calculation of the thermal relic abundance, reviewing and extending the formalism developed in [3] as we go along. Our starting point is the Boltzmann equation in a Friedmann-Robertson-Walker universe, Here t is cosmic time and H is the Hubble parameter, which, in a radiation-dominated universe, is given as a function of temperature by with g * (T ) the effective number of degrees of freedom and M P = 2.435 × 10 18 GeV the reduced Planck mass. Moreover, f χ is the dark matter phase space distribution which depends on the 3-momentum p χ (more precisely, it depends only on its modulus p χ by isotropy), E χ is the corresponding energy, and C[f χ ] is the collision operator governing the interactions between χ and the other particle species.
In the coannihilation phase at small mixing angle, as well as in the coscattering phase, the only processes contributing sizeably to C[f χ ] are inelastic scatterings off the thermal bath, where X and X are Standard Model particles and ψ can be either ψ 0 , ψ + or ψ − . Therefore, the collision operator C[f χ ] is given by with M the matrix element for the process (9), suitably summed and averaged over final and initial internal degrees of freedom, and dp i the Lorentz-invariant phase space measure We can safely assume that X and X are in chemical and kinetic equilibrium throughout. Moreover, we can assume that at least kinetic equilibrium is maintained for ψ because of efficient elastic and inelastic scattering processes with the thermal bath. However, we allow for the integrated number density to deviate from its equilibrium value, which in the Maxwell-Boltzmann approximation (f eq ψ = e −E ψ /T , neglecting quantum statistical factors) is given by Hence for ψ we use the standard ansatz where n ψ (t) is an unknown function to be determined. Detailed balance allows to replace, in Eq. (10), so that the collision term becomes Here we have again neglected quantum statistical factors, setting f eq X = e −E X /T . The function w(s) [22] is related to the unpolarized cross-section σ(s) by and can be obtained from the matrix element by integrating over the final state solid angle: By changing integration variables from cos θ X to s = m 2 χ + m 2 X + 2E χ E X − 2p χ p X cos θ X , and carrying out the φ X and p X integrations, Eq. (16) becomes where The dominant contribution to scattering in the singlet-triplet model is given by W -mediated processes where X and X are relativistic Standard Model fermions. 2 Neglecting Standard Model fermion masses one obtains where g 2 is the SU(2) gauge coupling and θ is the singlet-triplet mixing angle. Summing over the possible channels, i.e. setting ψ = ψ + and X = + , ν , u,d, c,s (where is any lepton) or the respective charge-conjugate states, gives an overall factor of 36 in the final expression for the collision term. The contribution of third-generation quarks will be Boltzmann-suppressed and therefore subdominant, since the top quark is non-relativistic at the temperatures of interest, provided that m is below the TeV scale. Other subdominant contributions come from processes where X and X are electroweak gauge bosons and Higgs bosons, as detailed in the Appendix.
We note that, if χ were guaranteed to be in kinetic equilibrium such that (7) could be further simplified by integrating over p χ on both sides, thus obtaining an ordinary differential equation for the integrated number density n χ (t), as in the standard formalism in usual coannihilation scenarios. However, in the coscattering phase of our model, kinetic equilibrium is lost at the same time as chemical equilibrium (and due to the decoupling of the same process), hence we need to solve Eq. (7) for all momentum modes separately, even if the final quantity we are interested in is the overall χ abundance.
To solve Eq. (7), we integrate over the solid angle of p χ and change variables from (t, p χ ) to In terms of these variables, where we have defined The Boltzmann equation Eq. (7) becomes where On each of the characteristic curves Eq. (24) becomes an ordinary differential equation Here the family of characteristics is parameterized by ρ such that q ρ (x 0 ) = ρ. The initial conditions for these ODEs are obtained by imposing that at early times, i.e. at some x 0 ≈ 1, χ is in equilibrium: The solution of Eq. (27) subject to these initial condition is f eq χ (y, q ρ (y)) B(y, q ρ (y)) exp − x y B(z, q ρ (z)) dz dy .
The right-hand side of Eq. (29) depends on the ψ + abundance n ψ (x), which is itself given by the solution of a Boltzmann equation whose collision operator depends on f χ . At sufficiently large temperatures, or equivalently at small x, one may set n ψ (x) ≈ n eq ψ (x) and use Eq. (29) to find the solution for f χ (x). In general, this is no longer true at larger x, so the coupled system of Boltzmann equations for the χ modes as well as for the momentum-integrated ψ number densities needs to be solved.
The Boltzmann equations governing the ψ abundances are most conveniently expressed in terms of the dimensionless quantities is the entropy density. 3 Following the steps in the standard coannihilation formalism [5] one obtains Here we have set Y eq ψ + = Y eq ψ 0 ≡ Y eq ψ , neglecting the small difference between m ψ 0 and m ψ + . It is understood that Y ψ − = Y ψ + , and the factors of 2 in the fourth and fifth line of Eqns. (32) are to account for terms where ψ − appears instead of ψ + . The thermally averaged annihilation cross-sections are defined as usual: Finally, σv ψ 0 X → ψ + X conversions are always efficient, which allows to work with a single effective ψ abundance Y ψ = Y ψ + = Y ψ 0 , similarly as in the standard coannihilation formalism. Summing the Boltzmann equations for Y ψ + , Y ψ − and Y ψ 0 , the ψ → ψ conversion terms cancel and one obtains Eqns. (27) and (35) constitute a coupled system of ordinary differential equations which we need to solve, subject to the initial conditions at some suitably small x 0 For a numerical solution of this system, we discretize q and work with N = 30 characteristic curves. The initial values of q at x = x 0 = 1 are chosen to lie between ρ = 0 and ρ = 50, at suitable points for efficient computation of the q-integral in Eq. (35) via N -point Gauss-Legendre quadrature. Since δ g is numerically small, the range of q covered by these curves diminishes by only about 10% between x = 1 and x = 400, see Eq. (26). We extract the annihilation cross-sections for the ψ-like states from micrOMEGAs [23]. Using CalcHEP [24] to compute the subdominant contributions to the χX → ψX scattering cross-sections (i.e. all contributions except those involving light fermions, for which we have the analytic expression Eq. (21)), we find that including these corrections changes the results at the sub-percent level for M = 500 GeV and by < 2% for M = 1000 GeV. To speed up the computation, we impose that ψ remains in equilibrium between x 0 = 1 and an intermediate temperature x int = 15, which allows to calculate f χ (x int , q ρ (x int )) directly from Eq. (29). For x > x int the full system of 31 differential equations is then solved using a standard adaptive fourth-order Runge-Kutta method. Fig. 1 shows the evolution of Y χ and Y ψ as a function of x for two representative points, one in the coannihilation phase and the other in the coscattering phase. In the coannihilation phase, the χ number density evolves along with the ψ number density, departing from equilibrium around x = x f ≈ 25, where ψψ → XX annihilations freeze out. At larger x it becomes approximately constant, while Y ψ continues to decline due to ψX → χX conversion which is still efficient (and energetically favored over the reverse process). In the coscattering phase, the dark matter abundance departs from its equilibrium value at relatively small x, whereas the ψ abundance remains in equilibrium until x = x f . Its slow decline at larger x is due to the now marginally efficient conversion processes. Eventually all remaining ψ-like states will decay into χ, but since these decays are three-body phase-space suppressed and therefore happen on a much larger timescale, we do not take them into account here. (The final relic abundance is nevertheless calculated from the total number density As is evident from the right panel of Fig. 1, the freeze-out of the scattering processes happens gradually, the deviation from the equilibrium number density being noticeable long before the number density ultimately settles and becomes constant. This is in marked contrast to the freeze-out of ψ annihilations which happens relatively abruptly. The left panel of Fig. 2 shows the shape of the momentum distributions q 2 f χ (q) at several x in the coscattering phase. One observes a marked departure from kinetic equilibrium starting This can also be seen from the right panel, which shows the rescaled collision operator B of Eq. (24) for the same data point as a function of q; one observes that the effective collision rate increases with q. When the momentum distribution is finally frozen (which is the case at x 100 for this choice of parameters), it bears little resemblance to a Maxwell-Boltzmann equilibrium distribution, leaning noticeably towards the higher modes to the right of the maximum. Consequently, the final relic abundance cannot be calculated reliably with the momentum-integrated Boltzmann equations by setting f χ (q, x) ∝ f eq χ (q, x). In the coscattering phase, using this simplification may change the predicted result for the relic density by factor of a few, and its use is justified only in the coannihilation phase. This is illustrated in Fig. 3, where we plot Y χ along with the prediction for Y χ obtained by two approximations. The first approximation (the dotted curve) is to impose f χ (q, x) ∝ f eq χ (q, x) and to solve a momentum-integrated Boltzmann equation not just for ψ but also for χ. The second approximation (the dashed curve) is to solve the momentum-dependent Boltzmann equation but to impose n ψ (x) = n eq ψ (x) for all x. Evidently, the use of neither of these approximations is justified in the coscattering phase. We have verified, however, that the χ abundance can be reliably calculated using the momentum-integrated Boltzmann equations in the coannihilation phase, as expected.
Thus, in the singlet-triplet model, early kinetic decoupling leaves a significant imprint on the final relic abundance. This is in contrast to the earlier study [3], which found numerically similar results for the relic densities predicted with either the full Boltzmann equations or the integrated ones, with a difference of the order of only 10%. However, the model considered in that study differs from ours in one crucial aspect, namely, the importance of decays and inverse decays ψ ↔ Xχ. In our model, these processes are phase-space suppressed because the mass difference M − m is always below m W , whereas in the model of [3] the mediator particle (the analogue of our ψ) can decay into χ and an on-shell b quark, rendering decays and inverse decays efficient until rather late times and thus helping to re-establish approximate kinetic equilibrium.
We illustrate the transition between the coannihilation and the coscattering phase in the singlet- In the right panel, the ratio Y χ /Y approx χ between the actual χ abundance and the abundance obtained with these two approximations.
triplet model in Fig. 4, where the curve indicates the parameters for which the observed dark matter relic density Ωh 2 = 0.120 [25] is reproduced. Here we are assuming m > 0, hence we can identify m χ = m and m ψ ± = m ψ 0 = M . In the coannihilation phase, the predicted relic density is almost independent of the mixing angle but depends sensitively on the mass difference M − m, whereas in the coscattering phase, it is the mixing angle which determines at what time the dark matter density departs from its equilibrium value.

Split supersymmetry with a bino LSP
The results of the previous section are valid for a general singlet-triplet model. However, when supposing that the UV completion is split supersymmetry, we can use them to obtain information on the value of the µ parameter. This is because the split supersymmetry couplingsg u,d andg u,d are given by the electroweak gauge couplings and by tan β at the scale M SUSY , Moreover, tan β and M SUSY are related by the condition λ H = 1 4 where λ H is the Standard Model Higgs quartic coupling. The couplingsg u,d andg u,d determine the Wilson coefficient λ of Eq.
Between |µ| and the electroweak scale, λ evolves according to its renormalization group equation (RGE), which at the one-loop level reads (39) .
Therefore, given either M SUSY or tan β, the mixing angle θ is fully determined by the bino, wino and higgsino masses ±m = m χ 0 1 , M = m χ 0 2 and |µ| = m χ 0 3,4 . Fixing tan β, we can now pick any point in the space of θ, M and m which gives rise to the observed relic density. We can then determine the corresponding µ by iteratively solving the RGEs for λ and the Standard Model couplings between the electroweak scale and |µ|, and the split supersymmetry RGEs [12] between |µ| and M SUSY , using the matching conditions Eqns. The result for tan β = 1 is shown in the right panel of Fig. 4. To obtain this curve we have used two-loop RGEs for the gauge couplings and one-loop RGEs with some dominant two-loop corrections for the others, with tree-level matching. While the precision of this calculation could be improved with state-of-the-art tools, we find that it is sufficient to estimate µ in this scenario, given that the RG evolution of the relevant couplings tends to be rather slow (a naive "tree-level" estimation of µ, using electroweak-scale values for all couplings, is within 5% of our loop-corrected result). For this choice of parameters, the scale M SUSY is given by the scale where the Standard Model Higgs quartic coupling λ H crosses zero. One finds that M SUSY tends to be about two orders of magnitude above µ, but since the RG evolution of λ H is subject to large uncertainties due to the uncertainty on the top Yukawa coupling, this result should be taken with a grain of salt.
It is interesting to note that a bino-like LSP coannihilating or coscattering with a wino-like neutralino and chargino remains an open and testable option for neutralino dark matter in split supersymmetry, besides the much-studied almost-pure bino and almost-pure higgsino cases which are not testable at present colliders (we remark that a sub-TeV bino-like LSP coannihilating with a higgsino is ruled out by now, see e.g. [26]). Indeed, our analysis shows that sub-TeV neutralino dark matter remains perfectly viable and necessarily comes with a coannihilation partner within the kinematic reach of the LHC, see Section 5.
A few comments about split supersymmetry with a very large µ term are finally in order. First, it is debatable to what extent the hierarchy |µ| |M 1,2 | is natural and whether or not that constitutes a problem. One might argue that the µ parameter is allowed by supersymmetry, so whatever dynamics generates the higgsino masses need not involve supersymmetry breaking, and the mass scales |µ| and |M i | need not be correlated. However, it is well known that in split supersymmetry the gaugino masses and µ are no longer separately protected against additive renormalization. In fact, a large µ parameter induces large one-loop threshold corrections to the gaugino masses at the higgsino scale: Here m A ≈ M SUSY is the mass scale of the heavy Higgs bosons. Hence, a hierarchy |µ|/|M 1,2 | 100 at moderate tan β requires some fine-tuning of the model parameters, since large radiative corrections to the gaugino masses must be cancelled against large tree-level values. We will not attempt to quantify the required fine-tuning here given that split supersymmetry is unnatural by construction anyway (not to mention that the dark matter relic density is also rather sensitive to the precise value of the mass difference M − m).
Second, note that the usual upper bound on the UV completion scale of split supersymmetry from vacuum stability [27] does not apply to the case where higgsinos are heavy. This is because, in the absence of dynamical higgsinos, the renormalization group running of the Higgs quartic coupling λ H is not modified at the one-loop level with respect to the Standard Model. Therefore this coupling can remain positive up to high scales, as is confirmed by our numerical analysis. In particular, a mass hierarchy M SUSY |µ| with |µ| ≈ 10 8 GeV is not excluded by vacuum stability.
Lastly, the gauge couplings in split supersymmetry with heavy higgsinos will not unify unless there are additional heavy states in incomplete GUT representations in the spectrum.

Constraints and signatures
As detailed e.g. in [9], next-to-minimal dark matter at small mixing angles evades the constraints from both direct and indirect dark matter detection experiments, since the dark matter particle is essentially singlet-like.
Big bang nucleosynthesis will place constraints on the singlet-triplet model if χ and ψ are almost mass degenerate, since for sufficiently small mass splittings and mixing angles, the ψ 0 → χ decay width will be highly suppressed and the ψ 0 lifetime can attain the range of seconds to minutes or more. Depending on the ψ 0 abundance after all scattering processes have frozen out, these late decays can conflict with the successful prediction for nuclear abundances. For example, taking θ = 1.22 × 10 −6 , M = 500 GeV and m = 490 GeV, the ψ 0 lifetime is about 10 s, which is still allowed by BBN constraints [28] given the predicted number density (see Fig. 2). However, for lifetimes 10 2 s, the constraints on the abundance become much more severe and the model will likely be excluded. A detailed analysis of the BBN constraints is beyond the scope of this work and left for a future study.
The possible collider signatures for singlet-triplet models have been analyzed in detail e.g. in [10,16,20,[29][30][31][32][33]. Therefore, here we will only briefly review the state of the most important LHC search for the model in the coscattering phase, which is the search for disappearing charged tracks. For this signature the dark matter particle χ plays practically no role, but it is the coannihilation partner ψ which is being probed. In detail, the charged states ψ ± are produced by the Drell-Yan process. For sufficiently small mixing angles θ 10 −4 they will preferentially decay into ψ 0 rather than into χ; ψ 0 is then stable on collider scales. The mass splitting m ψ + − m ψ 0 is induced by electroweak loops and is of the order of 160 MeV, which implies that ψ ± has a macroscopic lifetime because its decays are phase-space suppressed. It also implies that the only possible two-body decay is ψ ± → ψ 0 π ± with a very soft pion. The signature is therefore a charged ionization track left by ψ ± , with a typical length in the cm -m range, which at some point seemingly disappears (since the neutral ψ 0 is invisible and the pion is too soft to be identified).
In this channel, present data excludes triplet mass parameters below about 500 GeV [34,35], hence the benchmark points we have been using in Section 3 are still viable, if by a small margin. At 3000 fb −1 the projected LHC exclusion ranges up to M ≈ 900 GeV [36].

Summary
We have presented a precise calculation of the dark matter relic abundance from thermal freezeout in a scenario where the usual coannihilation formalism cannot be applied, as the relic density is set by coscattering. Our example model is singlet-triplet next-to-minimal dark matter (which can be realized in split supersymmetry), but the formalism could be applied to other similar models with little modification. The particle abundances in the dark matter sector are governed by the full momentum-dependent Boltzmann equations, which we have cast into a form amenable to a numerical solution, properly taking into account the effects of early kinetic decoupling. In the singlet-triplet model, the observed dark matter density can be reproduced by coannihilations for mixing angles θ 10 −5 and by coscattering for smaller mixing angles. The mass difference between the triplet-like and singlet-like states becomes smaller as the mixing angle is reduced. Our analysis also allows to chart the parameter space for bino-like dark matter in split supersymmetry with a very large µ parameter, where we find that coscattering starts setting in at values of |µ| 10 7 GeV. While we have not undertaken a full phenomenological study of the coscattering phase, we anticipate that the most promising experimental channel to constrain it will be LHC searches for disappearing charged tracks.

Appendix: Contributions to the collision operator for χ
In the singlet-triplet model, the dark matter particle is coupled to the Standard Model via the dimension-5 operators in Eq. (2). After electroweak symmetry breaking, they induce effective hhχ 0 χ 0 , hhχ 0 ψ 0 , hχ 0 χ 0 and hψ 0 χ 0 couplings (where h is the Higgs boson), and the last term gives rise to χ − ψ mixing and thus a χ 0 ψ ± W µ vertex. These are the only couplings between χ 0 and the Standard Model at the order 1/Λ (note that the χ 0 ψ 0 Z µ vertex is induced by a dimension-6 operator), and therefore the only ones relevant for the collision operator C[f χ ] if Λ is large.
Among the states in the spectrum, we have a mild mass hierarchy between the particles of the dark matter sector (χ 0 , ψ 0 , ψ ± ) and the Standard Model particles with masses around the electroweak scale (t, h, Z, W ). All other Standard Model particles can be treated as effectively massless. At the temperatures relevant to freeze-out, the dark matter sector states are of course non-relativistic, with their equilibrium number densities suppressed by Boltzmann factors ∼ e −m/T . We can therefore neglect all processes of the type χ 0 χ 0 → X and ψ ± χ 0 → X because they are both doubly Boltzmann-suppressed and Λ-suppressed. We can also neglect decays and inverse decays in the dark matter sector, since the mass splitting M − m is below the W mass, and therefore all two-body decays are kinematically forbidden.
What remains are 2 → 2 processes where dark matter sector particles scatter off the thermal bath. The dominant one among these is W -mediated χ 0 f → ψ ± f inelastic scattering where f and f are light Standard Model fermions (excluding third-generation quarks), since these are highly relativistic at the relevant temperatures. Note that Higgs-mediated χ 0 f → ψ 0 f or χ 0 f → χ 0 f scattering is m f -suppressed.
All processes involving t, h, Z or W in the initial or final state give subdominant contributions which, however, become relatively more important as the dark matter mass increases. For dark matter masses around 500 GeV, all these particles are already non-relativistic at x ≈ 20, whereas for dark matter masses around 1000 GeV, they are only on the brink of becoming non-relativistic at the temperatures of interest. For M = 500 GeV we have checked that including the subdominant processes changes the final relic abundance by less than 1%, and we consequently neglect them. For M = 1000 GeV we include the numerically most important ones among these corrections, specifically χ 0 t → ψ + b, χ 0b → ψ +t , χ 0 t → ψ 0 t, χ 0 W + → ψ + Z, χ 0 Z → ψ + W − , χ 0 W + → ψ + γ, χ 0 γ → ψ + W −, χ 0 W + → ψ + h, χ 0 h → ψ + W − and χ 0 h → ψ 0 h as well as the respective charge-conjugate processes. We find that their combined impact on the result is still at the < 2% level.